2.1. Problem description

The waterflooding problem is analysed within the framework of the 2-D model sketched in figure 1. It consists of an infinite poroelastic domain with a circular hole of radius $a$, constrained to deform under plane strain conditions and subjected to a far-field isotropic compressive stress $\sigma _{0}$ and a far-field pore pressure $p_{0}<\sigma _{0}$. The porous material is saturated by a Newtonian fluid. It is also assumed to be perfectly brittle but with a negligible toughness, i.e. $K_{Ic}=0$. There are six independent parameters to describe the Newtonian fluid and the poroelastic material, namely: dynamic viscosity $\mu$, Young's modulus $E$, Poisson's ratio $\nu$, Biot coefficient $\alpha$, permeability $k$ or mobility $\kappa =k/\mu$, and diffusivity $c$. Three derived parameters are introduced to simplify the problem formulation,

(2.1a–c)\begin{equation} E'=\frac{E}{1-\nu^{2}},\quad \eta=\frac{\alpha(1-2\nu)}{2(1-\nu)},\quad \mu'=12\mu, \end{equation}

where $E'$ denotes the plane strain Young's modulus, and $\eta \in [0,1/2]$ is a poroelastic stress coefficient (Detournay & Cheng Reference Detournay and Cheng1993).

Figure 1. (a) A schematic diagram of the bi-wing fracture near a borehole and (b) poroelastic medium $\mathcal {D}$, wellbore $\mathcal {W}$ and bi-wing crack $\mathcal {C}$ with corresponding field variables.

With the hole initially filled by the same fluid at pressure $p_{0}$, the system is initially equilibrated with a uniform pore pressure $p_{0}$. There is an initial elastic stress concentration at the hole boundary on account of the difference between the fluid pressure in the hole and the far-field stress. At time $t=0$, fluid is injected in the hole at a constant rate $Q_{0}$ (dimension $L^{2}T^{-1}$) causing a progressive change in the stress and pore pressure field in the vicinity of the hole that eventually lead to the initiation at the circular boundary of a symmetric bi-wing hydraulic fracture. The direction of fracture is associated with a small anisotropy of the far-field stress that causes the crack to propagate in the direction perpendicular to the least compressive far-field stress. Propagation of the fracture is tracked by the distance $\ell$ between the crack tip and the hole centre. This distance, a monotonic function of time $t$, will simply be referred to as the crack length.

The primary objective of the analysis is to determine the evolution of the borehole pressure $p_{w}$ and the crack length $\ell$, as well as the dependence of the functions $p_{w}(t)$ and $\ell (t)$ on the various parameters describing this system.

2.2. Assumptions

The mathematical model is constructed on the following two critical hypotheses: (i) the crack propagates in a region surrounding the borehole, where the pore pressure field is quasi-steady; (ii) fluid storage in the crack is negligible compared with the amount of fluid lost by leak-off. The justification for these two hypotheses lays in the presumed large permeability of the rock, which is assumed to be poorly consolidated. Such a rock is also assumed to have a negligible fracture toughness.

As a consequence of the first hypothesis, the diffusion equation governing the evolution of the pore pressure field degenerates into a Laplace equation in a region containing the crack. This degeneracy transforms the two-way poroelastic coupling between the stress and pore pressure fields to a one-way coupling; i.e. the pore pressure field can be solved first, and then acts as a forcing term in the elasticity equation governing the stress field. Although the solution still depends on time, due to the far-field asymptotic solution of the diffusion equation, the combined hypotheses lead to a history-independent solution. In other words, time $t$ becomes an independent parameter of the solution rather than a variable, as further discussed in § 3.3.

2.3. Field variables in Poroelastic domain $\mathcal {D}$ and crack $\mathcal {C}$

A 2-D cartesian coordinates system $(x,y)$ centred on the borehole is defined with the $x$-axis oriented along the crack. Two sub-domains are naturally introduced: the 2-D poroelastic domain $\mathcal {D}=\{(x,y)\,|\,x^{2}+y^{2}\geqslant a^{2}\}$ and the crack $\mathcal {C}=\{(x,y)\mid x\in [-\ell ,-a]\cup [a,\ell ], y=0\}$. Domain $\mathcal {D}$ is bounded by wellbore $\mathcal {W}=\{(x,y)\,|\,x^{2}+y^{2}=a^{2}\}$ and crack $\mathcal {C}$. The field variables defined on $\mathcal {D}$ are stress $\boldsymbol {\boldsymbol {\sigma }}$, displacement $\boldsymbol {\boldsymbol {u}}$, pore pressure $p$ and specific discharge $\boldsymbol {\boldsymbol {v}}$. On $\mathcal {C}$, the field variables are fluid pressure $p_{f}$, leak-off $g$, crack aperture $w$ and flux $q$. All these variables are functions of time $t$. The variables are constrained by the following continuity and jump conditions between the two sub-domains:

(2.2)\begin{equation} \left.\begin{gathered} p_{f}(x,t)=p(x,0,t),\quad p_{f}(x,t)={-}\sigma_{yy}(x,0,t),\quad \sigma_{xy}(x,0,t)=0,\\ w(x,t)=[{u_{y}(x,0,t)]},\quad g(x,t)=[{v_{y}(x,0,t)]}. \end{gathered}\right\} \end{equation}

Here $a< | x| < \ell$, and $[ {f(x,0,t)]}\equiv f(x,0^{+},t)-f(x,0^{-},t)$.

The field variables in $\mathcal {D}$ are governed by the theory of poroelasticity, while those in $\mathcal {C}$ are also governed by the lubrication equation and by the fracture propagation criterion. The complete set of governing equations and boundary conditions are presented in §§ 2.4–2.7.

2.4. Governing equations on $\mathcal {D}$

The equations of poroelasticity can be reduced to a set of two coupled partial differential equations that govern the displacement field $\boldsymbol {\boldsymbol {u}}(\boldsymbol {\boldsymbol {x}},t)$ in the porous solid and the pore pressure field $p(\boldsymbol {\boldsymbol {x}},t)$: namely, a Navier-type equation for $\boldsymbol {\boldsymbol {u}}(\boldsymbol {\boldsymbol {x}},t)$ with a body force term proportional to the pore pressure gradient, and a diffusion equation for $p(\boldsymbol {\boldsymbol {x}},t)$ with a source term proportional to the rate of change of $\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {\boldsymbol {u}}$ (Cheng Reference Cheng2016). The coupling term in the diffusion equation vanishes, however, when the solution reaches a steady state or when the displacement field is irrotational and the medium is infinite (Detournay & Cheng Reference Detournay and Cheng1993). As elaborated in more details in § 3.1, the first condition is indeed met in the problem in view of the a priori hypothesis that the crack is growing in a region where the flow is in a pseudo steady-state.

Thus, in the near-field the pore pressure is governed by Laplace equation

(2.3)\begin{equation} \kappa\nabla^{2}p={-}g(x,t)\delta(y), \end{equation}

where $\delta (y)$ denotes the Dirac delta function, and the leak-off $g(x,t)$ is part of the solution. While in the far field, $p(x,t)$ is given by the solution of

(2.4)\begin{equation} \frac{\partial p}{\partial t}-c\nabla^{2}p=Q_{0}\delta(\boldsymbol{\boldsymbol{x}})H(t), \end{equation}

where $H(t)$ is the Heaviside function. This asymptotic solution actually corresponds to the classical continuous source solution (Cheng Reference Cheng2016). The specific discharge $\boldsymbol {\boldsymbol {v}}$ is related to the pore pressure $p$ according to Darcy's law

(2.5)\begin{equation} \boldsymbol{\boldsymbol{v}}={-}\kappa\boldsymbol{\nabla} p. \end{equation}

The Navier equation for the displacement $\boldsymbol {\boldsymbol {u}}$ is given by

(2.6)\begin{equation} G\nabla^{2}\boldsymbol{\boldsymbol{u}}+\frac{G}{1-2\nu}\boldsymbol{\nabla}(\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\boldsymbol{u}})-\alpha\boldsymbol{\nabla} p=0, \end{equation}

where $G=E/2(1+\nu )$ is the shear modulus. The stress $\boldsymbol {\boldsymbol {\sigma }}$ is related to pore pressure $p$ and strain $\boldsymbol {\boldsymbol {\varepsilon }}=(\boldsymbol {\nabla }\boldsymbol {\boldsymbol {u}}+\boldsymbol {\boldsymbol {u}}\boldsymbol {\nabla })/2$ according to

(2.7)\begin{equation} \boldsymbol{\boldsymbol{\sigma}}+\alpha p\boldsymbol{\boldsymbol{I}}=2G\boldsymbol{\boldsymbol{\varepsilon}}+ \frac{2G\nu}{1-2\nu}(\boldsymbol{\boldsymbol{\varepsilon}}:\boldsymbol{\boldsymbol{I}})\boldsymbol{\boldsymbol{I}}, \end{equation}

