Initiation and Breakdown of an Axisymmetric Hydraulic Fracture Transverse to a Horizontal Wellbore

Abstract

We investigate the initiation and early-stage propagation of an axi-symmetric hydraulic fracture from a wellbore drilled in the direction of the minimum principal stress in an elastic and impermeable formation. Such a configuration is akin to the case of a horizontal well and a hydraulic fracture transverse to the well axis in an open hole completion. In addition to the effect of the wellbore on the elasticity equation, the effect of the injection system compressibility is also taken into account. The formulation accounts for the strong coupling between the elasticity equation, the flow of the injected fluid within the newly created crack and the fracture propagation condition. Dimensional analysis of the problem reveals that three dimensionless parameters control the entire problem: the ratio of the initial defect length over the wellbore radius, the ratio between the wellbore radius and a length-scale associated with the fluid stored by compressibility in the injection system during the well pressurization, and finally the ratio of the time-scale of transition from viscosity to toughness dominated propagation to the time-scale associated with compressibility effects. A fully coupled numerical solver is presented, and validated against solutions for a radial hydraulic fracture propagating in an infinite medium. The influence of the different parameters on the transition from the near-wellbore to the case of a hydraulic fracture propagating in an infinite medium is fully discussed.

1. Introduction

In this study, we are interested in the initiation of hydraulic fractures from an open-hole horizontal well. Horizontal wells are drilled preferably in the direction of the minimum horizontal stress in order to create hydraulic fractures perpendicular to the wellbore and therefore maximize the drainage from the reservoir. During the pressurization of the wellbore, tangential tensile stress is generated which can result in the initiation of longitudinal fractures parallel to the wellbore axis [1,2]. In an isotropic elastic medium, depending upon the stress field, the initiation of fractures transverse to the wellbore axis is favored by creating an initial flaw, i.e. an axisymmetric notch, of sufficient depth [3] (see Figure. 1 for a sketch). We focus solely on transverse hydraulic fractures in this contribution.

In the initial stage of propagation, these transverse fractures can be idealized as axisymmetric (radial) fractures around the wellbore until they hit a stress barrier or other type of heterogeneities. We investigate the initiation and propagation of such a fracture under constant injection of a Newtonian fluid, focusing on the case of “tight” rocks where leak-off is typically negligible. In particular, we are interested in clarifying when the effects associated with the near-wellbore region dissipate and no longer affect the hydraulic fracture propagation and how this transition takes place.

The initiation and the early stage of the propagation of such a hydraulic fracture is affected by two “near wellbore” effects: i) the finiteness of the wellbore and the initial flaw length and ii) a transient phenomenon associated with the release of the fluid stored by compressibility in the wellbore during the pressurization phase prior to the initiation of the fracture. This second effect is ultimately linked to the compressibility of the injection system (i.e. mostly the fluid volume stored by compressibility within the wellbore). These effects have been investigated for the case of plane-strain fractures in [4,5] and for the case of axisymmetric hydraulic fractures driven by an inviscid fluid (i.e. zero-viscosity) in [6]. In this paper, we extend these contributions to the case of viscous flow in axi-symmetric fractures.

A detailed dimensional analysis is performed indicating various time scales, length scales and dimensionless numbers controlling the problem. A numerical solver, along the lines of previous contributions [7,8], is presented. A series of numerical simulations are performed in order to study the transition from the near-wellbore propagation regime to the case of a radial hydraulic fracture in an infinite medium under constant injection rate.

2. Problem statement

Let us consider a horizontal well of radius $a$ drilled in the direction of the minimum horizontal stress ${\sigma }_h$. In this ideal configuration, we are interested in the initiation and propagation of a hydraulic fracture from a radial axisymmetric notch of length ${\mathcal{l}}_o$ transverse to the wellbore (Figure. 1). As previously mentioned, we assume that the radial notch length is sufficient to favor a transverse fracture as compared to a longitudinal one (see [3] for discussion on the effect of the stress field on the competition between transverse and longitudinal fracture initiation).

In the absence of any stress barriers or other heterogeneities, this axisymmetric transverse fracture of radius $R$ will transition toward a penny-shaped geometry when its radius is much greater than the wellbore radius (Figure. 1). We will denote as $p_f$ the fluid pressure inside the fracture (and in the wellbore), $p=p_f-{\sigma }_h$ the net pressure opening the fracture face and $w$ the corresponding fracture width. The axial stress variation azimuthally around the wellbore is given by [9]

$${\sigma }_a={\sigma }_h-2\nu ({\sigma }_H-{\sigma }_v)\frac{a^2}{r^2}\cos 2\theta$$E1

where $r$ and $\theta$ are the polar coordinates centered at the wellbore, ${\sigma }_H$ is the maximum horizontal stress, ${\sigma }_v$ is the vertical stress, $\nu$ is the poisson’s ratio and $a$ is the wellbore radius. The azimuthal average of the axial stress (i.e. integrating the above equation from $\theta =0$ to $\theta =2\pi$) reduces to ${\sigma }_h$. We will neglect the azimuthal variation of the axial stress close to the wellbore wall as a first approximation. The axial stress therefore reduces to the minimum horizontal stress. Under such an approximation, the model is truly axi-symmetric, i.e. independent of $\theta$.

Figure 1.

Sketch of a transverse fracture propagating from a horizontal wellbore (top), axisymmetric model (bottom)

2.1. Elasticity

The relation between the fracture width and the net loading acting on the crack is expressed by a hyper-singular boundary integral equation following the method of distributed dislocations [10]:

$$p\left(r\right)=p_f\left(r\right)-{\sigma }_h=\frac{E^{\mathrm{\text{'}}}}{2\pi }\underset{a}{\overset{a+\mathcal{l}}{\int }}J\left(r,r^{\mathrm{\text{'}}}\right)\frac{\partial w\left(r^{\mathrm{\text{'}}}\right)}{\partial r^{\mathrm{\text{'}}}}d{r}^{\mathrm{\text{'}}}$$E2

where the singular kernel $J\left(r,r^{\text{'}}\right)$ is given by [11]. It provides the normal stress due to an axisymmetric dislocation around a wellbore of radius $a$. Details of this kernel are recalled in Appendix A for completeness.

2.2. Fluid flow

We neglect the fluid compressibility compared to the fracture compliance as is usually done in modeling hydraulic fractures. The balance of mass then reduces to a strict volume balance:

$$\frac{\partial w}{\partial t}+\frac{1}{r}\frac{\partial }{\partial r}\left(rq\right)=0\mathrm{ }$$E3

Under the hypothesis of lubrication theory (low Reynolds number flow) for a Newtonian fluid, the fluid flux $q$ is given by Poiseuille law

$$q=-\frac{w^3}{{\mu }^{\mathrm{\text{'}}}}\frac{\partial p}{\partial r}$$E4

where ${\mu }^{\text{'}}=12\mu$ with $\mu$ the fluid viscosity.

2.3. Boundary conditions

During the propagation of a hydraulic fracture, a fluid lag may develop at the tip of the fracture [12,13]. Such a fluid lag is larger at early time during the propagation. In order to check whether the fluid lag should be taken into account, we can estimate the characteristic timescale associated with the disappearance of the fluid lag, which is equal to $t_{om}=\frac{E^{\text{'}2}{\mu }^{\text{'}}}{{\sigma }_o^3}$ (see [14,7]), where ${\sigma }_o$ denotes the far-field confining stress (here the minimum horizontal stress). This timescale is inversely proportional to ${\sigma }_o^3$ which means that for higher values of the confining stress this timescale is very small (typically the case in deep formations). For illustration, we use stress typical of an unconventional reservoir (e.g. Barnett shale) and typical rock parameters:

$$E=5.4\times {10}^6\mathrm{ }psi,\mathrm{ } \nu =0.21,\mathrm{ } {\sigma }_o={\sigma }_h=3390\mathrm{ }psi\mathrm{ }\left(\text{at}\text{the}\text{depth}\text{of}\mathrm{ }5000ft\right)\mathrm{ },{\mathrm{ }K}_{Ic}=1500\mathrm{ }psi/\sqrt{in}.\mathrm{ }$$E5

For the case of a slick water stimulation (viscosity of 1 cP), the characteristic time $t_{om}$ is equal to $0.0014 sec$. For a gel treatment (tangent viscosity of approximately 100 cP), this characteristic time still remains small, $t_{om}=0.14 sec$. The transient effect associated with the disappearance of the fluid lag can thus be ignored for the conditions typically encountered in slick water fracturing of deep horizontal wells (i.e. in the Barnett shale). In the remaining of this paper, we assume that the fluid front coincides with the fracture front.

The conditions at the tip of the fracture are thus:

$$q\left(r=\mathcal{l}+a\right)=0\mathrm{ } w\left(r=\mathcal{l}+a\right)=0.$$E6

The fluid flux entering the fracture is equal to the total fluid injection rate (for example the injection rate at the wellhead) minus the fluid volume stored in the well due to its compressibility (essentially the compressibility of the fluid inside the wellbore). Therefore, from the fluid mass conservation in the wellbore (between the injection point and the fracture inlet), one can write the following boundary condition at the inlet of the fracture:

$$q=\frac{1}{2\pi a}\left(Q_o-U\frac{\partial p_b}{\partial t}\right),\mathrm{ }$$E7

where $U$ ($barrels/psi$) is the injection system compressibility and $p_b$ denotes the wellbore pressure. Pressure continuity is ensured at the fracture in-let by the following condition

$$p_f\left(r=a\right)=p_b.$$E8

It is important to note that no friction pressure drop (e.g. perforation drop) is taken into account in this injection boundary condition.

2.4. Initiation and fracture propagation condition

Prior to any opening of the initial defect of length ${\mathcal{l}}_o$, the pressurization rate is uniform ($q=0$in Equation (7)). We thus start the simulation only when the fluid pressure in the wellbore and in the notch has reached the minimum horizontal stress, which is the pressure at which this initial defect starts to open: we will denote this start-up time $t_o$. For modeling purposes, at $t_o$, the initial condition will be taken as a vanishingly small net pressure (i.e. a fluid pressure just slightly above the minimum horizontal stress). Due to the continuous fluid injection considered here, the wellbore pressure keeps increasing. The fracture will initiate its growth once the stress intensity factor reaches its critical value $K_{Ic}$, i.e. the fracture toughness of the rock. Once fracture initiation has occurred, we assume that the fracture propagates under quasi-static equilibrium such that the stress intensity factor is always equal to its critical value. For a pure mode I fracture considered here, this condition can be expressed as an asymptote on the fracture opening near the tip:

$$w~\frac{K^{\mathrm{\text{'}}}}{E^{\mathrm{\text{'}}}}{\left(\left(\mathcal{l}+a\right)-r\right)}^{1/2}\mathrm{ } 1-\frac{r}{\left(\mathcal{l}+a\right)}\ll 1\mathrm{ }.$$E9

where $K^{\mathrm{\text{'}}}=\sqrt{32/\pi }{K}_{Ic}$

3. Scaling

Let us scale the variables involved in the problem in order to grasp the effects of various physical phenomena acting at various scales. We introduce a characteristic length scale $L_{\star }$, characteristic fracture width $w_{\star }$, characteristic pressure $p_{\star }$, characteristic fluid flux $q_{\star }$, and a characteristic time $t_{\star }$ and scale the variables as follow:

$$r=L_{\star }\rho ,\mathrm{ } \mathcal{l}=L_{\star }\gamma ,\mathrm{ } a=L_{\star }{\mathcal{G}}_a,\mathrm{ } p=p_{\star }\mathrm{\Pi },\mathrm{ } w=w_{\star }\mathrm{\Omega },\mathrm{ } q=q_{\star }\mathrm{\Psi },\mathrm{ } t=t_{\star }\tau ,$$E10

$\rho$, $\gamma$and ${\mathcal{G}}_a$ denotes respectively the dimensionless coordinate along the fracture, the dimensionless fracture length and the dimensionless wellbore radius. The dimensionless opening, net pressure and fluid flux are denoted as $\mathrm{\Omega }$, $\mathrm{\Pi }$ and $\mathrm{\Psi }$ respectively.

Using the above scaling, the governing equations are converted to dimensionless form where the different scales are yet to be defined:

• Elasticity operator

$$\mathrm{\Pi }=\frac{{\mathcal{G}}_e}{2\pi {\mathcal{G}}_a}\left(\underset{1}{\overset{1+\raisebox{1ex}{$\gamma $}\!\left/ \!\raisebox{-1ex}{${\mathcal{G}}_a$}\right.}{\int }}\mathrm{J}\left(\rho ,{\rho }^{\mathrm{\text{'}}}\right)\frac{\partial \mathrm{\Omega }\left({\rho }^{\mathrm{\text{'}}}\right)}{\partial {\rho }^{\mathrm{\text{'}}}}d{\rho }^{\mathrm{\text{'}}}\right)$$E11

where ${\mathcal{G}}_e=\frac{E^{\text{'}}{w}_{\star }}{p_{\star }{L}_{\star }}$ is a dimensionless group associated with the elasticity operator and ${\mathcal{G}}_a=\frac{a}{L_{\star }}$ is a dimensionless group associated with the effect of the wellbore radius.

• Continuity equation

$$\frac{\partial \mathrm{\Omega }}{\partial \tau }+{\mathcal{G}}_r\frac{1}{\rho }\frac{\partial }{\partial \rho }\left(\rho \mathrm{\Psi }\right)=0\mathrm{ }$$E12

where ${\mathcal{G}}_r=\frac{t_{\star }{q}_{\star }}{w_{\star }{L}_{\star }}$ is a dimensionless group associated with fluid conservation within the fracture.

• Poiseuille law for lubrication

$$\mathrm{\Psi }=-\frac{{\mathrm{\Omega }}^3}{{\mathcal{G}}_m}\frac{\partial {\mathrm{\Pi }}_f}{\partial \rho }$$E13

where ${\mathcal{G}}_m=\frac{{\mathrm{\mu }}^{\mathrm{\text{'}}}{q}_{\star }{L}_{\star }}{w_{\star }^3{p}_{\star }}$ is a dimensionless group associated with fluid viscosity.

• Inlet boundary condition

$$\mathrm{\Psi }=\frac{1}{2\pi }\left({{\mathcal{G}}_q\mathrm{\Psi }}_o-{\mathcal{G}}_U\frac{\partial {\mathrm{\Pi }}_b}{\partial \tau }\right)$$E14

where ${\mathcal{G}}_q=\frac{Q_o}{L_{\star }{q}_{\star }}$ is a dimensionless group associated with the effect of the injection rate, and ${\mathcal{G}}_U=\frac{U{p}_{\star }}{L_{\star }{q}_{\star }{t}_{\star }}$ is a dimensionless group associated with the injection system compressibility.

• The LEFM propagation condition reduces to

$$\mathrm{\Omega }~{\mathcal{G}}_k{\left(\left(\gamma +{\mathcal{G}}_a\right)-\rho \right)}^{1/2},\mathrm{ } \text{for}\left(\gamma +{\mathcal{G}}_a\right)-\rho \ll 1\mathrm{ }$$E15

where ${\mathcal{G}}_k=\frac{K^{\text{'}}\sqrt{L_{\star }}}{E^{\text{'}}{w}_{\star }}$ is a dimensionless group associated with the rock fracture toughness.

First, it is natural to set the dimensionless groups associated with the injection rate ${\mathcal{G}}_q$ (we inject fluid) and elasticity ${\mathcal{G}}_e$ to unity. In setting ${\mathcal{G}}_e$ to 1, we account for the fact that the crack opening is typically much smaller than the fracture length and that the net pressure is also much smaller than the rock elastic modulus. The dimensionless group associated with fluid conservation ${\mathcal{G}}_r$ is also set to unity to account for the fact that injected fluid volume remains in the fracture in the absence of leak-off.

In the case of the propagation of a radial hydraulic fracture in an infinite medium, energy dissipation is attributed to two competing mechanisms i.e. viscous forces associated with fluid flow within the crack and the creation of new fracture surfaces (i.e. fracture toughness) [15]. At the beginning of the propagation, i.e. for small fracture radius, viscous forces are the dominant dissipative process, and fracture toughness can be neglected in such a viscosity dominated regime of propagation. A self-similar solution has been obtained in [15] for that case, and will be denoted as the $M$-vertex solution. As time increases, fracture energy slowly takes over viscous forces as the main dissipative mechanism. Ultimately, at large time, the fracture propagates in the so-called toughness dominated regime of propagation, where viscosity can be neglected. Here again, an analytical solution exists [16], and will be referred as the $K$-vertex solution. In an infinite medium, the radial hydraulic fracture therefore transition from the viscosity ($M$) to the toughness ($K$) regime of propagation.

This picture is modified when accounting for near-wellbore effects. These effects will eventually dissipate for fracture length much larger than the wellbore radius. This transition from the near-wellbore to the infinite medium solution is of particular interest. The effect of the wellbore and of the system compressibility will affect the system response at early time, i.e. when the radius of the fracture is comparable to that of the wellbore and when the system compressibility still has an effect. At large time, the transient associated with the fracture breakdown and the release of the fluid stored by compressibility prior to the crack initiation will become insignificant: i.e. the fluid flux entering the crack will then be equal to the injected flow rate. The solution will thus behave as the infinite medium solution [15,16,17].

It is therefore interesting to introduce two different scaling. The first scaling relates to the case where the system compressibility and toughness are important, i.e. at early time or for small fractures. We will denote such a scaling as the Compressibility-Toughness scaling and denote it as$UK$. This scaling is based on the characteristic time of transition from the compressibility effects ($U$) to the infinite medium solution corresponding to toughness dissipation ($K$). The characteristic scales in that Compressibility-Toughness scaling are obtained as:

$$L_{uk}={\left(E^{\mathrm{\text{'}}}U\right)}^{1/3},\mathrm{ } p_{uk}=\frac{K^{\mathrm{\text{'}}}}{E^{\mathrm{\text{'}}1/6\mathrm{ }}{U}^{1/6\mathrm{ }}},\mathrm{ } w_{uk}=\frac{K^{\mathrm{\text{'}}}{U}^{1/6\mathrm{ }}}{E^{\mathrm{\text{'}}5/6\mathrm{ }}},\mathrm{ } q_{uk}=\frac{Q_o}{E^{\mathrm{\text{'}}1/3\mathrm{ }}{U}^{1/3\mathrm{ }}},t_{uk}=\frac{K^{\mathrm{\text{'}}}{U}^{\mathrm{\text{'}}5/6}}{E^{\mathrm{\text{'}}1/6}{Q}_o}\mathrm{ }$$E16

In particular, the dimensionless time in that scaling is$\tau =\frac{t}{t_{uk}}$. The corresponding dimensionless numbers are:

$${{\mathcal{G}}_m}_{uk}=\frac{E^{\mathrm{\text{'}}8/3}{Q}_o{\mu }^{\mathrm{\text{'}}}}{K^{\mathrm{\text{'}}4}{U}^{1/3\mathrm{ }}},\mathrm{ } {{\mathcal{G}}_q}_{uk}=1,\mathrm{ } {{\mathcal{G}}_U}_{uk}=1,{{\mathcal{G}}_k}_{uk}=1,{{\mathcal{G}}_a}_{uk}=\frac{a}{{\left(E^{\mathrm{\text{'}}}U\right)}^{1/3}}$$E17

The second scaling of interest corresponds to the case where the transient effects associated with the wellbore and injection compressibility have vanished: i.e. when the model reduces to the case of radial fracture propagating in an infinite medium – we will call this scaling the Viscosity-Toughness scaling and denote it as $MK$. This scaling is based on the characteristic time of transition from the viscosity dissipation ($M$) to the toughness dissipation ($K$). The corresponding characteristic scales are [15]:

$$L_{mk}=\frac{E^{\mathrm{\text{'}}3}{Q}_o{\mu }^{\mathrm{\text{'}}}}{K^{\mathrm{\text{'}}4}},\mathrm{ }{p}_{mk}=\frac{K^{\mathrm{\text{'}}3}}{E^{\mathrm{\text{'}}3/2\mathrm{ }}\sqrt{Q_o{\mu }^{\mathrm{\text{'}}}}},w_{mk}=\sqrt{\frac{E^{\mathrm{\text{'}}}{Q}_o{\mu }^{\mathrm{\text{'}}}}{K^{\mathrm{\text{'}}2}}},q_{mk}=\frac{K^{\mathrm{\text{'}}4}}{E^{\mathrm{\text{'}}3}{\mu }^{\mathrm{\text{'}}}},t_{mk}=\frac{E^{\mathrm{\text{'}}13/2}{Q}_o^{3/2}{\mu }^{\mathrm{\text{'}}5/2}}{K^{\mathrm{\text{'}}9}}\mathrm{ }$$E18

The dimensionless time is here defined as $\tilde{\tau }=\frac{t}{t_{mk}}$ and the different dimensionless numbers in that scaling are given as

$${{\mathcal{G}}_m}_{mk}=1,\mathrm{ } {{\mathcal{G}}_q}_{mk}=1,\mathrm{ } {{\mathcal{G}}_U}_{mk}=\frac{K^{\mathrm{\text{'}}12}U}{E^{\mathrm{\text{'}}8}{Q}_o^3{{\mu }^{\mathrm{\text{'}}}}^3},{{\mathcal{G}}_k}_{mk}=1,{{\mathcal{G}}_a}_{mk}=\frac{a{K}^{\mathrm{\text{'}}4}}{E^{\mathrm{\text{'}}3}{Q}_o{\mu }^{\mathrm{\text{'}}}}$$E19

In order to grasp the transition from the early time where the near-wellbore effects are important to the large-time solution of propagation in an infinite medium, we introduce $\chi$ as the ratio of the timescales associated with the previously defined scalings:

$$\chi =\frac{t_{mk}}{t_{uk}}=\frac{E^{\mathrm{\text{'}}20/3}{Q}_o^{5/2}{\mu }^{\mathrm{\text{'}}5/2}}{K^{\mathrm{\text{'}}10}{U}^{5/6\mathrm{ }}}.\mathrm{ }$$E20

Large values of $\chi$corresponds to high viscosity ${\mu }^{\text{'}}$ and/or low injection system compressibility $U$. In this case, the transition from the near wellbore solution to the infinite medium solution occurs prior to the transition from the viscosity ($M$) to the toughness ($K$) dominated regime of propagation of an infinite radial hydraulic fracture. Smaller values of $\chi$ correspond to lower viscosity/higher injection system compressibility. In this case the transition from near wellbore solution to the infinite medium solution occurs in the $K$ regime of the infinite medium solution.

Now introducing $\mathcal{A}={{\mathcal{G}}_a}_{uk}=\frac{a}{{\left(E^{\text{'}}U\right)}^{1/3}}$, the ratio between the wellbore radius and a lengthscale associated with the volume stored in the injection system by compressibility, the dimensionless numbers for the two scalings previously presented can be written as:

$${{\mathcal{G}}_m}_{uk}={\mathrm{\chi }}^{2/5},\mathrm{ } {{\mathcal{G}}_a}_{uk}=\mathcal{A}$$E21

and

$${{\mathcal{G}}_U}_{mk}={\mathrm{\chi }}^{-6/5},{{\mathcal{G}}_a}_{mk}=\mathcal{A}{\chi }^{-2/5}$$E22

Correspondence between the two scalings can also be obtained as the function of $\chi$:

$$\frac{L_{mk}}{L_{uk}}={\chi }^{2/5},\mathrm{ } \frac{w_{mk}}{w_{uk}}={\chi }^{1/5},\mathrm{ } \frac{p_{mk}}{p_{uk}}={\chi }^{-1/5},\mathrm{ } \frac{q_{mk}}{q_{uk}}={\chi }^{-2/5},$$E23

$$\frac{{{\mathcal{G}}_m}_{mk}}{{{\mathcal{G}}_m}_{uk}}={\chi }^{-2/5},\mathrm{ } \frac{{{\mathcal{G}}_U}_{mk}}{{{\mathcal{G}}_U}_{uk}}={\chi }^{-6/5},\mathrm{ } \frac{{{\mathcal{G}}_a}_{mk}}{{{\mathcal{G}}_a}_{uk}}={\chi }^{-2/5}.$$E24

In addition to the above mentioned factors, the ratio of dimensionless initial defect length ${\gamma }_o$ to the dimensionless wellbore radius $\mathcal{A}$ also effects the hydraulic fracture initiation. It can be concluded here that the problem of axisymmetric hydraulic fracture depends only on three dimensionless parameters: the timescale ratio $\chi$, dimensionless wellbore radius $\mathcal{A}$ and the ratio of the initial defect length to the wellbore radius $\frac{{\mathcal{l}}_{\mathrm{o}}}{a}=\frac{{\mathrm{\gamma }}_{\mathrm{o}}}{\mathcal{A}}$.

3.1. Field and laboratory conditions: Scaling

Laboratory experiments are typically performed in order to study one particular aspect of a problem independently. Results of these experiments are then complemented by numerical and theoretical studies and ultimately verified by field experiments. Conditions in the laboratory experiments should be controlled in such a manner that they represent as close as possible a scaled version of the field conditions to be investigated. The scaling and the dimensionless parameters, described in the previous section, are thus critical in identifying the key parameters to simulate the right physics in the laboratory [18].

We will compare the different scales (in the Compressibility-Toughness scaling $UK$) of the problem for “typical” laboratory and field conditions in order to illustrate their differences and the need to carefully design laboratory experiments. It is worthwhile to recall that these scaling are based on the assumption of an axisymmetric fracture initiating from a circular notch. Parameters representative from the Barnett shale (5) are again considered for both the field conditions and the laboratory sample for illustration.

A typical wellbore diameter ($2a$) in the field is $8.75 in$ and a typical wellbore diameter in the laboratory is $1 in$. The constant flow rate $Q_o$ in the field is about $20 barrels/min$ and in the laboratory setting is around $5 ml/min$. The pressurization rate $\beta$ before breakdown in the field will be taken as $60 psi/s$ due to higher pumping rates in the field whereas in the laboratory it is about $1.2 psi/s$. The wellbore compressibility $U$ results from the compressibility of the fluid in the wellbore and the injection lines, as well as the compressibility of the wellbore and injection lines themselves. It is expressed as the ratio between the constant flow rate and the pressurization rate $U=Q_o/\beta$. The typical fluid used in the field is slick water [19] with a viscosity of $1 cP$. We consider glycerin as the fluid used in the laboratory with viscosity $1000 cP$.

The compressibility length scale $L_{uk}$ and time scale $t_{uk}$ corresponding to these parameters are displayed in Table 1. It is shown in the previous section that the dimensionless parameters depend upon the two dimensionless numbers i.e. $\mathcal{A}$ and $\chi$. These dimensionless numbers are also given in Table 1. While $L_{uk}$ and $t_{uk}$ give the length and the timescale associated with the compressibility effects, the value of $\chi$ determine the infinite medium propagation regime after the dissipation of compressibility effects. Large value of $\chi$ i.e. ($\chi >1$) means that the fracture will propagate in the viscosity regime whereas, smaller values of $\chi$ i.e. ($\chi <1$) means that the fracture will propagate in the toughness regime after the dissipation of compressibility effects.

		Field	Laboratory
Length scale of transition from the compressibility effects to the infinite medium propagation	$L_{uk}$(ft)	$56$	$2.4$
Timescale of transition from the compressibility effects to the infinite medium propagation	$t_{uk}$(sec)	$3$	$743$
Ratio of the timescales	$\chi$	$14$	$0.0036$
Dimensionless wellbore radius	$\mathcal{A}$	$0.0065$	$0.017$

Table 1.

Characteristic length scales and dimensionless parameters for the field and the laboratory conditions.

It is obvious from Table 1, that in the field, the fracture propagates in the viscosity dominated regime ($\chi >1$) whereas in the laboratory for the parameters chosen here, the fracture propagates in the toughness dominated regime ($\chi <1$) after the dissipation of the early-time compressibility effects. For the field conditions, the compressibility length scale is equal to $56 ft$. Which means that the propagation in the field is dominated by the compressibility effect until the fracture reaches about $56 ft$ (i.e. 150 times the wellbore radius) in to the formation. Even though this is a large value as compared to the wellbore radius, the length scale is still small as compared to the final fracture length which may be of the order of $800 \text{to} 1000 ft$ in the field. For laboratory conditions, the length scale is $2.4 ft$ (i.e. about 60 times the wellbore radius) which is smaller as compared to the field conditions, still the length-scale is large enough as compared to the specimen size that entire fracture propagation is dominated by the injection system compressibility.

The previous example has emphasized the differences with field conditions for a particular set of experimental parameters. However, these laboratory parameters can be appropriately adjusted in order to study a given regime of propagation. There can be different goals for an experimental campaign. For example, if the goal is to study hydraulic fracture propagation then the compressibility effects must be reduced in order to speed up the convergence to the infinite medium solution. These effects can be reduced by manipulating the material of the test block, using smaller injection lines and by using needle control valves as it was done for example in [20,21].

4. Numerical algorithm

The governing equations are solved in their dimensionless form described in Section 3. The elasticity equation is discretized by Displacement Discontinuity Method (DDM) using piecewise constant elements with the tip element correction [22] for better accuracy. The fluid flow is discretized by the finite volume method.

At time $t_n$, the fracture length is given as ${\gamma }_n$, fracture opening as ${\mathrm{\Omega }}_n$ and net pressure ${\mathrm{\Pi }}_n$. The algorithm consists of two nested loops. In the outer loop (time-stepping loop), a fracture increment $\mathrm{\Delta }\gamma =\mathrm{\Delta }\rho$ is specified and the time step $\mathrm{\Delta }{\mathrm{\tau }}_{\mathrm{n}}$, for which the tip asymptotic condition (15) is satisfied, is found iteratively. In the inner loop (i.e. the Reynolds solver), the coupled system of elasticity and lubrication equations are solved in order to find the fracture opening and pressure profiles corresponding to a given fracture length and (trial) time step. Details of the time stepping loop are given in Appendix B and the Reynolds Solver is described in Appendix C.

5. Results and discussion

The problem of the initiation of an axisymmetric hydraulic fracture from a wellbore depends only on three dimensionless parameters i.e. the timescale ratio $\chi$, dimensionless wellbore radius $\mathcal{A}$ and the ratio of the initial defect length and the wellbore radius $\frac{{\mathcal{l}}_{\mathrm{o}}}{a}=\frac{{\mathrm{\gamma }}_{\mathrm{o}}}{\mathcal{A}}$. The parameter space for these three dimensionless quantities is now explored. We aim to see their effect on the transition to the analytical solution of a penny-shaped hydraulic fracture propagating in an infinite medium [15,17] as well as on the breakdown pressure (i.e. the maximum pressure recorded).

In Figure. 2, the evolution of fracture length $\gamma$, inlet opening $\mathrm{\Omega }\left(0\right)$ and wellbore pressure ${\mathrm{\Pi }}_{wb}$ are plotted for different values of the timescale ratio $\chi$ for values of $\mathcal{A}=0.1$ and $\frac{{\gamma }_o}{\mathcal{A}}=0.3$. This is done in order to investigate the transition to the infinite medium solution. These results are presented in the Viscosity-Toughness ($MK$) scaling. The infinite medium solution goes from the viscosity dominated dissipation ($M$ vertex solution) to the toughness dominated dissipation ($K$ vertex solution). This transition is called the $MK$ edge of the semi-infinite medium solution. It can be observed in Figure. 2 that fracture length, opening and wellbore pressure of our numerical solution accounting for the wellbore effects converges at large time to the $MK$ edge solution in an infinite medium. The timescales ratio $\chi$ governs where the transient effects end along the MK edge of the infinite medium solution. For larger value of $\chi$, the transient effects end on the viscosity asymptote of the infinite medium solution, whereas for small value of $\chi$, the transient effect ends on the toughness asymptote of the infinite medium solution. The fact that the infinite medium solution is recovered at large time ultimately validates the large time behavior obtained by our numerical algorithm.

Figure 2.

Effect of $\chi$ on the convergence of fracture length $\mathrm{\gamma }$, inlet opening $\mathrm{\Omega }$ and wellbore pressure ${\mathrm{\Pi }}_{wb}$ to the infinite medium propagation solution for ${\gamma }_o=0.003, \mathcal{A}=0.1$. The results are displayed in the $MK$ scaling. The solid blue lines correspond to the simulation results, the dashed blue line is the $\mathcal{M}$ vertex solution, the dashed red line is the $\mathcal{K}$ vertex solution and the solid magenta line is the $\infty -$ medium solution.

It is also observed that the pressure and opening do not converge monotonically to the asymptotic solution. This is a characteristic feature of the near wellbore solution also obvious from Figure. 3 where the results are displayed in the Compressibility-Toughness ($UK$) scaling. Such a non-monotonic behavior is associated with fracture breakdown and is more pronounced for larger injection system compressibility, i.e. when the release of the stored fluid in the newly created fracture is more sudden.

The effects of $\chi$ on the breakdown pressure (defined as the maximum wellbore pressure), fracture length and effective flux entering the fracture are shown in Figure. 3 in the Compressibility-Toughness ($UK$) scaling. It has been observed that due to the strong fluid-solid coupling, the pressure in the wellbore keeps rising even after the fracture has already initiated. The pressure at which a fracture starts to propagate is called the initiation pressure and the highest pressure recorded is called the breakdown pressure. This difference in the initiation and the breakdown pressure has been observed theoretically [5,4] as well as experimentally [23] and it depends upon the fluid viscosity, injection rate and the system compressibility. It is observed from Equation (20) that low value of the timescale ratio $\chi$ corresponds to low viscosity and high injection system compressibility. It can be seen in Figure. 3 that the initiation pressure remains similar while higher breakdown pressures are obtained for higher values of $\chi$ and lower breakdown pressures are obtained for lower values of $\chi$. The breakdown is also much more abrupt for the case of low values of $\chi$. Such an abrupt breakdown corresponds to the unstable crack growth in the limiting case of a inviscid fluid ($\chi <{10}^{-4}$). Our numerical results are plotted along the inviscid solution of [6] in Figure. 4, where we can see that the numerical solution converges to the inviscid fluid solution for small $\chi$. Similar results for the case of a plane-strain hydraulic fracture are reported in [5].

Figure 3.

Effect of $\chi$ on the wellbore pressure ${\mathrm{\Pi }}_{wb}$, fracture length $\mathrm{\gamma }$ and effective flux entering the fracture for ${\gamma }_o=0.003, \mathcal{A}=0.1$ in $UK$ scaling. The blue coloured lines represent $\chi ={10}^{-4}$, magenta colour lines represent $\chi ={10}^3$ and the golden colour lines represent $\chi ={10}^4$

Figure 4.

Unstable crack growth in the case of inviscid fluid for $\chi =1.5\times {10}^{-4},\mathcal{ }\mathcal{A}=0.1$ in $UK$ scaling. The blue coloured lines represent ${\gamma }_o=0.005$, magenta colour lines represent ${\gamma }_o=0.02$ and the golden colour lines represent ${\gamma }_o=0.32$. The dashed magenta line represents the $K$-vertex solution whereas the dashed green line represents the inviscid fluid solution of [6].

Lastly, the value of $\frac{{\gamma }_{ϵ}}{\mathcal{A}}$ is plotted against the values of $\chi$ in Figure. 5, where ${\gamma }_{ϵ}$ is defined as the frature length $\gamma$ for which the infinite medium solution has been reached within a given tolerence ${ϵ}_{\infty }$ where

$${ϵ}_{\mathrm{\infty }}=\frac{\left|{\gamma }_{\mathrm{\infty }}-\gamma \right|}{\left|{\gamma }_{\mathrm{\infty }}\right|},$$E25

here ${\gamma }_{\infty }$ is the fracture radius for the penny shaped fracture propagation in an infinite medium. It can be seen from Figure. 5, that ${\gamma }_{ϵ}$ varies, for that case, from 11 to 13 times the wellbore radius for different values of $\chi$. This variation shows that there is minimal effect of $\chi$ on the length to wellbore ratio at which the hydraulic fracture is no longer affected by near-wellbore effects.

Let us now compare the effect of different ratios of $\frac{{\gamma }_o}{\mathcal{A}}=\frac{{\mathcal{l}}_o}{a}$ i.e. the ratio of the initial defect length to the wellbore radius. In Figure. 6, the dimensionless wellbore pressure, fracture length and effective flux entering the fracture are plotted for various ratios $\frac{{\gamma }_o}{\mathcal{A}}$ for $\chi =1$ and $\mathcal{A}=0.1$ in the Compressibility-Toughness ($UK$) scaling. It can be seen that higher breakdown pressures are obtained for lower ratios of $\frac{{\gamma }_o}{\mathcal{A}}$ and lower breakdown pressures are obtained for higher values of $\frac{{\gamma }_o}{\mathcal{A}}$. For lower values of $\frac{{\gamma }_o}{\mathcal{A}}$, there is a sudden drop in pressure after the breakdown similar to the behavior observed for a low viscosity fluid/ highly compressible injection system.

Figure 5.

Effect of $\chi$ on the convergence to infinite space solution for $\frac{{\gamma }_o}{\mathcal{A}}=0.003$ and $\mathcal{A}=0.1$.

It can be seen from Figure. 6, that there is no significant difference in convergence to the infinite space solution for different ratios of $\frac{{\gamma }_o}{\mathcal{A}}$.

Figure 6.

Effect of $\frac{{\gamma }_o}{\mathcal{A}}$ on the dimensionless wellbore pressure ${\mathrm{\Pi }}_{wb}$, fracture length $\mathrm{\gamma }$ and effective flux entering the fracture for $\chi =1, \mathcal{A}=0.1$ in $UK$ scaling. The blue coloured lines represent $\frac{{\gamma }_o}{\mathcal{A}}=1$, golden colour lines represent $\frac{{\gamma }_o}{\mathcal{A}}=0.05$ and the magenta colour lines represent $\frac{{\gamma }_o}{\mathcal{A}}=0.01$.

Figure 7.

Effect of $\mathcal{A}$ on the wellbore pressure ${\mathrm{\Pi }}_{wb}$, fracture length $\mathrm{\gamma }$ and effective flux entering the fracture for $\chi =1, \frac{{\gamma }_0}{\mathcal{A}}=0.1$ in $UK$ scaling. The blue coloured lines represent $\mathcal{A}=1.2$, magenta colour lines represent $\mathcal{A}=0.4$ and the golden colour lines represent $\mathcal{A}=0.05$

Finally, in Figure. 7, the effect of the dimensionless wellbore radius is considered for a fixed ratio $\frac{{\gamma }_0}{\mathcal{A}}=0.1$ and $\frac{{\gamma }_0}{\mathcal{A}}=0.01$. The results are displayed in the Compressibility-Toughness ($UK$) scaling. It is observed that the breakdown pressure is higher for a smaller dimensionless wellbore radius. It is to be noticed that the convergence of only one solution to the infinite medium solution is shown in the figure. All other solutions also converge to the infinite medium solution but at a larger time which is not shown here in order to focus on the breakdown phase. The effect of $\mathcal{A}$ on the convergence to the infinite space solution is displayed on Figure. 8. It can be seen that for dimensionless wellbore radius greater than 0.2, ${\gamma }_{ϵ}$ is only about 6 to 8 times $\mathcal{A}$. For wellbore radius smaller than 0.2, an exponential increase is observed on the transitional fracture length. The case of small $\mathcal{A}$ corresponds to conditions where a large amount of fluid can be stored in the wellbore during the pressurization phase, i.e $a\ll {\left(E^{\text{'}}U\right)}^{\frac{1}{3}}$. For those cases, the release of this stored volume of fluid in the fracture after its initiation results in an effective flow rate entering the crack much larger the injected one. This is the causes of such a large transition length to the infinite medium propagation for small values of $\mathcal{A}$.

Figure 8.

Effect of $\mathcal{A}$ on the convergence to infinite space solution for $\frac{{\gamma }_o}{\mathcal{A}}=0.1$. Two convergence criteria are used i.e ${ϵ}_{\infty }=0.1$ and ${ϵ}_{\infty }=0.01$.

6. Conclusions

The initiation of axisymmetric hydraulic fractures from a horizontal wellbore has been investigated with an emphasis on near wellbore effects. Through a detailed dimensional analysis, two characteristic timescales were identified. The first characteristic time $t_{mk}$ defines the timescale of transition from viscosity dominated propagation to the toughness dominated propagation and the second characteristic time $t_{uk}$ defines the timescale of transition from near wellbore effects to the infinite medium propagation. The ratio of these timescales $\chi$ (see Equation (20)) increases with increasing fluid viscosity and decreases with increasing injection system compressibility. The large time behavior of the numerical algorithm was verified by the convergence of the numerical solution to the propagation solution in an infinite medium for a sufficiently large fracture compared to the wellbore size. The effect of the timescale ratio $\chi$ on the convergence to the infinite medium solution was investigated. It was found that the numerical solution converges to the infinite medium solution for each value of $\chi$. The value of $\chi$ dictates on which infinite medium regime of propagation, the transient solution converges to. The solution converges to the toughness dominated propagation regime for small values of $\chi$, whereas it converges to the viscosity dominated regime for large values of $\chi$.

It was also found that the near wellbore effects are present up to a fracture length about 12 times the wellbore radius for different values of $\chi$. This shows that the transitional length is not affected by the value of $\chi$. The variation of the ratio of initial defect length to the wellbore radius $\frac{{\gamma }_o}{\mathcal{A}}=\frac{{\mathcal{l}}_o}{a}$ has a minimal effect on convergence to the infinite medium solution. In contrast, a small dimensionless wellbore radius $\mathcal{A}$ has a profound effect on convergence to the infinite medium solution. For small dimensionless wellbore radius $\mathcal{A}=\frac{a}{{\left(E^{\text{'}}U\right)}^{1/3}}$, a large increase in the transitional length is observed. This behavior is due to the injection system compressibility. The larger the compressibility, the larger the volume of fluid stored during pressurization and ultimately, the larger is the effective flux of fluid entering the fracture at breakdown compared to the nominal injection rate. This has the consequence of delaying the transition toward the solution of a hydraulic fracture propagating in an infinite medium under constant injection rate. It is important to note, however, that the presence of valves/perforations in the injection system will help dissipate the energy associated with such compressibility effects observed for small $\mathcal{A}$.

Appendix A

The edge dislocation kernel $J\left(r,r^{\text{'}}\right)$ mentioned in Equation (2) is given by [11] as follows

$$J(r,r^{\mathrm{\text{'}}})=r\left[R\left(r,r^{\mathrm{\text{'}}}\right)+S\left(r,r^{\mathrm{\text{'}}}\right)\right]$$E26

$$R\left(r,r^{\mathrm{\text{'}}}\right)=\left\{\begin{array}{c}\frac{1}{{\left(r^2-r^{\mathrm{\text{'}}2}\right)}^2}E\left(\frac{r}{r^{\mathrm{\text{'}}}}\right), r<r^{\mathrm{\text{'}}}\\ \frac{r}{r^{\mathrm{\text{'}}}}\frac{1}{r^2-r^{\mathrm{\text{'}}2}}E\left(\frac{r}{r^{\mathrm{\text{'}}}}\right)-\frac{1}{r{r}^{\mathrm{\text{'}}}}K\left(\frac{r}{r^{\mathrm{\text{'}}}}\right), r>r^{\mathrm{\text{'}}}\end{array}\right.$$E27

$$S\left(r,r^{\mathrm{\text{'}}}\right)=\underset{0}{\overset{\mathrm{\infty }}{\int }}\left[P\left(\xi ,x^{\mathrm{\text{'}}}\right)\alpha \left(\xi ,x\right)+Q\left(\xi ,x^{\mathrm{\text{'}}}\right)\beta \left(\xi ,x\right)\right]/\mathrm{\Delta }\left(\xi \right)d\xi ,$$E28

where $E$ and $K$ are complete elliptic integrals of first and second kind respectively and

$$P\left(\xi ,r^{\mathrm{\text{'}}}\right)={\xi }^2{I}_0{K}_1\left(\xi r^{\mathrm{\text{'}}}\right)-{\xi }^2{r}^{\mathrm{\text{'}}}I{K}_0\left(\xi r^{\mathrm{\text{'}}}\right),$$E29

$$\alpha \left(\xi ,r\right)=-\left[2\left(2-\nu \right)+{\xi }^2\right]K_1{K}_0\left(\xi r\right)/\xi +\xi x{K}_0{K}_1\left(\xi r\right)+r{K}_1{K}_1\left(\xi r\right)-3K_0{K}_0\left(\xi r\right),$$E30

$$Q\left(\xi ,r^{\mathrm{\text{'}}}\right)=-\left[2\left(1-\nu \right)+{\xi }^2\right]I_1{K}_1\left(\xi r^{\mathrm{\text{'}}}\right)+\xi I_0{K}_1\left(\xi r^{\mathrm{\text{'}}}\right)-\xi r^{\mathrm{\text{'}}}I{K}_0\left(\xi r^{\mathrm{\text{'}}}\right)+{\xi }^2{r}^{\mathrm{\text{'}}}I{K}_0\left(\xi r^{\mathrm{\text{'}}}\right),$$E31

$$\beta \left(\xi ,r\right)=2K_1{K}_0\left(\xi r\right)+\xi r{K}_1{K}_1\left(\xi r\right)+\xi K_0{K}_0\left(\xi r\right),$$E32

where $I_n$ and $K_n$ are the modified Bessel functions of first and second kind respectively where $n=\mathrm{0,1}$. In Equation (26), $S$ is a semi-infinite integral given by Equation (28) which has a very slow rate of convergence. In order to improve the convergence, we follow a method similar to [11]: $S$ is taken out of the integral and evaluated separately in closed form as follows

$$S_{\left\{new\right\}}\left(r,r^{\mathrm{\text{'}}}\right)=\underset{0}{\overset{\mathrm{\infty }}{\int }}\left|\left[P\left(\xi ,r^{\mathrm{\text{'}}}\right)\alpha \left(\xi ,r\right)+Q\left(\xi ,r^{\mathrm{\text{'}}}\right)\beta \left(\xi ,r\right)\right]/\mathrm{\Delta }\left(\xi \right)-A^{\mathrm{*}}\left(\xi ,r,r^{\mathrm{\text{'}}}\right)\right|d\xi +\underset{0}{\overset{\mathrm{\infty }}{\int }}{A}^{\mathrm{*}}\left(\xi ,r,r^{\mathrm{\text{'}}}\right)d\xi ,$$E33

where $A^{*}\left(\xi ,r,r^{\text{'}}\right)$ is the third order Taylor expansion of $\frac{\left(P\alpha +Q\beta \right)}{\mathrm{\Delta }}$ at infinity. The integral of $A^{*}\left(\xi ,r,r^{\text{'}}\right)$ can be obtained analytically.

Appendix B

The time step is computed by imposing the LEFM tip asymptote (15) in a weak form in the tip element. In the case of negligible toughness, care should be taken as the governing equations degenerate in the tip region. An algorithm based on the LEFM asymptote then requires very fine mesh for good convergence (see [24] for discussion).

The asymptotic volume of the tip element corresponding to the LEFM crack opening asymptote Eq. (15) is given as

$$\widehat{\mathrm{V}}={\mathcal{G}}_{\mathrm{k}}\underset{0}{\overset{\mathrm{\Delta }\rho }{\int }}{\widehat{\rho }}^{1/2}d\widehat{\rho }=\frac{2}{3}{\mathcal{G}}_{\mathrm{k}}\mathrm{\Delta }{\mathrm{\rho }}^{3/2}$$E34

where $\widehat{\rho }=\left(\gamma +\mathcal{A}\right)-\rho$ is the distance from the fracture tip and $\mathrm{\Delta }\mathrm{\rho }$ is the size of the tip element. The new time step $\mathrm{\Delta }{\tau }^{n+1}$ is found by finding the zero of the function $F$ defined as the mismatch between the current tip volume and the volume consistent with the LEFM asymptote:

$$F={\mathrm{V}}_{\mathrm{N}}-\widehat{\mathrm{V}}$$E35

where ${{\mathrm{V}}_{\mathrm{N}}=\mathrm{\Delta }\mathrm{\rho }\times \mathrm{\Omega }}_{\mathrm{N}}$ is the volume of the last element of the fracture. For each trial value of the time step, the coupled fluid-solid equations are solved in order to obtain a new estimate of the opening in the tip element ${\mathrm{\Omega }}_{\mathrm{N}}$. A secant method is used to find the zero of F. The iterative loop is repeated until the time step converges and the zero of (35) is found. The convergence criteria is

$$\frac{{‖\mathrm{\Delta }{\tau }^{n+1}-\mathrm{\Delta }{\tau }^n‖}_2}{{‖\mathrm{\Delta }{\tau }^n‖}_2}<ϵ,\mathrm{ } \text{with}\mathrm{ } ϵ={10}^{-4}.\mathrm{ }$$E36

Appendix C

The coupled problem of fluid-solid coupling inside the fracture is described as: for a given fracture length $\mathrm{\gamma }$, and time $\mathrm{\tau }+\mathrm{\Delta }{\tau }^n$, find the opening ${\mathrm{\Omega }}^{\mathrm{n}}+\mathrm{\Delta }{\mathrm{\Omega }}^n$. The continuity equation (12) is discretized using the Finite Volume Method as follows:

$$\mathrm{\Delta }{\mathrm{\Omega }}_{\mathrm{i}}=\frac{{\mathcal{G}}_r\mathrm{\Delta }\tau }{{\rho }_i\mathrm{\Delta }\rho }\left({\rho }_{i-1/2}{\mathrm{\Psi }}_{i-1/2}-{\rho }_{i+1/2}{\mathrm{\Psi }}_{i+1/2}\right)$$E37

where ${\rho }_i$ is the element midpoint coordinate and ${\rho }_{i-1/2}$ and ${\rho }_{i+1/2}$ are the element end point coordinates. The flux entering an element is denoted by ${\mathrm{\Psi }}_{i-1/2}$ while the outgoing flux is denoted by ${\mathrm{\Psi }}_{i+1/2}$. These fluxes are computed using the Poiseuille equation as follows:

$${\mathrm{\Psi }}_{\mathrm{i}-1/2}=-K_{i-1/2}\frac{{\mathrm{\Pi }}_{f\left(i\right)}-{\mathrm{\Pi }}_{f(i-1)}}{\mathrm{\Delta }\rho }$$E38

$${\mathrm{\Psi }}_{\mathrm{i}+1/2}=-K_{i+1/2}\frac{{\mathrm{\Pi }}_{f(i+1)}-{\mathrm{\Pi }}_{f\left(i\right)}}{\mathrm{\Delta }\rho }$$E39

Where the entrance hydraulic conductivity of an element $K_{i-1/2}$ is given by

$$K_{i-1/2}=\frac{1}{{\mathcal{G}}_m}{\left(\frac{{\mathrm{\Omega }}_i+{\mathrm{\Omega }}_{i-1}}{2}\right)}^3.$$E40

The flux entering the fracture from the wellbore is given by the inlet boundary condition

$${\mathrm{\Psi }}_{1/2}=\frac{1}{2\pi }\left({{\mathcal{G}}_q\mathrm{\Psi }}_o-{\mathcal{G}}_U\frac{\partial {\mathrm{\Pi }}_b}{\partial \tau }\right)$$E41

while no flow is assumed out of the fracture tip

$${\mathrm{\Psi }}_{\mathrm{N}}=0.$$E42

The elasticity equation (11) is discretized by the Displacement Discontinuity Method (DDM) to the following form

$${\mathit{A}}_{ij}\left({\mathrm{\Omega }}_j^0+\mathrm{\Delta }{\mathrm{\Omega }}_j\right)+{\mathrm{\Pi }}_h={\mathrm{\Pi }}_{f\left(i\right)}$$E43

where ${\mathit{A}}_{ij}$ is the stiffness matrix which is dense and symmetric for a regular mesh. Now putting the expression for ${\mathrm{\Pi }}_{f\left(i\right)}$ from Equation (43) in (38) and (39) and then $\mathrm{\Psi }$ from Equations (38) and (39) in Equation (37) we get the following expression after rearranging the terms

$$\mathrm{\Xi }\left(\mathrm{\Delta }\mathrm{\Omega }\right)\cdot \mathrm{\Delta }\mathrm{\Omega }=\mathrm{\Gamma }\left(\mathrm{\Delta }\mathrm{\Omega }\right)$$E44

where $\mathrm{\Xi }$ is expressed as

$$\mathrm{\Xi }={\delta }_{ij}-B_{ij}$$E45

where ${\delta }_{ij}$ is the Kronecker delta and

$$B_{ij}=\alpha \left\{\begin{array}{c}\frac{{\mathcal{G}}_U\mathrm{\Delta }\rho }{2\pi \mathrm{\Delta }\tau {\rho }_i}{A}_{\left(i,j\right)}+K_{e\left(i\right)}{A}_{\left(i+1,j\right)}-K_{e\left(i\right)}{A}_{\left(i,j\right)}, i=1 j=1,N\\ K_{w\left(i-1\right)}{A}_{\left(i-1,j\right)}-\left(K_{w\left(i-1\right)}+K_{e\left(i\right)}\right)A_{\left(i,j\right)}+K_{e\left(i\right)}{A}_{\left(i+1,j\right)}, i=2,N-1, j=1,N\\ K_{w\left(i-1\right)}{A}_{\left(i-1,j\right)}-K_{w\left(i-1\right)}{A}_{\left(i,j\right)}, i=N, j=1,N\end{array}\right.$$E46

$$\mathrm{\Gamma }=\alpha \left\{\begin{array}{c}\frac{{\mathcal{G}}_q{\mathrm{\Psi }}_o\mathrm{\Delta }\rho }{2\pi {\rho }_i}{A}_{\left(i,j\right)}+K_{e\left(i\right)}\left(A_{\left(i+1,j\right)}-A_{\left(i,j\right)}\right)\cdot {\mathrm{\Omega }}_j^0, i=1, j=1,N\\ B_{ij}\cdot {\mathrm{\Omega }}_j^0, i=2,N-1, j=1,N\\ {\mathrm{\Psi }}_m\Delta \rho +B_{Nj}\cdot {\mathrm{\Omega }}_j^0, i=N, j=1,N\end{array}\right.$$E47

where $\mathrm{\alpha }=\frac{{\mathcal{G}}_r\mathrm{\Delta }\tau }{\mathrm{\Delta }{\rho }^2}$ and

$$K_{\mathrm{e}}={\frac{{\rho }_{i-1/2}}{{\rho }_i}K}_{i-1/2},\mathrm{ } i=2,N\mathrm{ }$$E48

$$K_{\mathrm{e}}={\frac{{\rho }_{i+1/2}}{{\rho }_i}K}_{i+1/2},\mathrm{ } i=1,N-1$$E49

The nonlinear system (45) is solved by fixed point iteration

$$\mathrm{\Xi }\left(\mathrm{\Delta }{\mathrm{\Omega }}^{\mathrm{k}}\right)\cdot \mathrm{\Delta }{\mathrm{\Omega }}^{\mathrm{k}+1}=\mathrm{\Gamma }\left(\mathrm{\Delta }{\mathrm{\Omega }}^{\mathrm{k}}\right)$$E50

with under-relaxation

$$\mathrm{\Delta }{\mathrm{\Omega }}^{\mathrm{k}+1}=\left(1-{\eta }^{\gamma }\right)\mathrm{\Delta }{\mathrm{\Omega }}^{\mathrm{k}}+{\eta }^{\gamma }\mathrm{\Delta }{\mathrm{\Omega }}^{\mathrm{k}+1}$$E51

where $\eta$ is the relaxation parameter. The convergence criteria is

$$\frac{{‖\mathrm{\Delta }{\mathrm{\Omega }}^{k+1}-\mathrm{\Delta }{\mathrm{\Omega }}^k‖}_2}{{‖\mathrm{\Delta }{\mathrm{\Omega }}^{k+1}‖}_2}<ϵ,\mathrm{ } \text{with}\mathrm{ } ϵ={10}^{-5}.\mathrm{ }$$E52