where $\boldsymbol {\boldsymbol {I}}$ denotes the second-order identity tensor.

The system (2.3)–(2.7) represent the complete set of equations governing the fields $\boldsymbol {\boldsymbol {u}}(\boldsymbol {\boldsymbol {x}},t)$, $p(\boldsymbol {\boldsymbol {x}},t)$, $\boldsymbol {\boldsymbol {\sigma }}(\boldsymbol {\boldsymbol {x}},t)$, $\boldsymbol {\boldsymbol {v}}(\boldsymbol {\boldsymbol {x}},t)$ in $\mathcal {D}$.

2.5. Boundary conditions on $\mathcal {W}$

The injection pressure $p_{w}$ and the fraction $(1-\varPhi$) of the injection rate $Q_{0}$ directly entering the rock through the borehole wall (to be discussed near (3.21)) are both a priori unknown functions of time $t$. The borehole pressure $p_{w}$ represents a boundary condition on wellbore $\mathcal {W}$ for both the stress and the pore pressure,

(2.8a,b)\begin{equation} \boldsymbol{\boldsymbol{\sigma}}(\boldsymbol{\boldsymbol{x}},t)\boldsymbol{\cdot}\boldsymbol{\boldsymbol{n}}={-}p_{w}\boldsymbol{\boldsymbol{n}},\quad p(\boldsymbol{\boldsymbol{x}},t)=p_{w},\quad \boldsymbol{\boldsymbol{x}}\in\mathcal{W}, \end{equation}

where $\boldsymbol {\boldsymbol {n}}$ denotes the unit normal external vector on $\mathcal {W}$. The rate of fluid volume passing through the borehole wall and the flux at the crack inlet are constrained by the total injection $Q_{0}$,

(2.9)\begin{equation} -\oint_{\mathcal{W}}\boldsymbol{\boldsymbol{v}}(\boldsymbol{\boldsymbol{x}},t)\boldsymbol{\cdot}\boldsymbol{\boldsymbol{n}}\,{\rm d} s+2q(a,t)=Q_{0}, \end{equation}

where the symmetry condition has been taken into account in (2.9).

2.6. Governing equations on $\mathcal {C}$

The equations governing the fluid flow in the crack $\mathcal {C}$ are Poiseuille's law

(2.10)\begin{equation} q={-}\frac{w^{3}}{\mu'}\frac{\partial p_{f}}{\partial x}, \end{equation}

and the continuity equation

(2.11)\begin{equation} \frac{\partial q}{\partial x}+g=0, \end{equation}

with the storage term neglected in accord with the simplifying assumption stated earlier. Combining (2.11) and (2.10) yields the Reynolds lubrication equation

(2.12)\begin{equation} \frac{1}{\mu'}\frac{\partial}{\partial x}\left({w^{3}\frac{\partial p_{f}}{\partial x}}\right)=g. \end{equation}

Two additional boundary conditions are required for the lubrication equation. One condition imposes the pressure at the fracture inlet, and the other one a vanishing flux at the crack tip (Detournay & Peirce Reference Detournay and Peirce2014),

(2.13)$$\begin{gather} p_{f}({\pm} a,t)=p_{w}, \end{gather}$$

(2.14)$$\begin{gather}q({\pm}\ell(t),t)=0. \end{gather}$$

2.7. Tip asymptotics

According to linear elastic fracture mechanics (LEFM), the crack aperture asymptotically behaves near the advancing tip as

(2.15)\begin{equation} w(x)\propto(\ell- |x| )^{{3}/{2}},\quad x\to\pm\ell, \end{equation}

if the fracture toughness $K_{Ic}=0$ and provided that the fluid pressure is finite at the tip. This latter condition is indeed fulfilled as the fluid pressure in the crack is continuous with the pore pressure field, which must be regular as it is governed by the Laplace equation in a finite region (Evans Reference Evans2010). (Since a 2-D point source leads to a logarithmically singular pore pressure, the pore pressure induced by a distributed leak-off at the crack tip is thus regular as can be confirmed by convolving the source function with the leak-off.) Hence, the fluid pressure $p_{f}$ at the tip can be expanded as

(2.16)\begin{equation} p_{f}(x)=c_{0}+c_{1}({\ell- |x| })+O ({\ell- |x| })^{2},\quad x\to\pm\ell, \end{equation}

which satisfies the requirement of the crack aperture tip asymptotic solution (2.15).

Substituting the above tip asymptotics for $p_{f}$ and $w$ into the lubrication equation (2.12), shows that the leak-off $g$ near the tip should behave as

(2.17)\begin{equation} g\propto(\ell- |x| )^{{7}/{2}},\quad x\to\pm\ell. \end{equation}

The tip asymptotic solution outlined above differs from tip asymptotic solutions for hydraulic fractures propagating in permeable media, constructed on the assumptions that leak-off is either governed by the one-dimensional Carter law (Lenoach Reference Lenoach1995; Garagash, Detournay & Adachi Reference Garagash, Detournay and Adachi2011; Detournay Reference Detournay2016) or by the diffusion equation (Detournay & Garagash Reference Detournay and Garagash2003). Solutions built on assuming Carter leak-off (Howard & Fast Reference Howard and Fast1957) predict $p_{f}(x)$ to be singular (with the singularity depending on whether the fracture propagates in the viscosity- or the toughness-dominated regime) and leak-off rate to have a square root singularity. On the other hand, asymptotic solutions obtained on the basis of the diffusion equation require the existence of a lag region, which is filled by pore fluid circulating in and out of the cavity.

3.1. Pore pressure field

The pore pressure field is formulated as a superposition of three particular solutions, each satisfying the condition that the pore pressure on the borehole boundary $\mathcal {W}$ is uniform,

(3.1)\begin{equation} p(x,y,t)=p_{0}+p_{i}(x,y,t)+p_{l}(x,y,t), \end{equation}

where $p_{i}(x,y,t)$ is the pore pressure induced by continuous injection from the borehole in the absence of a fracture, and $p_{l}(x,y,t)$ the pore pressure field associated with leak-off from the fracture. The field $p_l$ does not contribute to the total flow rate entering the porous medium $\mathcal {D}$, as explained later.

The field $p_{i}$ is given by the classical source solution of the diffusion equation (Carslaw & Jaeger Reference Carslaw and Jaeger1959)

(3.2)\begin{equation} p_{i}(x,y,t)=\frac{Q_{0}}{4{\rm \pi}\kappa}E_{1}\left({\frac{r^{2}}{4ct}}\right). \end{equation}

Inside the quasi-stationary region which expands as $\chi \sqrt {ct}$, the diffusion equation effectively degenerates to Laplace equation (2.3). In this region the exponential integral function $E_1$ simplifies to (Abramowitz & Stegun Reference Abramowitz and Stegun1972)

(3.3)\begin{equation} p_{i}(x,y,t)\approx{-}\frac{Q_{0}}{2{\rm \pi}\kappa}\left({\ln\frac{r}{2\sqrt{ct}}+\frac{\gamma}{2}}\right),\quad r<\chi\sqrt{ct}, \end{equation}

with $r=\sqrt {x^{2}+y^{2}}$ and $\gamma =0.577216\cdots$ denoting the Euler gamma constant. The number $\chi \approx 0.35$ defines the conditions for which the asymptotic solution (3.3) applies within an error less than 1 %.

The field $p_{l}(x,y,t)$ breaks the axial symmetry of the injection-induced pore pressure $p_{i}(x,y,t)$ by accounting for leak-off from the fracture. This field, which is also assumed to satisfy the Laplace equation, is constructed by distributing fluid sources along the crack $\mathcal {C}$ and image sinks inside the borehole $\mathcal {W}$ so that there is no net fluid injection in the poroelastic media $\mathcal {D}$. The locations of the image sinks inside $\mathcal {W}$ are chosen so that the pore pressure $p_{l}$ is uniform on $\mathcal {W}$. Thus, the field $p_{l}(x,y,t)$ is expressed as a convolution integral of the leak-off $g$ with the singular kernel $\tilde {P}(x,y,s,a)$ (Gao & Detournay Reference Gao and Detournay2020a), i.e.

(3.4)\begin{equation} p_{l}(x,y,t)=\frac{1}{2{\rm \pi}\kappa}\int_{a}^{\ell}g(s,t)\tilde{P}(x,y,s,a)\,{\rm d}s. \end{equation}

The kernel $\tilde {P}(x,y,s,a)$ accounts for the problem symmetry with respect to the $y$-axis by taking the form

(3.5)\begin{equation} \tilde{P}(x,y,s,a)=P(x,y,s,a)+P({-}x,y,s,a), \end{equation}

with $P$ representing the pore pressure field generated by a source located at $(s,0)$ and an image sink positioned at $(a^{2}/s,0)$, $s>a$. On $y=0$, the kernel $P$ is simply given by

(3.6)\begin{equation} P(x,0,s,a)={-}\ln |x-s|+\ln \left| {x-\frac{a^{2}}{s}}\right| . \end{equation}

Expression (3.4) for the pore pressure field $p_{l}$ ensures that there is no contribution from this field to the total flow rate entering the domain in $\mathcal {D}$. In other words,

(3.7)\begin{equation} \oint_{\mathcal{W}}\boldsymbol{\boldsymbol{v}}_{l}\boldsymbol{\cdot}\boldsymbol{\boldsymbol{n}}\,{\rm d}s=2q(a,t). \end{equation}

With the introduction of the image sink, the same flux entering the fracture inlet is coming back through the borehole boundary $\mathcal {W}$, so that $(1-\varPhi )Q_0$ is leaking through $\mathcal {W}$ and $\varPhi Q_0$ through the walls of crack $\mathcal {C}$. After superposing all the fields, the rate of fluid volume injected into the media is $Q_0$. Here $\varPhi$ denotes the fraction of injected fluid leaking through $\mathcal {C}$; an expression to calculate $\varPhi$ is given in (3.21).

In summary, the pore pressure field induced by injection $Q_{0}$ and leak-off $g$ is found by combining (3.1), (3.2) and (3.4). On the crack ($a\leqslant |x| \leqslant \ell , y=0$) the fluid pressure reads as

(3.8)\begin{equation} p_{f}(x,t)=p_{0}-\frac{Q_{0}}{2{\rm \pi}\kappa}\left( {\ln\frac{ |x| }{2\sqrt{ct}} + \frac{\gamma}{2}}\right)+\frac{1}{2{\rm \pi}\kappa}\int_{a}^{\ell}g(s,t)\tilde{P}(x,0,s,a)\,{\rm d}s. \end{equation}

To simplify the notation, the second parameter of $\tilde {P}$ is omitted in the rest of the paper, if it is zero, i.e. $\tilde {P}(x,s,a)\equiv \tilde {P}(x,0,s,a)$.

3.2. Stress field

Taking advantage of the linearity of the elasticity equations (2.6) and (2.7), the stress field $\boldsymbol {\boldsymbol {\sigma }}$ in $\mathcal {D}$ can also be constructed by superposition of particular solutions. Here, only the stress $\sigma _{yy}$ on the crack $\mathcal {C}$ is of concern; it is decomposed into four parts: (i) the in-situ stress $-\sigma _{0}$, (ii) the stress $\sigma _{yy}^{f}$ induced by the crack aperture $w$, (iii) the poroelastic stress $\sigma _{yy}^{p}$ induced by the pore pressure change $p-p_{0}$, and (iv) a stress $\sigma _{yy}^{c}$ resulting from enforcing the stress boundary condition (2.8a). Thus,

(3.9)\begin{equation} \sigma_{yy}(x,0,t)={-}\sigma_{0}+\sigma_{yy}^{f}+\sigma_{yy}^{p}+\sigma_{yy}^{c}. \end{equation}

The crack induced stress $\sigma _{yy}^{f}$ can be expressed as (Hills et al. Reference Hills, Kelly, Dai and Korsunsky1996)

(3.10)\begin{equation} \sigma_{yy}^{f}(x,t)=\frac{E'}{4{\rm \pi}}\int_{a}^{\ell}w(s,t)\frac{\partial}{\partial s}\tilde{H}(x,s,a)\,{\rm d}s, \end{equation}

where

(3.11)\begin{equation} \tilde{H}(x,s,a)\equiv H(x,s,a)+H({-}x,s,a), \end{equation}

noting also that $\tilde {H}(x,a,a)=0$ and $w(\ell )=0$. The kernel function $H(x,s,a)$,

(3.12)\begin{align} H(x,s,a) &= \frac{s^{2}-a^{2}}{sx^{2}}-\frac{a^{2}\left(s^{2}-a^{2}\right)}{s\left(sx-a^{2}\right)^{2}} \left(\frac{s^{2}}{a^{2}}-\frac{s^{2}-a^{2}}{sx-a^{2}}\right) \nonumber\\ &\quad -\frac{s}{sx-a^{2}}+\frac{a^{2}}{x^{3}}+\frac{1}{x-s}+\frac{1}{x}, \end{align}

represents the stress $\sigma _{yy}(x,0)$ induced by a unit normal dislocation located at $(s,0)$ along the $y$-axis that satisfies the boundary condition $\boldsymbol {\boldsymbol {\sigma }}\boldsymbol {\cdot }\boldsymbol {\boldsymbol {n}}=0$ on borehole $\mathcal {W}$ (Dundurs & Mura Reference Dundurs and Mura1964; Hills et al. Reference Hills, Kelly, Dai and Korsunsky1996).

Echoing the decomposition of $p$ in § 3.1, the poroelastic stress $\boldsymbol {\boldsymbol {\sigma }}^{p}$ is expressed as the sum of an injection-induced stress $\boldsymbol {\boldsymbol {\sigma }}_{i}^{p}$ and the leak-off induced stress $\boldsymbol {\boldsymbol {\sigma }}_{l}^{p}$. On crack $\mathcal {C}$, $\boldsymbol {\boldsymbol {\sigma }}_{i}^{p}$ is given by (Cheng & Detournay Reference Cheng and Detournay1998; Gao & Detournay Reference Gao and Detournay2020a)

(3.13)\begin{equation} \left.\begin{gathered} \sigma_{xx,i}^{p}(x,0,t) ={-}2\eta\frac{Q_{0}}{4{\rm \pi}\kappa} \left[{-\frac{\gamma}{2}-\ln\frac{ |x|}{2\sqrt{ct}}+ \frac{1}{2}}\right]={-}\eta p_{i}(x,0,t)-\eta\frac{Q_{0}}{4{\rm \pi}\kappa},\\ \sigma_{yy,i}^{p}(x,0,t) ={-}2\eta\frac{Q_{0}}{4{\rm \pi}\kappa} \left[{-\frac{\gamma}{2}-\ln\frac{ |x| }{2\sqrt{ct}}- \frac{1}{2}}\right]={-}\eta p_{i}(x,0,t)+\eta\frac{Q_{0}}{4{\rm \pi}\kappa},\\ \sigma_{xy,i}^{p}(x,0,t) =0. \end{gathered}\right\} \end{equation}

It can be shown, using the poroelastic steady-state source solution, that the leak-off induced stress $\boldsymbol {\boldsymbol {\sigma }}_{l}^{p}$ on the crack $\mathcal {C}$ is given by (Gao & Detournay Reference Gao and Detournay2020a)

(3.14a,b)\begin{equation} \sigma_{xx,l}^{p}(x,0,t)=\sigma_{yy,l}^{p}(x,0,t)={-}\eta p_{l}(x,0,t),\quad \sigma_{xy,l}^{p}(x,0,t)=0. \end{equation}

Combining (3.13), (3.14a,b) and (3.1) leads to the following expressions for the induced normal and shear stresses $\sigma _{yy}^{p}$ and $\sigma _{xy}^{p}$ on $\mathcal {C}$:

(3.15a,b)\begin{equation} \sigma_{yy}^{p}(x,0)={-}\eta(p_{f}(x)-p_{0})+\eta\frac{Q_{0}}{4{\rm \pi}\kappa},\quad \sigma_{xy}^{p}(x,0)=0,\quad a\leqslant|x| \leqslant\ell. \end{equation}

The normal and shear stresses $\{\sigma _{rr}^{p},\sigma _{r\theta }^{p}\}$ acting on borehole boundary $\mathcal {W}$ need also to be evaluated to determine $\sigma _{yy}^{c}$ in (3.9), as shown below. While it is clear that the borehole injection-induced component $\boldsymbol {\boldsymbol {\sigma }}_{i}^{p}$ of $\boldsymbol {\boldsymbol {\sigma }}^{p}$ satisfies $\boldsymbol {\boldsymbol {\sigma }}_{i}^{p}(\boldsymbol {\boldsymbol {x}})\boldsymbol {\cdot }\boldsymbol {\boldsymbol {n}}=\sigma _{xx,i}^{p}(a,0)\boldsymbol {\boldsymbol {n}}$ on $\boldsymbol {\boldsymbol {x}}\in \mathcal {W}$, we assume that the leak-off induced normal stress on $\mathcal {W}$ is also uniform and given by $\boldsymbol {\boldsymbol {\sigma }}_{l}^{p}(\boldsymbol {\boldsymbol {x}})\boldsymbol {\cdot }\boldsymbol {\boldsymbol {n}}\approx \sigma _{xx,l}^{p}(a,0)\boldsymbol {\boldsymbol {n}}=-\eta p_{l}(a,0)\boldsymbol {\boldsymbol {n}}$ on $\boldsymbol {\boldsymbol {x}}\in \mathcal {W}$. With this assumption and (3.13), the expressions for $\sigma _{rr}^{p}$ and $\sigma _{r\theta }^{p}$ acting on $\mathcal {W}$ simplify to

(3.16a,b)\begin{equation} \sigma_{rr}^{p}(\boldsymbol{\boldsymbol{x}})={-}\eta(p_{w}-p_{0})-\eta\frac{Q_{0}}{4{\rm \pi}\kappa},\quad \sigma_{r\theta}^{p}(\boldsymbol{\boldsymbol{x}})=0,\quad \boldsymbol{\boldsymbol{x}}\in\mathcal{W}. \end{equation}

The approximation adopted for $\boldsymbol {\boldsymbol {\sigma }}_{l}^{p}$ recognizes that $\boldsymbol {\boldsymbol {\sigma }}_{l}^{p}$ is negligible compared with $\boldsymbol {\boldsymbol {\sigma }}_{i}^{p}$ for a short crack, as only a small portion of fluid is leaking from the crack; while for a long crack, the actual stress boundary conditions on borehole $\mathcal {W}$ are effectively irrelevant at the scale of the fracture.

Finally, the fourth term $\sigma _{yy}^{c}$ in (3.9) is obtained by enforcing the boundary condition (2.8a). After considering the in-situ stresses $\sigma _0$ and the poroelasticity induced normal stress $\sigma _{rr}^{p}$ on borehole $\mathcal {W}$ in (3.16a), $\sigma _{yy}^{c}$ is given by (Cheng Reference Cheng2016)

(3.17)\begin{equation} \sigma_{yy}^{c}(x)=\left[{p_{w}-\sigma_{0}-\eta(p_{w}-p_{0})-\eta\frac{Q_{0}}{4{\rm \pi}\kappa}}\right]\frac{a^{2}}{x^{2}}. \end{equation}

In summary, the normal stress $\sigma _{yy}$ on crack $\mathcal {C}$ is obtained by substituting (3.10), (3.15a) and (3.17) into (3.9) to give

(3.18)\begin{align} \sigma_{yy}(x,0) &={-}\sigma_{0}+\left[{p_{w}-\sigma_{0}-\eta(p_{w}-p_{0})-\frac{\eta Q_{0}}{4{\rm \pi}\kappa}}\right] \frac{a^{2}}{x^{2}} \nonumber\\ &\quad -\eta(p_{f}(x)-p_{0})+\frac{\eta Q_{0}}{4{\rm \pi}\kappa}+\frac{E'}{4{\rm \pi}} \int_{a}^{\ell}w(s)\frac{\partial}{\partial s}\tilde{H}(x,s,a)\,{\rm d}s. \end{align}

3.2.1. Influence of a far-field deviatoric stress

The assumption of an isotropic far-field stress $\sigma _0$ could be relaxed by introducing the minimum and maximum in-situ stresses $\sigma _h$ and $\sigma _H$ in the $y$- and $x$-directions, respectively. This is equivalent to adding an independent parameter, the deviatoric stress, $S_0=(\sigma _H-\sigma _h)/2$, into the model.

However, the only difference introduced by $S_0$ is adding the term

(3.19)\begin{equation} -S_{0}\left({1-3\frac{a^{2}}{x^{2}}}\right)\frac{a^{2}}{x^{2}} \end{equation}

into (3.17) and (3.18), as well as changing $\sigma _0$ to $\sigma _h$ in (3.18).

This additional term (3.19) affects the fracture initiation pressure $p_{wi}$, but it is not changing the analysis otherwise. In fact, by substituting $\ell =a$ and the Terzaghi effective stress $\sigma _{yy}(a,0)+p_b=0$ into (3.18), the fracture initiation pressure now reads as

(3.20)\begin{equation} p_{wi}=\frac{\sigma_h-S_0-\eta p_0}{1-\eta}. \end{equation}

This result is the well-known Haimson–Fairhurst (H–F) breakdown criterion (Haimson & Fairhurst Reference Haimson and Fairhurst1967) with zero tensile strength. Since $S_0$ essentially affects only the fracture initiation pressure, we have assumed that $S_0=0$.

3.3. Reduced system of equations

As shown above, the problem can be entirely formulated in terms of the crack length $\ell (t)$, the crack aperture $w(x,t)$, the fracture pressure $p_{f}(x,t)$ and the leak-off $g(x,t)$. This reduction to variables defined only along crack $\tilde {\mathcal {C}}$ is achieved by (i) expressing the pore pressure field $p$ and the stress field $\boldsymbol {\boldsymbol {\sigma }}$ in domain $\mathcal {D}$ as convolution integrals of distributed singularities over crack $\tilde {\mathcal {C}}$, (ii) accounting for the problem symmetry, and (iii) enforcing the continuity conditions (2.2) $p=p_{f}=-\sigma _{yy}$ on $\tilde {\mathcal {C}}$.

The system consisting of the two integral equations (3.8) and (3.18), together with the lubrication equation (2.12) and boundary conditions (2.13)–(2.15) is closed. Since time $t$ does not appear in a differential operator, $t$ is actually a parameter and not a variable. Although the problem is formally history-independent, the variation of the solution with time should be consistent; in particular, the crack length $\ell$ is expected to be a monotonically increasing function of time. To emphasize the demotion of $t$ from a variable to a parameter, the dependence on $x$ and $t$ of the field variables defined on $\tilde {\mathcal {C}}$ is now denoted as $(x;t)$.

Once $\ell (t)$, $w(x;t)$, $p_{f}(x;t)$ and $g(x;t)$ have been solved at time $t$, the field variables in $\mathcal {D}$ can be determined using convolution integrals. Also, the flooding efficiency $\varPhi (t)\equiv 2q(a;t)/Q_{0}$, defined as the portion of injected fluid leaking to the rock through the crack, can be expressed as

(3.21)\begin{equation} \varPhi(t)=\frac{2}{Q_{0}}\int_{a}^{\ell}g(x;t)\,{\rm d}x, \end{equation}

after making use of the lubrication equation (2.12) and the crack tip boundary condition (2.14).

6.1. Radial-flow regime ($R$-regime)

We start the asymptotic analysis by presenting solutions for pressure $\varPi$ and crack length $\varLambda$ in the rock-flow regime, or $R$-regime, which corresponds to the radial-flow edge in the conceptual phase diagram sketched in § 5. First, we show that after rescaling, the solution only depends on two parameters: namely, time ${{\bar {\tau }}}$ and poroelastic coefficient $\eta$. Indeed, rescaling the variables according to

(6.1a–f)\begin{equation} {{\bar{\tau}}}=\alpha_k^{{-}2}\tau,\quad {{\bar{\varPi}}}=\varPi,\quad {{\bar{\varOmega}}}=\alpha_k^{{-}1}\varOmega,\quad {{\bar{\varGamma}}}=\alpha_k^{{-}1}\varGamma,\quad {{\bar{\varLambda}}}=\varLambda,\quad {{\bar{\varPhi}}}=\alpha_k^{{-}2}\varPhi, \end{equation}

leads to this alternative formulation of the governing equations,

(6.2)$$\begin{gather} \alpha_\tau^{2}\frac{\partial}{\partial \xi}\left({{{\bar{\varOmega}}}^{3}\frac{\partial{{\bar{\varPi}}}}{\partial \xi}}\right)={{\bar{\varGamma}}}, \end{gather}$$

(6.3)$$\begin{gather}(1-\eta){{\bar{\varPi}}}={-}\left[{(1-\eta){{\bar{\varPi}}}_{w}-\frac{\eta}{2{\rm \pi}}}\right]\frac{\alpha_\tau^{2}}{\xi^{2}}-\frac{\alpha_\tau}{4{\rm \pi}} \int_{\alpha_\tau}^{1}{{\bar{\varOmega}}}(\zeta,{{\bar{\tau}}})\frac{\partial}{\partial \zeta}\tilde{H}({\xi,\zeta,\alpha_\tau})\,{\rm d}\zeta, \end{gather}$$

(6.4)$$\begin{gather}{{\bar{\varPi}}}={-}\frac{1}{2{\rm \pi}}\ln\left({2{{\bar{\varLambda}}}{\xi}}\right)+\frac{\alpha_k^{2}}{2{\rm \pi}\alpha_\tau} \int_{\alpha_\tau}^{1}{{\bar{\varGamma}}}(\zeta,{{\bar{\tau}}})\tilde{P}({\xi,\zeta,\alpha_\tau})\,{\rm d}\zeta, \end{gather}$$

and of the flooding efficiency

(6.5)\begin{equation} {{\bar{\varPhi}}}={-}2\alpha_\tau{{\bar{\varOmega}}}^{3}(\alpha_\tau){{\bar{\varPi}}}'(\alpha_\tau). \end{equation}

The boundary conditions are similar to (4.19a,b) and omitted here. Note that the scaled crack length ${{\bar {\varLambda }}}$ is hidden in the borehole radius $\alpha _\tau$ in the governing equations, and is a part of the solution.

Since the second term on the right-hand side of (6.4) is negligible compared with the first term on account that $\alpha _k\ll 1$, (6.4) reduces to

(6.6)\begin{equation} {{\bar{\varPi}}}(\xi;{{\bar{\tau}}},\eta)={-}\frac{1}{2{\rm \pi}}\ln\left[{2{{\bar{\varLambda}}}({{\bar{\tau}}},\eta){\xi}}\right], \end{equation}

and on the borehole boundary $\xi =\alpha _\tau$, the pressure is simply given by

(6.7)\begin{equation} {{\bar{\varPi}}}_{w}({{\bar{\tau}}})={-}\frac{1}{2{\rm \pi}}\ln{\frac{2}{\sqrt{{{\bar{\tau}}}}}}. \end{equation}

Noting that $\alpha _\tau =({{\bar {\varLambda }}}\sqrt {{{\bar {\tau }}}})^{-1}$, the above equations prove that the solution in the $R$-regime only depends on ${{\bar {\tau }}}$ and $\eta$.

Fracture initiation at the borehole wall corresponds to $\alpha _\tau =1$ according to its definition in (4.2), since fracture toughness is zero and $\ell =a$ and ${{\bar {\varOmega }}}=0$. Substituting these conditions into the elasticity equation (6.3) yields the fracture initiation pressure

(6.8)\begin{equation} {{\bar{\varPi}}}_{wi}=\frac{\eta}{4{\rm \pi}(1-\eta)}. \end{equation}

This result is consistent with the H–F breakdown criterion (Haimson & Fairhurst Reference Haimson and Fairhurst1967) for the particular case of negligible tensile strength and isotropic far-field stress. As explained in § 3.2.1, in the case of anisotropic far-field stress there would be an additional term proportional to the far-field stress deviator $S_0$ in expression (6.8) for the fracture initiation pressure. Time ${{\bar {\tau }}}_{i}$ and fracture length ${{\bar {\varLambda }}}_{i}$ at initiation are then deduced from (6.8), (6.7) and definition (6.1a–f) of ${{\bar {\tau }}}$,

(6.9a,b)\begin{equation} {{\bar{\tau}}}_{i}=4\exp\left({{\frac{\eta}{1-\eta}}}\right),\quad {{\bar{\varLambda}}}_{i}=\frac{1}{2}\exp\left({{-\frac{\eta}{2(1-\eta)}}}\right), \end{equation}

with ${{\bar {\tau }}}_{i}\in [4,4e]$ and ${{\bar {\varLambda }}}_{i}\in [1/(2\sqrt {e}),1/2]$ since $\eta \in [0,1/2]$.

Finally, the fracture aperture ${{\bar {\varOmega }}}(\xi ;{{\bar {\tau }}},\eta )$ and crack length ${{\bar {\varLambda }}}({{\bar {\tau }}},\eta )$ in the $R$-regime are determined by the elasticity singular integral equation (6.3) with ${{\bar {\varPi }}}(\xi )$ given by (6.6), together with propagation criterion (4.19a,b). These equations are solved using the numerical technique described in Appendix A. The computed fracture length ${{\bar {\varLambda }}}({{\bar {\tau }}},\eta )$ is plotted in figure 3 for $\eta =\{0,0.25,0.5\}$.

6.2. Small-time asymptote ($E$-vertex)

In the early time following fracture initiation, when $\ell -a\ll a$, the hydraulic fracture can be treated as an edge crack in a semi-infinite plane with uniform net pressure. Note that the crack opening $w$ at the mouth of an edge crack is given by $w=\beta {\rm \Delta} p\ell '/E'$ (Stallybrass Reference Stallybrass1970), with $\beta =2.91$, ${\rm \Delta} p$ and $\ell '$ denoting the net pressure and the length of the edge crack, respectively. Hence, the small-time crack opening at the borehole wall reads as

(6.10)\begin{equation} {{\bar{\varOmega}}}_{wE}({{\bar{\tau}}},\eta)=\beta\left({{{\bar{\varLambda}}}\sqrt{{{\bar{\tau}}}}-1}\right){\rm \Delta}{{\bar{\varPi}}}_{wE}, \end{equation}

with ${\rm \Delta} {{\bar {\varPi }}}_{wE}=2(1-\eta ){{\bar {\varPi }}}_{w}-\eta /(2{\rm \pi} )$, where ${{\bar {\varPi }}}_{w}$ is determined with (6.7). The flooding efficiency ${{\bar {\varPhi }}}_{E}$ at the $E$-vertex is then deduced from (6.5) to be

(6.11)\begin{equation} {{\bar{\varPhi}}}_{E}({{\bar{\tau}}},\eta)=\frac{\beta^{3}}{\rm \pi}\left({{{\bar{\varLambda}}}\sqrt{{{\bar{\tau}}}}-1}\right)^{3}{\rm \Delta}{{\bar{\varPi}}}_{wE}^{3}. \end{equation}

6.3. Intermediate-time asymptote ($I$-vertex)

Provided that $t_{d}\ll t\ll t_{k}$ or, equivalently, ${{\bar {\tau }}}\gg 1$ and $\tau \ll 1$, the solution trajectory passes through the neighbourhood of the $I$-vertex. The existence of an intermediate-time asymptote hinges, therefore, on the condition that $t_{d}/t_{k}=\alpha _{k}^{2}\lll 1,$ which is effectively met for $t_{d}/t_{k}<5\times 10^{-3}$ according to numerical simulations. At the $I$-vertex, the fracture length is large compared with the borehole radius, but the fracture is still hydraulically invisible. Since $\ell \gg a$, the stress intensity factor can be calculated using the weight function for a Griffith crack (Bueckner Reference Bueckner1970),

(6.12)\begin{equation} K_{I}=2\sqrt{\rm \pi}\int_{0}^{1}\frac{{{\bar{\varPi}}}(\xi)}{\sqrt{1-\xi^{2}}}{\rm d}\xi. \end{equation}

With ${{\bar {\varPi }}}(\xi )$ given by (6.6) on account that the fracture does not affect the porous media flow, (6.12) yields

(6.13)\begin{equation} K_{I}={-}\frac{\sqrt{\rm \pi}}{2}\ln{{\bar{\varLambda}}}. \end{equation}

Thus, the assumption of zero fracture toughness implies that ${{\bar {\varLambda }}}=1$ at the $I$-vertex, as illustrated in figure 3 when ${{\bar {\tau }}}\to \infty$. The fluid pressure in the crack is then given by

(6.14)\begin{equation} {{\bar{\varPi}}}_{I}={-}\frac{1}{2{\rm \pi}}\ln 2{\xi}. \end{equation}

Furthermore, since $\alpha _\tau \ll 1$, the elasticity equation (6.3) simplifies to

(6.15)\begin{equation} (1-\eta){{\bar{\varPi}}}_{I}(\xi)={-}\frac{1}{4{\rm \pi}{{\bar{\tau}}}^{{1}/{2}}} \int_{{-}1}^{1}\frac{{{\bar{\varOmega}}}_{I}(\zeta)}{(\xi-\zeta)^{2}}\mathrm{d}\zeta, \end{equation}

which can be inverted to yield (Sneddon & Lowengrub Reference Sneddon and Lowengrub1969)

(6.16)\begin{equation} {{\bar{\varOmega}}}_{I}(\xi)=\frac{4}{\rm \pi}(1-\eta){{\bar{\tau}}}^{{1}/{2}} \int_{0}^{1}{{\bar{\varPi}}}_{I}(\zeta)\ln\left|\frac{\sqrt{1-\xi^{2}}+ \sqrt{1-\zeta^{2}}}{\sqrt{1-\xi^{2}}-\sqrt{1-\zeta^{2}}}\right|\mathrm{d}\zeta. \end{equation}

The explicit form of the fracture aperture ${{\bar {\varOmega }}}_{I}(\xi ;{{\bar {\tau }}},\eta )$ at the $I$-vertex can then be deduced from (6.16) with ${{\bar {\varPi }}}_{I}$ in (6.14),

(6.17)\begin{equation} {{\bar{\varOmega}}}_{I}=\sqrt{\frac{\rm \pi}{8}}(1-\eta){{\bar{\tau}}}^{{1}/{2}}\omega(\xi), \end{equation}

with

(6.18)\begin{equation} \omega(\xi)=\sqrt{1-\xi^{2}}-\frac{\rm \pi}{2}| \xi |+\xi\arctan\left({\frac{\xi}{\sqrt{1-\xi^{2}}}}\right), \end{equation}

and leak-off ${{\bar {\varGamma }}}_{I}(\xi ;{{\bar {\tau }}},\eta )$ is then obtained by integrating the lubrication equation (6.2) using (6.14) and (6.17),

(6.19)\begin{equation} {{\bar{\varGamma}}}_{I}(\xi;{{\bar{\tau}}},\eta)=\sqrt{\frac{\rm \pi}{2^{11}}}(1-\eta)^{3}{{\bar{\tau}}}^{{1}/{2}}\upsilon(\xi), \end{equation}

with

(6.20)\begin{equation} \upsilon(\xi)=\frac{\omega^{2}(\xi)}{\xi^{2}}\left({\rm \pi}\xi+\sqrt{1-\xi^{2}}-2\xi\arctan \frac{\xi}{\sqrt{1-\xi^{2}}}\right). \end{equation}

Finally, flooding efficiency ${{\bar {\varPhi }}}_{I}({{\bar {\tau }}},\eta )$ is deduced from (6.5),

(6.21)\begin{equation} {{\bar{\varPhi}}}_{I}=\sqrt{\frac{\rm \pi}{2^{9}}}(1-\eta)^{3}{{\bar{\tau}}}^{{3}/{2}}. \end{equation}

It can readily be confirmed that ${{\bar {\varOmega }}}_{I}\stackrel {\xi \to 1}{\sim }(1-\xi )^{{3}/{2}}$ and ${{\bar {\varGamma }}}_{I}\stackrel {\xi \to 1}{\sim }(1-\xi )^{{7}/{2}}$, consistent with the tip asymptotics discussed in § 2.7.

The crack aperture at the borehole, ${{\bar {\varOmega }}}_{wI}({{\bar {\tau }}},\eta )={{\bar {\varOmega }}}_{I}(\alpha _\tau ;{{\bar {\tau }}},\eta )$, is obtained by setting $\xi =\alpha _\tau ={{\bar {\tau }}}^{-({1}/{2})}$ in (6.17),

(6.22)\begin{equation} {{\bar{\varOmega}}}_{wI}=\sqrt{\frac{\rm \pi}{8}}(1-\eta){{\bar{\tau}}}^{{1}/{2}}, \end{equation}

while the leak-off at the crack inlet ${{\bar {\varGamma }}}_{wI}({{\bar {\tau }}},\eta )={{\bar {\varGamma }}}_{I}(\alpha _\tau ;{{\bar {\tau }}},\eta )$ is

(6.23)\begin{equation} {{\bar{\varGamma}}}_{wI}=\sqrt{\frac{\rm \pi}{2^{11}}}(1-\eta)^{3}{{\bar{\tau}}}^{{3}/{2}}. \end{equation}

The different time exponent in (6.23) for ${{\bar {\varGamma }}}_{wI}({{\bar {\tau }}},\eta )$ and in (6.19) for ${{\bar {\varGamma }}}_{I}(\xi ;{{\bar {\tau }}},\eta )$ is a consequence of the singular behaviour of the leak-off function near the origin, $\varGamma _{I}\sim \xi ^{-2}$, and the movement of borehole boundary $\alpha _\tau ={{\bar {\tau }}}^{-({1}/{2})}$ in the $\xi$-space. Singularity $\varGamma _{I}\sim \xi ^{-2}$ itself results from ${{\bar {\varGamma }}}\sim \partial ^{2}\varPi /\partial \xi ^{2}$ on account that ${{\bar {\varOmega }}}$ is uniform near $\xi =0$, and from the assumption that ${{\bar {\varPi }}}\sim \ln 2\xi$ holds. But this expression for ${{\bar {\varPi }}}$ is valid only provided that the porous media flow equation (6.4) is dominated by the borehole injection at the $I$-vertex. In other words, the leak-off induced term in (6.4),

(6.24)\begin{equation} G(\xi;\alpha_k,{{\bar{\tau}}},\eta)\equiv\sqrt{{{\bar{\tau}}}}\alpha_k^{2} \int_{\alpha_\tau}^{1}{{\bar{\varGamma}}}(\zeta,{{\bar{\tau}}},\eta)\tilde{P}({\xi,\zeta,\alpha_\tau})\,{\rm d}\zeta, \end{equation}

should be negligible compared with the injection term $\ln 2{\xi }$. However, $G(\xi )$ grows unbounded when $\alpha _\tau \to 0$ if ${{\bar {\varGamma }}}\sim \xi ^{-2}$. In fact, expanding $G(\xi )$ near the origin yields

(6.25)\begin{equation} G(\xi)\approx\sqrt{\frac{\rm \pi}{2}}\frac{\alpha_k^{2}}{24\xi^{2}}(1-\eta)^{3}\sqrt{{{\bar{\tau}}}},\quad \alpha_\tau\leqslant\xi\ll 1, \end{equation}

implying that $G(\alpha _\tau )\sim \alpha _k^{2}{{\bar {\tau }}}^{{3}/{2}}$. This result indicates the formation of a boundary layer at $\xi =\alpha _\tau$ when ${{\bar {\tau }}}\gg 1$. Under these conditions, the above expressions for ${{\bar {\varOmega }}}_{I}$, ${{\bar {\varPi }}}$ and ${{\bar {\varGamma }}}_{I}$ have to be understood as outer solutions. In fact, numerical simulations show that for ${{\bar {\tau }}}\gtrapprox 1.47/\alpha _k^{{4}/{3}}$ (corresponding to $\tau \gtrapprox 1.47\alpha _k^{{2}/{3}}$ and $\varPhi _{I}\gtrapprox 0.0698(1-\eta )^{3}$) causes $G(\xi )$ to be large enough near the boundary to break the dominance of the injection-induced component $\ln {2\xi }$ in the pressure profile. Growth of leak-off ${{\bar {\varGamma }}}$ with time causes a progressive increase of flooding efficiency ${{\bar {\varPhi }}}$, and eventually the departure of the solution from the $I$-vertex.

We have not made any attempt to construct an explicit solution in the boundary layer to be matched with the outer solution. Rather the solution was computed numerically using an element size small enough to capture the boundary layer. Numerically computed leak-off $\varGamma (\xi )$ and pressure $\varPi (\xi )$ can be found in figure 4 for $\alpha _k=10^{-7}$, $\tau =10^{-2}$ and $\eta =0$. The leak-off profile $\varGamma (\xi )$ illustrated in figure 4(a) confirms the existence of an inner boundary layer and of an intermediate asymptotic region where $\varGamma =\sqrt {{\rm \pi} \tau /2^{11}}\xi ^{-2}$. The tip asymptote $\varGamma =\sqrt {{\rm \pi} \tau }(1-\xi )^{{7}/{2}}/12$ can also be seen in this figure. Moreover, numerical simulations indicate that the boundary layer corresponds approximately to $\xi \in [\alpha _\tau ,0.562\tau ]$ with $\tau \in [1.47\alpha _k^{{2}/{3}},0.178]$, noting that the intermediate asymptote disappears for $\alpha _k>0.042$. The profile of $\varPi (\xi )$ in figure 4(b) shows a combination of the boundary layer and the outer solution ${{\bar {\varPi }}}(\xi )\sim \ln {2\xi }$.

Figure 4. Profiles of $\alpha _k=10^{-7}$, $\tau =10^{-2}$ and $\eta =0$ for (a) leak-off $\varGamma (\xi )$ and (b) fluid pressure $\varPi (\xi )$. The near borehole and near tip asymptotic solutions of $\varGamma (\xi )$ in the $I$-vertex, i.e. $\sqrt {{\rm \pi} /2048}\xi ^{-2}$ and $\sqrt {{\rm \pi} }(1-\xi )^{{7}/{2}}/12$, are compared in (a). The asymptotic solution of $\varPi (\xi )$, $-1/(2{\rm \pi} ) \ln {2\xi }$, is compared in (b). A boundary layer in $\xi <0.005$ is observed.

6.4. Large-time asymptote (fracture-flow regime, $F$-vertex)

At large time, the fracture becomes conductive enough that all the injected fluid leaks into the rock through the fracture walls, i.e. $\varPhi _{F}\simeq 1$, with the subscript $F$ used to denote a quantity defined at the $F$-vertex. Furthermore, numerical results confirm that $\varLambda _{F}=1$ when $\alpha _\tau \to 0$ ($\tau \to \infty$), implying that the crack grows as $\ell \propto \sqrt {t}$ like the radius of the steady-state region. In this section we formulate a time series expansion of the solution near the $F$-vertex. However, we show that this solution does not comply with the tip conditions for leak-off $\varGamma$ and aperture $\varOmega$ discussed in § 2.7, which leads to the existence of a boundary layer at the crack tip. The obtained asymptotic large-time solution, referred to as the $F0$-solution, is thus only self-similar in the bulk.

6.4.1. $F$-vertex solution and near $F$-vertex solution

The large-time limits $\varLambda \to 1$, $\varPhi \to 1$ and $\alpha _\tau \to 0$ result in a simplification of the system of (4.16)–(4.18). In particular, imposing that $2\sqrt {\tau }\int _{0}^{1}\varGamma _{F}\,{\rm d} \xi =1$, or $\varPhi _{F}=1$, enables the removal of the term proportional to $\ln \xi$ from the porous media flow equation (4.18). The system of equations (4.16)–(4.18) is therefore simplified to

(6.26)\begin{equation} \left.\begin{gathered} \frac{1}{\tau}\frac{\partial}{\partial \xi}\left({\varOmega_{F}^{3}\frac{\partial \varPi_{F}}{\partial \xi}}\right)=\varGamma_{F},\\ (1-\eta)\varPi_{F}={-}\frac{1}{2{\rm \pi}\sqrt{\tau}}\int_{0}^{1}\varOmega_{F}(\zeta,\tau) \frac{\xi^{2}+\zeta^{2}}{(\xi^{2}-\zeta^{2})^{2}}{\rm d}\zeta,\\ \varPi_{F}={-}\frac{\ln 2}{2{\rm \pi}}-\frac{\sqrt{\tau}}{2{\rm \pi}}\int_{0}^{1} \varGamma_{F}(\zeta,\tau)\ln | {\xi^{2}-\zeta^{2}} | \,{\rm d}\zeta. \end{gathered}\right\} \end{equation}

Replacing $\varGamma _{F}$, $\varPi _{F}$ and $\varOmega _{F}$ in the above system of equations by a regular asymptotic expansion of the form $\mathcal {S}_{F}=\mathcal {S}_{F0}\tau ^{-\beta _{0}}+\mathcal {S}_{F1}\tau ^{-\beta _{1}}$, and balancing terms of the same orders leads to the expansion

(6.27)\begin{equation} \left.\begin{gathered} \varGamma_{F}(\xi,\tau) =\varGamma_{F0}(\xi)\tau^{-({1}/{2})}+\varGamma_{F1}(\xi) \tau^{-({3}/{4})}+O(\tau^{{-}1}),\\ \varOmega_{F}(\xi,\tau) =\varOmega_{F0}(\xi)\tau^{{1}/{4}}+ \varOmega_{F1}(\xi)\tau^{0}+O(\tau^{-({1}/{4})}),\\ \varPi_{F}(\xi,\tau) =\varPi_{F0}(\xi)\tau^{-({1}/{4})}+\varPi_{F1}(\xi) \tau^{-({1}/{2})}+O(\tau^{-({3}/{4})}), \end{gathered}\right\} \end{equation}

with the zero-order solution governed by

(6.28)\begin{equation} \left.\begin{gathered} \frac{\partial}{\partial \xi}\left({\varOmega_{F0}^{3}\frac{\partial \varPi_{F0}}{\partial \xi}}\right)=\varGamma_{F0},\\ (1-\eta)\varPi_{F0}={-}\frac{1}{2{\rm \pi}}\int_{0}^{1}\varOmega_{F0}(\zeta) \frac{\xi^{2}+\zeta^{2}}{(\xi^{2}-\zeta^{2})^{2}}{\rm d}\zeta,\\ \int_{0}^{1}\varGamma_{F0}(\zeta){\ln |{\xi^{2}-\zeta^{2}}| } \,{\rm d}\zeta={-}\ln{2}, \end{gathered}\right\} \end{equation}

and for the first-order correction by

(6.29)\begin{equation} \left.\begin{gathered} \frac{\partial}{\partial \xi}\left({\varOmega_{F0}^{3}\frac{\partial \varPi_{F1}}{\partial \xi}+3\varOmega_{F0}^{2} \varOmega_{F1}\frac{\partial \varPi_{F0}}{\partial \xi}}\right)=\varGamma_{F1},\\ (1-\eta)\varPi_{F1}={-}\frac{1}{2{\rm \pi}}\int_{0}^{1}\varOmega_{F1}(\zeta) \frac{\xi^{2}+\zeta^{2}}{(\xi^{2}-\zeta^{2})^{2}}{\rm d}\zeta,\\ \varPi_{F0}={-}\frac{1}{2{\rm \pi}}\int_{0}^{1}\varGamma_{F1}(\zeta){\ln |{\xi^{2}-\zeta^{2}} | }\,{\rm d}\zeta. \end{gathered}\right\} \end{equation}

The last equation of the system (6.28) can be directly solved for $\varGamma _{F0}(\xi )$ to yield

(6.30)\begin{equation} \varGamma_{F0}=\frac{1}{{\rm \pi}\sqrt{1-\xi^{2}}}. \end{equation}

The first two equations in (6.28) to be solved for $\varOmega _{F0}(\xi )$ and $\varPi _{F0}(\xi )$ are actually similar to those ruling the propagation of a plane strain hydraulic fracture in the viscosity/leak-off dominated regime, with Carter leak-off. This solution – referred to as the $\tilde {M}$-solution (Adachi & Detournay Reference Adachi and Detournay2008; Hu & Garagash Reference Hu and Garagash2010) – is characterized by negligible storage of fluid in the fracture and negligible toughness. Its construction assumes, according to Carter's model (Howard & Fast Reference Howard and Fast1957), that the leak-off rate is inversely proportional to the time elapsed since the time of the first exposure, i.e. $g(x,t)\sim (t-t_{0}(x))^{-({1}/{2})}$, where $t_{0}$ is the time when the crack tip was at position $x$. The similarity between the $F0$-solution and the $\tilde {M}$-solution stems from the formal equivalence between expression (6.30) for leak-off $\varGamma _{F0}$ combined with a growth law $\ell \propto \sqrt {t}$ and Carter's leak-off model, even though the underlying physics is quite different.

Taking into account the difference in the scaling factors between the equations governing the $F0$- and $\tilde {M}$-solutions, the zero-order $\varOmega _{F0}(\xi ;\eta )$ and $\varPi _{F0}(\xi ;\eta )$ can be written as

(6.31a,b)\begin{equation} \varOmega_{F0}=(1-\eta)^{{1}/{4}}{\rm \pi}^{-({1}/{4})}\varOmega_{\tilde{M}}(\xi),\quad \varPi_{F0}=(1-\eta)^{-({3}/{4})}{\rm \pi}^{-({1}/{4})}\varPi_{\tilde{M}}(\xi), \end{equation}

where aperture $\varOmega _{\tilde {M}}$ and pressure $\varPi _{\tilde {M}}$ can be expressed as Frobenius expansions in terms of Gegenbauer polynomials (Adachi & Detournay Reference Adachi and Detournay2008). In particular, the leading terms of crack aperture and fluid pressure at the borehole are given by

(6.32a,b)\begin{equation} \varOmega_{wF0}=1.447(1-\eta)^{{1}/{4}}\tau^{{1}/{4}},\quad \varPi_{wF0}=0.424(1-\eta)^{-({3}/{4})}\tau^{-({1}/{4})}. \end{equation}

The first-order solution $\{\varGamma _{F1},\varPi _{F1},\varOmega _{F1}\}$ is then calculated by numerically solving the linear system of (6.29) using expressions (6.31a,b) for $\varOmega _{F0}(\xi )$ and $\varPi _{F0}(\xi )$. Since the lubrication equation in the first order is linearized by the time series expansion, it is easily solved. It can be verified that

(6.33a–c)\begin{equation} \varGamma_{F1}\propto(1-\eta)^{-({3}/{4})},\quad \varOmega_{F1}\propto(1-\eta)^{-({1}/{2})},\quad \varPi_{F1}\propto(1-\eta)^{-({3}/{2})}. \end{equation}

6.4.2. Crack tip boundary layer at the $F$-vertex

The zero-order leak-off (6.30) is singular at the crack tip $\varGamma _{F0}\sim (1-\xi )^{-({1}/{2})}$ and, thus, in contradiction with the expected non-singular crack tip asymptotic behaviour $\varGamma \sim (1-\xi )^{{7}/{2}}$; see (2.17). Furthermore, as discussed in § 2.7, the fracture propagation criterion (4.19a,b) requires crack aperture $\varOmega \sim (1-\xi )^{{3}/{2}}$ and fluid pressure $\varPi$ to be finite at the tip. However, these two requirements are not fulfilled by the zero-order solution (6.31a,b), which is characterized by tip asymptotes $\varOmega _{F0}\sim (1-\xi )^{{5}/{8}}$ and $\varPi _{F0}\sim (1-\xi )^{-({3}/{8})}$ (Adachi & Detournay Reference Adachi and Detournay2008). These discrepancies point to the existence of a boundary layer at the crack tip near the $F$-vertex, and also suggest that the tip asymptotes of the zero-order solution actually represent the intermediate asymptotes of the $F$-vertex solution.

Figure 5 compares the leak-off profile $\varGamma (\xi )$ numerically computed for $\tau =10^{4}$ and $\eta =0$ with the zero-order solution $\varGamma _{F0}\tau ^{-({1}/{2})}$ and the first two orders time expansion $\varGamma _{F0}\tau ^{-({1}/{2})}+\varGamma _{F1}\tau ^{-({3}/{4})}$. This figure confirms that the $F$-vertex asymptotic solution with the two leading terms can be understood as an (approximated) outer solution to the problem and that a boundary layer indeed exists at the crack tip. Numerical simulations indicate that the size of the boundary layer shrinks with increasing time as $\tau ^{-({1}/{4})}$.

Figure 5. Leak-off profile $\varGamma (\xi )$ obtained by $F$-vertex asymptotic solutions and numerical solutions for $\tau =10^{4}$ and $\eta =0$.

The numerical solutions calculated at $\tau =\{10^{6},10^{7}\}$ with $\eta =0$ are shown as functions of the distance from the crack tip ($1-\xi$) in the log-log plots of figure 6. To remove the asymptotic dependence on time of the solution, $\varGamma \tau ^{{1}/{2}}$ and $\varOmega \tau ^{-({1}/{4})}$ are plotted instead of $\varGamma$ and $\varOmega$. Figures 6(a) and 6(b) show the transition between the tip asymptotic solutions of the boundary layer and of the zero-order $F$-vertex series expansion, for $\varGamma$ and $\varOmega$, respectively. The results confirm that the autonomous tip asymptotes of the $F0$-solution indeed act as intermediate asymptotes, which shield the global solution from the tip boundary layer at large time. The numerical pressure profiles shown in figure 6(c) also indicate that pressure $\varPi$ is finite at the crack tip. To conclude this analysis of the large-time solution, it is worth pointing out that the existence of a tip boundary layer at the $F$-vertex can be traced to the conflict between the singular leak-off at the crack tip $\varGamma _{F0}\sim (1-\xi )^{-({1}/{2})}$ and the non-singular expected asymptote $\varGamma \sim (1-\xi )^{{7}/{2}}$.

Figure 6. Numerical results for (a) leak-off $\varGamma \tau ^{{1}/{2}}$, (b) aperture $\varOmega \tau ^{-({1}/{4})}$ and (c) fluid pressure $| {\varPi } |$ near the tip at the large time ($F$-vertex) with logarithmic scale axes.

7.1. Numerical solution near vertices

First, we compare the numerical solutions computed for $\alpha _k=10^{-6}$, $\tau =\{1.15\times 10^{-11},10^{-7},10^{4}\}$ and $\eta =\{0,0.5\}$ with the asymptotic solutions derived in § 6. The borehole radius $\alpha _{k}$ and time $\tau$ were selected to ensure that the numerical solutions are in the neighbourhood of the three vertices $E$, $I$ and $F$. The two values of $\eta$ bracket the range of poroelastic effects, which are inexistent if $\eta =0$ but maximum if $\eta =1/2$.

Figure 7 illustrates the fluid-pressure profile $\varPi$ and crack-aperture profile $\varOmega$ on the crack $\tilde {\mathcal {C}}$, at the three vertices. These quantities have been scaled by their borehole values. The scaled borehole radius is of order $O(0.1)$ at the $E$-vertex, but $\alpha _\tau \ll 1$ at the $I$- and $F$-vertices. In particular, $\alpha _\tau =0.37$ for $\eta =0$ and $\alpha _\tau =0.71$ for $\eta =1/2$ at the $E$-vertex, and, thus, the $E$-vertex solution in figure 7 starts from the middle of the plots. The poroelastic coefficient $\eta$ has hardly any effect on the $I$- and $F$-vertex scaled solutions shown in figure 7, as these solutions are proportional to a power of $(1-\eta )$, which is cancelled when divided by $\varOmega _{w}$ and $\varPi _{w}$. However, near the $E$-vertex, there is a large dependence of $\alpha _\tau$ on $\eta$, due to its influence on crack length ${{\bar {\varLambda }}}$.

Figure 7. Scaled (a) crack aperture $\varOmega (\xi )/\varOmega _{w}$ and (b) pressure $\varPi (\xi )/\varPi _{w}$ near the $E$-, $I$- and $F$-vertices with poroelastic coefficient $\eta =\{0,0.5\}$.

The pore pressure field is significantly different at the three vertices as confirmed by the contour plots of $\varPi (\xi ,\zeta )$ in figure 8; there is indeed a radial-flow regime at the $E$- and $I$-vertices and a fracture-flow regime at the $F$-vertex.

Figure 8. Pressure field $\varPi (\xi ,\zeta )$ for time $\tau =\{10^{-11},10^{-4},10^{6}\}$ with parameters $\alpha _k=10^{-6}$ and $\eta =0$. Borehole boundary and bi-wing crack are plotted with thick solid lines, although the scaled boreholes in (b,c) are too small to be observed. (a) $\tau =10^{-11}$, E-vertex, (b) $\tau =10^{-4}$, I-vertex and (c) $\tau =10^{6}$, F-vertex.

7.2. Numerical solutions

Numerical simulations have been conducted to track the evolution of the waterflooding process from fracture initiation to the large-time self-similar fracture-flow regime. The histories of $\varPi_{w1}$, $\varOmega _{w}$, $\varLambda$ and $\varPhi$, computed for $\alpha _k=\{10^{-3},10^{-1},10\}$ and $\eta =\{0,0.25,0.5\}$, are shown in figure 9. Asymptotic solutions for $\eta =0.25$ at the $E$-, $I$- and $F$-vertices show good agreement with the numerical results; see figure 9.

Figure 9. Evolution of (a) borehole pressure $\varPi_{w1}$, (b) crack aperture at inlet $\varOmega _{w}$, (c) scaled crack length $\varLambda$ and (d) flooding efficiency $\varPhi$ for $\alpha _k=\{10^{-3},10^{-1},10\}$ and $\eta =\{0,0.25,0.5\}$. Asymptotic solutions in $E$-, $I$- and $F$-vertices for $\eta =0.25$ are shown in dashed lines.

The numerical results confirm that the borehole pressure increases with time in the radial-flow regime (6.7) but decreases with time in the fracture-flow regime (6.32b); a peak pressure therefore exists at an intermediate time between the two asymptotic solutions. This peak pressure is not related to fracture initiation (borehole breakdown), but reflects the transition between the two flow regimes. The borehole pressure is not affected by the poroelastic coefficient $\eta$ in the radial-flow regime, but is amplified by $(1-\eta )^{-({1}/{4})}$ in the fracture-flow regime, as predicted by the asymptotic solutions.

The crack aperture at the borehole inlet is shown in figure 9(b). Asymptotic solutions for the $E$-vertex (6.10), $I$-vertex (6.22) and $F$-vertex (6.32a) are also drawn in this figure. Due to the pore pressure induced volumetric strain, the fracture aperture decreases as the poroelastic coefficient $\eta$ increases. The numerical results show that the $I$-vertex regime is inexistent for $\alpha _k=10$. The absence of the $I$-vertex regime is also confirmed by the solution trajectory in the $\{\varPhi ,\varLambda \}$ phase figure for $\alpha _k=10$, which can be seen to be deflected from the $I$-vertex in figure 2.

Figure 9(c) illustrates the increase of crack length $\varLambda$ with $\tau$, especially in the radial-flow regime, from the fracture initial value (6.9b) to the KGD crack edge as $\varLambda =1$. The solution $\varLambda (\tau /\alpha _k^{2})$ for the radial-flow regime shown in figure 3 matches the numerical result well when $\alpha _k$ is small. However, a larger error is found for $\alpha _k=10$, as expected from the deflection of the solution trajectory for large $\alpha _k$ from the radial-flow edge shown in figure 2. The evolution of flooding efficiency $\varPhi$ plotted in figure 9(d) illustrates the transition from radial flow ($\varPhi =0$) to fracture flow ($\varPhi =1$).

7.3. Peak pressure

The dependence of the peak pressure on $\alpha _k$ and $\eta$ is illustrated in figure 10(a), with computations conducted for $\alpha _k=\{10^{-4},10^{-3},\ldots ,10\}$ and $\eta =\{0,0.25,0.5\}$. The corresponding time $\tau _{peak}$ to reach the peak pressure is shown in figure 10(b). Figure 11(a) shows the variation of crack length $\varLambda _{{peak}}(\alpha _k;\eta )$, reached when the injection pressure is maximum, with $\alpha _k$ for $\eta =\{0,0.25,0.5\}$, while figure 11(b) shows the corresponding ratio $\ell _{{peak}}/a =\varLambda _{peak} \sqrt {\tau _{peak}}/\alpha _k$. The points $\{\varPhi ,\varLambda \}$ pertaining to the maximum injection pressure are also marked in figure 2.

Figure 10. (a) Peak borehole pressure $\varPi_{w}$ varies with $\alpha _k=\{10^{-4},10^{-3},\ldots ,10\}$ and $\eta =\{0,0.25,0.5\}$; (b) corresponding time $\tau$ to reach the peak pressure.

Figure 11. (a) Scaled crack length $\varLambda$ and (b) the ratio of crack length compared with borehole radius $\ell /a$ when the borehole pressure reaches peak value.