2.1. Governing equations for a PKN hydraulic fracture in a permeable medium

Since the solid medium is assumed to be homogeneous, the fracture will grow symmetrically about the injection point. We choose to locate the origin of the $x,y$ coordinate system at this point (see figure 1). The primary unknowns in a HF problem are the field variables comprising the fracture aperture $w(x,y,t)$ : $z=\pm {w}({x,y},t)/2$ and the fluid pressure ${p}_f({x,y},t)$ within the fracture (or the net pressure ${p}({x,y},t)={p}_f({x},t)-\sigma _0$ ), both functions of $x$ , $y$ and $t$ , and the fracture half-length $\ell (t)$ , which, along with the constant fracture height $H$ , define the evolving fracture geometry $\mathcal{F}(t)= [-\ell (t),\ell (t)]\times [-H/2,H/2]$ . The solution depends on the volumetric injection rate $Q(t)$ and four alternate material parameters defined in order to keep subsequent formulae uncluttered, namely: the plane strain modulus $E'$ , the alternate viscosity ${\mu} '$ , the alternate fracture toughness $K'$ , and alternate Carter leak-off coefficient $C'$ , which are defined as follows:

(2.1) \begin{equation} E^{\prime }=\frac {E}{1-\nu ^{2}},\text{ }{\mu} ^{\prime }=12{\mu} ,\text{ }K^{\prime }=\bigg( \frac {32}{\unicode{x03C0} }\bigg) ^{1/2}K_{Ic},\text{ }C^{\prime }=2C_{L}. \end{equation}

2.1.1. Elasticity

The equations relating the aperture $w$ to the net pressure $p$ on the walls of a planar fracture occupying the rectangular region $(x,y)\in [-\ell ,\ell ]\times [-H/2,H/2]$ in an infinite, homogeneous, elastic solid can be condensed into a hypersingular integral equation of the form (see Appendix A.2 and Crouch & Starfield Reference Crouch and Starfield1983; Hills et al. Reference Hills, Kelly, Dai and Korsunsky1996)

(2.2) \begin{equation} p(x,y,t) =-\frac {E^{\prime }}{8\unicode{x03C0} }\int _{-\ell }^{\ell }\int _{-H/2}^{H/2}\frac {w(x^{\prime },y^{\prime },t)\textrm{d}x^{\prime }\textrm{d}y^{\prime }}{\left [ (x^{\prime }-x)^{2}+(y^{\prime }-y)^{2}\right ] ^{3/2}} .\end{equation}

If the aspect ratio of the rectangular fracture is large, i.e. $\ell \gg H$ , the typical assumption that is made in the PKN approximation (Perkins & Kern Reference Perkins and Kern1961; Nordgren Reference Nordgren1972) is that a state of plane strain prevails for vertical cross-sections. In a state of plane strain $w$ is independent of $x^{\prime }$ and (2.2), in the limit $\ell \rightarrow \infty$ , is reduced (see A.3.1 of Appendix A) to that of a plane strain fracture in the vertical direction with semi-height $H/2$ . Now, under the assumption that the pressure in such vertical cross-sections is constant, (Adachi & Peirce Reference Adachi and Peirce2008) used the elliptic shape of the Sneddon solution for a plane strain fracture subject to a constant pressure (Sneddon Reference Sneddon1995) (see also A.4) to motivate the following ansatz for the fracture aperture $w(x,y,t)$ in which the $x$ and $y$ variables are assumed to be separated into the following product:

(2.3) \begin{equation} w(x,y,t)=w(x,0,t)\sqrt {1-(2y/H)^{2}} .\end{equation}

Now applying the vertical averaging operator $\mathcal{A}(\boldsymbol{\cdot }):=\frac {1}{H}\int _{-H/2}^{H/2}\,\textrm{d}y$ to (2.3) and defining $\bar {w}:=\mathcal{A}(w)$ we obtain

(2.4) \begin{equation} w(x,0,t)=\frac {4}{\unicode{x03C0} }\bar {w}(x,t) .\end{equation}

Substituting the ansatz (2.3) into (2.2) and evaluating the inner iterated integral explicitly, makes it possible to reduce the planar integral equation to a 1-D integral equation (Adachi & Peirce Reference Adachi and Peirce2008) as follows:

(2.5) \begin{align} p(x,t) &=-\frac {E^{\prime }}{2\unicode{x03C0} ^{2}}\int _{-\ell }^{\ell }\bar {w}(x^{\prime },t)\left [ \int _{-H/2}^{H/2}\frac {\sqrt {1-(2y^{\prime }/H)^{2}}\textrm{d}y^{\prime }}{\left [ (x^{\prime }-x)^{2}+y^{\prime }{}^{2}\right ] ^{3/2}}\right ] \textrm{d}x^{\prime } ,\\[-8pt]\nonumber \end{align}

(2.6) \begin{align} &=\frac {\bar {E}}{\unicode{x03C0} H}\int _{-\ell }^{\ell }\bar {w}(x^{\prime },t)\frac {{\textrm d}G(2(x^{\prime }-x)/H)}{{\textrm d}x^{\prime }}\textrm{d}x^{\prime } , \end{align}

where

(2.7) \begin{align} \bar {E}=({2}/{\unicode{x03C0} })E^{\prime }, G(s)=\frac{\sqrt{1+s^{2}}}{s}E \!\left(\frac{1}{1+s^{2}} \right) \end{align}

and

(2.8) \begin{align} E(m)=\int\limits_{0}^{\unicode{x03C0} /2}\sqrt {1-m\sin ^{2}\theta }{\textrm d}\theta \end{align}

is the complete elliptic integral of the second kind. This non-local approximation to the PKN geometry has also been extended (Dontsov & Peirce Reference Dontsov and Peirce2015a ) to include slender geometries such as the so-called P3D geometry (Adachi et al. Reference Adachi, Detournay and Peirce2010). We observe that the kernel $G$ of the integral equation (2.6) has the following asymptotic behaviour (Adachi & Peirce Reference Adachi and Peirce2008):

(2.9) \begin{align} &G(\rho )\overset {\rho \rightarrow 0}{=}\left ( \frac {1}{\rho }-\frac {1}{2}\rho \ln \rho +O(\rho )\right ) ,\ G(\rho )\overset {\rho \gg 1}{\rightarrow }\frac {\unicode{x03C0} }{2}\left ( \textrm{sign}(\rho )+\frac {1}{4\rho ^{2}}+O(\rho ^{-4})\right ) \text{ and }\nonumber \\&\qquad\qquad G^{\prime }(\rho )\overset {\rho \gg 1}{\rightarrow }\unicode{x03C0} \delta (\rho ) .\end{align}

Using the scaling property of the delta function $\delta (cs)= ({\delta (s)}/{| c|})$ and the asymptotic behaviour (2.9) $_c$ , we observe that, for points $|x-\ell |\gg H/2$ away from the tip region, which we define as $|x-\ell |\lt {H}/{2}$ , the non-local elasticity equation (2.6) reduces to the classic local pressure-width equation

(2.10) \begin{equation} p=\bar {E}\frac {\bar {w}}{H} .\end{equation}

We will revisit this formal calculation in § 3.2.2, in which we consider the asymptotic behaviour of the pressure within and away from the tip region.

2.1.2. Lubrication

The Reynolds lubrication equation is obtained by combining Poiseuille’s law and the continuity equation to yield

(2.11) \begin{equation} \frac {\partial w}{\partial t}=\frac {1}{{\mu} ^{\prime }}\boldsymbol{\nabla }\boldsymbol{\cdot }( w^{3}\boldsymbol{\nabla }p) -g(x,y,t)+Q(t)\delta (x,y) .\end{equation}

We assume that the velocity of fluid leak-off to the porous medium $g(x,y,t)$ is captured by Carter’s leak-off model $g(x,y,t)= ({C^{\prime }}/{\sqrt {t-t_{0}(x,y)}})$ , where $t_{0}(x,y)$ denotes the time of first exposure of point $(x,y)$ to the fracturing fluid (Carter Reference Carter, Howard and Fast1957). The point source is represented by the $\delta$ -function in (2.11) and we will be considering a propagation phase during which the fluid is injected at a constant volumetric rate $Q_0$ followed by a shut-in phase initiated at time $t_s$ , after which time there is no further fluid injected into the fracture. Thus the source function $Q(t)$ can be expressed as

(2.12) \begin{align} Q(t)=\begin{cases} Q_{0} & 0\lt t\lt t_{s}\\ 0 & t\;{\geqslant}\; t_{s}. \end{cases} \end{align}

Applying the vertical averaging operator $\mathcal{A}(\boldsymbol{\cdot })$ to (2.11) we obtain

(2.13) \begin{equation} \frac {\partial \bar {w}}{\partial t}=\frac {1}{\bar {{\mu} }}\frac {\partial }{\partial x}\left ( \bar {w}^{3}\frac {\partial p}{\partial x}\right ) -g(x,t)+Q(t)\frac {\delta (x)}{H} ,\end{equation}

where $\bar {{\mu} }=({\unicode{x03C0} ^{2}}/{12}){\mu} ^{\prime }=\unicode{x03C0} ^{2}{\mu}$ and $g(x,t)= {C^{\prime }}/({\sqrt {t-t_{0}(x)}})$ . We note that, in contrast to (2.11), the $\delta$ -function here implicitly represents a uniformly distributed line source over the full height of the channel, i.e. covering the interval $ y\in [-H/2,H/2]$ , but with an area injection rate $Q(t)/H$ scaled to be equivalent to the volumetric injection rate of the point source given in (2.11).

Combining the averaged lubrication equation (2.13) with the local elasticity equation (2.10), we obtain the classical porous medium partial differential equation (PDE) that governs the evolution of the PKN fracture (Kemp Reference Kemp1989; Kovalyshen & Detournay Reference Kovalyshen and Detournay2010)

(2.14) \begin{align} \frac {\partial \bar {w}}{\partial t}=\frac {\bar {E}}{4\bar {{\mu} }{H}}\frac {\partial ^{2}\bar {w}^{4}}{\partial x^{2}}-g(x,t)+Q(t)\frac {\delta (x)}{H}. \end{align}

2.1.3. Initial, boundary and propagation conditions at the moving front $\ell (t)$

Initial conditions: the initial conditions are formally given by

(2.15) \begin{equation} \ell =0,\quad \bar {w}=0,\quad p=0,\quad \text{at}\quad t=0\text{.} \end{equation}

Boundary conditions: for coalescent fluid and fracture fronts the boundary conditions at the crack tips $x=\pm \ell$ are given by zero fracture aperture and zero flux conditions (Detournay & Peirce Reference Detournay and Peirce2014)

(2.16) \begin{equation} \bar {w}=0,\quad \bar {w}^{3}\frac {\partial p}{\partial x}=0,\quad \text{at}\quad x=\pm \ell . \end{equation}

For a symmetric fracture the $\delta$ -function source can be replaced by an equivalent flux boundary condition at $x=0$

(2.17) \begin{equation} \frac {1}{\bar {{\mu} }} \bar {w}^{3}\frac {\partial p}{\partial x}(0^{+},t)=\frac {Q(t)}{2H}, \end{equation}

so that only half the fracture needs to be modelled.

Given these initial and boundary conditions, the lubrication equation, assuming Carter leak-off (2.13), can, after integrating both in time and space and exploiting symmetry, be expressed alternatively as the following global continuity equation:

(2.18) \begin{equation} 2\int _{0}^{\ell (t)}\bar {w}(x,t)\,\textrm{d}x+4C^{\prime }\int _{0}^{\ell (t)}\sqrt {t-t_{o}(x)}\,\textrm{d}x=V_{f}(t), \end{equation}

where

(2.19) \begin{align} V_{f}=\begin{cases} Q_{o}t/H & 0\lt t\lt t_{s}\\ V_{o}=Q_{o}t_{s}/H & t\;{\geqslant}\; t_{s}. \end{cases} \end{align}

This equation simply establishes that the total volume of fluid injected at time $t$ is equal to the volume of fluid contained in the crack plus the total volume of fluid lost to the permeable rock.

Propagation condition: since the fracture is assumed to propagate in limit equilibrium, LEFM (Rice 1968; Spence & Sharp Reference Spence and Sharp1985) implies that the alternate mode-I stress intensity factor $K$ satisfies the following inequality constraint:

(2.20) \begin{equation} K:=\bar {E}\lim _{x\rightarrow \pm \ell }\frac {\bar {w}(x)}{\sqrt {( \ell \mp x) }}\;{\leqslant}\; \bar {K}, \end{equation}

where $\bar {K}= ({8}/{\unicode{x03C0} })^{1/2} K_{Ic}$ . Equality in (2.20) occurs when the fracture is propagating, i.e. $V\gt 0$ , while strict inequality in (2.20) is associated with arrest characterised by $V=0, \ 0 \;{\leqslant}\; K \lt \bar {K}$ , and recession is characterised by $V\lt 0$ and $K=0$ .

We note that in the PKN formulation we have introduced the barred quantities $\bar {E}, \ \bar {{\mu} }$ , and $ \bar {K}$ to keep the PKN equations uncluttered by numerical factors. For convenience, the alternate PKN parameters are summarised

(2.21) \begin{equation} \bar {E}=\frac {2}{\unicode{x03C0} }E^{\prime }=\frac {2}{\unicode{x03C0} }\frac {E}{1-\nu ^{2}}, \text{ }\bar {{\mu} }=\frac {\unicode{x03C0} ^{2}}{12}{\mu} ^{\prime }=\unicode{x03C0} ^{2}{\mu} , \text{ }\ \bar {K}=\left (\frac {8}{\unicode{x03C0} }\right )^{1/2} K_{Ic} .\end{equation}

3.1. Stretched coordinate system

In the analysis that follows, it is convenient to introduce a stretched coordinate system $s(t)=x/\ell (t)$ . The unknown functions in this coordinate system will be represented by $\bar {w}({s},t)$ and ${p}({s},t)$ . It is also convenient to introduce a coordinate $\hat {x}$ located at the tip and pointing inwards toward the centre of the HF (see figure 1). The corresponding stretched coordinate located at the tip $s=1$ is represented by $\hat {s}=1-s$ and the unknown functions in this case are $\hat {w}(\hat {s},t)$ and $\hat {p}(\hat {s},t)$ .

After integration by parts, the elasticity equation (2.6), in the tip-based stretched coordinate system $\hat {s}$ , can be written in the form

(3.1) \begin{equation} \hat {p}(\hat {s},t)=-\frac {\bar {E}}{\unicode{x03C0} H}\int \limits _0^2\frac {{\textrm d}\hat {w}(\hat {s}^{\prime },t)}{{\textrm d}\hat {s}^{\prime }}G\left(\frac {\hat {s}^{\prime }-\hat {s}}{\varepsilon }\right){\textrm d}\hat {s}^{\prime }=I_{\mathcal{A}}+I_{\mathcal{B}} ,\end{equation}

where $\varepsilon =H/(2\ell ) \ll 1$ . For the purposes of the subsequent asymptotic analysis, it is convenient to decompose the domain of integration in (3.1) into the union of two disjoint subsets as follows: $(0,2)=\mathcal{A}\cup \mathcal{B}$ . The set $\mathcal{A}$ is used to denote the set of points $\hat {s}$ remote from the source point $\hat {s}^{\prime }$ , i.e. $| \hat {s}^{\prime }-\hat {s}| \gt \varepsilon$ , and the contribution from this region to the integral on the right of (3.1) is denoted by $I_{\mathcal{A}}$ . The set $\mathcal{B}$ is used to denote the set of points within an $\varepsilon$ -neighbourhood of the source point, i.e. $| \hat {s}^{\prime }-\hat {s}| \lt \varepsilon$ , and the contribution from this region to the integral on the right of (3.1) is denoted by $I_{\mathcal{B}}$ .

The lubrication equation (2.13) in the stretched tip coordinate system becomes

(3.2) \begin{equation} \frac {\partial \hat {w}}{\partial t}+(1-\hat {s})\frac {\dot {\ell }}{\ell }\frac {\partial \hat {w}}{\partial \hat {s}}=\frac {1}{\bar {{\mu} }\ell ^{2}}\frac {\partial }{\partial \hat {s}}\left ( \hat {w}^{3}\frac {\partial \hat {p}}{\partial \hat {s}}\right ) -\hat {g}(\hat {s},t)+Q(t)\frac {\delta (\hat {s}-1)}{H} ,\end{equation}

where $\hat {g}(\hat {s},t) = {C^{\prime }}/({\sqrt {t-t_{0}(\ell (1-\hat {s}))}})$ and $0 \;{\leqslant}\; \hat {s}\;{\leqslant}\; 1$ .

3.2. Asymptotic analysis of the elasticity operator

In this subsection we investigate the impact that the distinct far- and near-field asymptotic behaviours of the elasticity kernel $G(s)$ , provided in (2.9), have on the action of the non-local integral operator when the receiving point $\hat {s}$ is remote from as well as near to the tip region. Although the aperture and pressure are functions of both $\hat {s}$ and $t$ , for brevity, in this subsection we will suppress the explicit dependence on $t$ .

3.2.1. Asymptotic expansions for the elasticity operator

Assuming $\mathcal{A}=\{ \hat {s}^{\prime }\in (0,2):| \hat {s}^{\prime }-\hat {s}| \gt \varepsilon \}$ and $\mathcal{B}=\{ \hat {s}^{\prime }\in (0,2):| \hat {s}^{\prime }-\hat {s}| \lt \varepsilon\}$ , we decompose the crack region $(0,2)=\mathcal{A}\cup \mathcal{B}$ into the union of a set $\mathcal{A}$ within which the far-field expansion (2.9) $_b$ holds

(3.3) \begin{equation} I_{\mathcal{A}}\sim -\frac {\bar {E}}{2H}\int \limits _{\mathcal{A}}\frac {{\textrm d}\hat {w}(\hat {s}^{\prime })}{{\textrm d}\hat {s}^{\prime }}\bigg[ \textrm{sign}\bigg(\frac {\hat {s}^{\prime }-\hat {s}}{\varepsilon }\bigg)+\frac {\varepsilon ^{2}}{4(\hat {s}^{\prime }-\hat {s})^{2}}+\ldots \bigg] {\textrm d}\hat {s}^{\prime } ,\end{equation}

and a set $\mathcal{B}$ within which the near-field expansion (2.9) $_a$ holds

(3.4) \begin{equation} I_{\mathcal{B}}\sim -\frac {\bar {E}}{\unicode{x03C0} H}\mathop {\int \!\!\!\!\!\!-}\limits _{\mathcal{B}}\frac {{\textrm d}\hat {w}(\hat {s}^{\prime })}{{\textrm d}\hat {s}^{\prime }}\bigg[ \frac {\varepsilon }{\hat {s}^{\prime }-\hat {s}}-\frac {1}{2}\frac {(\hat {s}^{\prime }-\hat {s})}{\varepsilon }\ln \frac {| \hat {s}^{\prime }-\hat {s}| }{\varepsilon }+\ldots \bigg] {\textrm d}\hat {s}^{\prime } ,\end{equation}

where the $\int \!\!\!\!\!-$ symbol implies that the integral should be interpreted in a Cauchy principal value sense, which is necessary because of the Cauchy behaviour of the leading term in the asymptotic expansion of the kernel.

3.2.2. Outer expansion

To investigate the asymptotic behaviour of $\hat {p}$ for receiving points $\hat {s}$ away from the tip region, i.e. $\hat {s}\gt \varepsilon$ , we choose $\mathcal{A}=(0,\hat {s}-\varepsilon )\cup (\hat {s}+\varepsilon ,2)$ and $\mathcal{B}=(\hat {s}-\varepsilon ,\hat {s}+\varepsilon )$ . Integrating the first term in (3.3) directly and evaluating the subsequent terms in the expansion as in (Adachi & Peirce Reference Adachi and Peirce2008), we obtain

(3.5) \begin{eqnarray} I_{\mathcal{A}} &\sim &\frac {\bar {E}}{H}\left \{ \frac {1}{2}\left [ \hat {w}( \hat {s}-\varepsilon )+\hat {w}(\hat {s}+\varepsilon )\right ] -\int \limits _{ \mathcal{A}}\frac {{\textrm d}\hat {w}(\hat {s}^{\prime })}{{\textrm d}\hat {s}^{\prime }}\left [ \frac {\varepsilon ^{2}}{8(\hat {s}^{\prime }-\hat {s})^{2}}+\ldots \right ] {\textrm d} \hat {s}^{\prime }\right \}\nonumber \\ &\sim &\frac {\bar {E}}{H}\bigg\{ \hat {w}(\hat {s})-\frac {1}{4}\varepsilon ^{2}\ln \varepsilon \frac {{\textrm d}^{2}\hat {w}(\hat {s})}{{\textrm d}\hat {s}^{2}}+\ldots \bigg\}, \end{eqnarray}

and the leading behaviour of (3.4) can be shown to be (Adachi & Peirce Reference Adachi and Peirce2008)

(3.6) \begin{align} I_{\mathcal{B}}\sim -\frac {\bar {E}}{\unicode{x03C0} H}\bigg\{ \varepsilon ^{2}\bigg( \frac {41}{18}+\frac {2}{3}\ln 2\bigg) \frac {{\textrm d}^{2}\hat {w}(\hat {s})}{{\textrm d}\hat {s}^{2}}+\ldots \bigg\} .\end{align}

Combining these two asymptotic expansions we observe that the asymptotic expansion for the pressure at receiving points away from the tip region is dominated by the far-field behaviour within the set $\mathcal{A}$ and is given by

(3.7) \begin{equation} \hat {p}(\hat {s})=\frac {\bar {E}}{H}\bigg[ \hat {w}(\hat {s})-\frac {1}{4}\varepsilon ^{2}\ln \varepsilon \frac {d^{2}\hat {w}(\hat {s})}{{\textrm d}\hat {s}^{2}}+O(\varepsilon ^{2})\bigg] .\end{equation}

We note that the leading behaviour of this so-called outer expansion for the pressure is the same as that for the local elasticity operator (2.10), which was derived formally using the properties of the $\delta$ -function.

3.2.3. Inner expansion

For receiving points $\hat {s}$ within the tip region $\hat {s}\lt \varepsilon$ , we choose $\mathcal{A}=(\hat {s}+\varepsilon ,2)$ and $\mathcal{B}=(0,\hat {s}+\varepsilon )$ and assume that $\hat {w}(\hat {s})=\hat {A}(t)$ $\hat {s}^{\lambda }$ . From the analysis in (Adachi & Peirce Reference Adachi and Peirce2008), it follows that the leading behaviour of $I_{\mathcal{A}}$ is given by

(3.8) \begin{align} I_{\mathcal{A}}\sim O(\varepsilon ^{2}), \end{align}

and, defining the transformation $\hat {\xi }=\hat {s}/\varepsilon$ as in (Adachi & Peirce Reference Adachi and Peirce2008), the expression for $I_{\mathcal{B}}$ becomes

(3.9) \begin{eqnarray} I_{\mathcal{B}} &\sim &-\frac {\bar {E}\hat {A}\lambda \varepsilon ^{\lambda }}{ \unicode{x03C0} H}\mathop {\int \!\!\!\!\!\!-}\limits _{0}^{1+\hat {\xi }}{\hat {\xi }}^{\prime (\lambda -1)}\left [ \frac {1}{\hat {\xi }^{\prime }-\hat {\xi }}-\frac {1}{2}( \hat {\xi }^{\prime }-\hat {\xi }) \ln | \hat {\xi }^{\prime }-\hat {\xi }| +\ldots \right ] {\textrm d} \hat {\xi }^{\prime }\nonumber \\ &\sim &-\frac {\bar {E}\hat {A}\lambda \varepsilon ^{\lambda }}{\unicode{x03C0} H}\left \{\!\! \begin{array}{c} -\unicode{x03C0} \cot ( \unicode{x03C0} \lambda) \hat {\xi }^{\lambda -1}\quad \text{for }\dfrac { 1}{2}\;{\leqslant}\; \lambda \lt 1 \\[8pt] -\ln | \hat {\xi }|\quad \text{for }\lambda =1 \end{array} \right \} +C. \end{eqnarray}

Thus for receiving points within the tip region, the leading behaviour for the pressure corresponding to the aperture power law is of the form

(3.10) \begin{equation} \hat {p}(\hat {s})\sim \frac {\bar {E}\hat {A}}{2\ell }\left \{ \begin{array}{c} \lambda \cot \left ( \unicode{x03C0} \lambda \right ) \hat {s}^{\lambda -1}\quad \text{for } \dfrac {1}{2}\;{\leqslant}\; \lambda \lt 1 \\[8pt] \dfrac {1}{\unicode{x03C0} }\ln |\hat {s}|\quad \text{for }\lambda =1\end{array}\right \} +C .\end{equation}

We note that (3.10) is the PKN equivalent of the plane strain eigenfunction presented in appendix C of Spence & Sharp (Reference Spence and Sharp1985) and derived in Appendix A (see (A35)).

3.3. Vertex asymptotes for the coupled equations

Throughout this section we will assume that the leading behaviour for the aperture is a power law $\hat {w}(\hat {s},t)\stackrel {\hat {s}\rightarrow 0}{\sim } \hat {A}(t) \hat {s}^{\lambda }$ . Substituting this aperture power law and the corresponding leading behaviour for the tip pressure (3.10) into the tip lubrication equation (3.2) leads to the following asymptotic relations:

(3.11) \begin{align} \frac {1}{2} \;{\leqslant}\; \lambda \lt 1:\dot {\hat {A}}\hat {s}^{\lambda }+\frac {\dot {\ell }}{\ell }\lambda \hat {A}\hat {s}^{\lambda -1}&\sim \frac {\hat {A}^{4}\bar {E}}{\bar {{\mu} }\ell ^{3}}\lambda (\lambda -1)(2\lambda -1)\cot ( \unicode{x03C0} \lambda ) \hat {s}^{4\lambda -3}- \hat {g} ,\\[-9pt]\nonumber \end{align}

(3.12) \begin{align} \lambda =1:\dot {\hat {A}}\hat {s}+\frac {\dot {\ell }}{\ell }\hat {A}&\sim \frac {\hat {A}^{4}\bar {E}}{\unicode{x03C0} \bar {{\mu} }\ell ^{3}}\hat {s}-\hat {g} . \end{align}

We note that if $\lambda \gt 1$ the aperture power law and the corresponding tip pressure cannot satisfy the lubrication equation and the elasticity equation simultaneously. The plain strain state plays a fundamental role in the tip asymptotics of any planar fracture with a smooth boundary in a 3-D elastic medium. Indeed, in the limit as the receiving point approaches the fracture boundary (see Peirce & Detournay (Reference Peirce and Detournay2008) and Appendix A.2), the planar elasticity integral equation (A27) reduces to that of a semi-infinite fracture in a state plane strain (A33). Comparing (3.10) $_a$ with (A35), we observe that the non-local PKN asymptote (3.10) $_a$ can be made formally identical to the plane strain asymptote (A35), by replacing $\bar {E}$ in (3.10) by $E^{\prime }/2$ . Moreover, if we replace $\bar {{\mu} }$ in the PKN lubrication equation (2.13) by ${\mu} ^{\prime }$ then we obtain the plane strain lubrication equation – but for the difference in the source term, which is of no consequence in this asymptotic analysis. Thus by making the substitutions

(3.13) \begin{equation} \bar {E} \rightarrow \frac {E^{\prime }}{2}\quad \text{ and }\quad \bar {{\mu} }\rightarrow {\mu} ^{\prime }, \end{equation}

the non-local PKN asymptotic relations (3.11)–(3.12) become precisely those for a HF in a state of plane strain (see for example (Peirce & Detournay Reference Peirce and Detournay2022c )). Conversely, to transition from the plane strain asymptotic relations to those of the non-local PKN model, the following substitutions are required:

(3.14) \begin{equation} E^{\prime } \rightarrow 2 \bar {E} \quad\text{ and } {\mu} ^{\prime } \rightarrow \bar {{\mu} } .\end{equation}

The ability to transition from the plane strain asymptotic relations to those of the non-local PKN fracture, makes it possible to use the established plane strain propagation asymptotes (Detournay Reference Detournay2004, Dontsov & Peirce, Reference Dontsov and Peirce2015b ; Detournay Reference Detournay2016) to determine the corresponding PKN asymptotes by making the substitutions (3.14). Not only do these substitutions make it possible to obtain the so-called vertex asymptotes for the PKN from those for plane strain, but they can also be used to determine the corresponding PKN edge and tri-process multiscale asymptotes. Indeed, the multiscale asymptotes (Garagash et al. Reference Garagash, Detournay and Adachi2011; Dontsov & Peirce Reference Dontsov and Peirce2015b ), established for the plane strain case, can be used directly to obtain the PKN asymptotes by making use of the substitutions (3.14).

3.3.1. Propagation asymptotes ( $K=\bar {K},V\gt 0$ )

Viscous asymptote: in the absence of leak-off and in an elastic medium with close to zero toughness, the dominant balance occurs between the second and third terms in (3.11). Matching these terms, we find $\lambda =2/3$ and, solving for $\hat {A}$ , we obtain the viscous asymptote

(3.15) \begin{equation} \bar {w}_{m}(x)=\bar {\beta }_{m}\left ( \frac {\bar {{\mu} }V}{\bar {E}}\right ) ^{1/3}(\ell -x)^{2/3} ,\end{equation}

where $\bar {\beta }_{m}=3^{5/6}$ . Note that applying the substitutions (3.14) to the plane strain viscous asymptote $w_{m}(x)=\beta _{m} ({{\mu} ^{\prime }V}/{E^{\prime }} )^{1/3}(\ell -x)^{2/3}$ , where ${\beta }_{m}=2^{1/3} 3^{5/6}$ (see Spence & Sharp Reference Spence and Sharp1985; Desroches et al. Reference Desroches, Detournay, Lenoach, Papanastasiou, Pearson, Thiercelin and Cheng1994), we recover (3.15).

Now if we combine the expression for the average width $\bar {w}$ defined in (2.4) with the viscous asymptote $w_{m}$ , we find that $\bar {w}=({\unicode{x03C0} }/{4})w_{m}=({3}/{4})^{1/3}\bar {w}_{m}(x)$ , which amounts to a $9\,\%$ difference between these two estimates for the average asymptote $\bar {w}$ . This is not surprising considering the exchange of limits implied by that taking the asymptotic limit of the averaged equations to obtain $\bar {w}_{m}$ , compared with the latter evaluation of $\bar {w}$ by applying the averaging factor $({\unicode{x03C0} }/{4})$ to the plane strain viscous asymptote.

Leak-off asymptote: if leak-off is dominant, we replace $\hat {g}$ by the asymptotic behaviour in the tip region $ \hat {g} \sim ({C^{\prime }\dot {\ell }^{1/2}})/{\ell ^{1/2}\hat {s}^{1/2}}$ , and find that the dominant balance is between the third and last terms in (3.11). It follows that $\lambda =5/8$ and, solving for $\hat {A}$ , we obtain

(3.16) \begin{equation} \bar {w}_{\tilde {m}}(x)=\beta _{\tilde {m}}\left ( \frac {\bar {{\mu} }C^{\prime }V^{1/2}}{\bar {E}}\right ) ^{1/4}(\ell -x)^{5/8} ,\end{equation}

where $\beta _{\tilde {m}}={4}/{[15 ( \sqrt {2}-1 ) ]^{1/4}}$ . Note that if we apply the substitutions (3.14) to the classic Lenoach tip asymptote $w_{\tilde {m}}(x)=\beta _{\tilde {m}} ({2{\mu} ^{\prime }C^{\prime }V^{1/2}}/{E^{\prime }} )^{1/4}(\ell -x)^{5/8}$ see (Lenoach Reference Lenoach1995), we obtain (3.16). If we combine the expression for the average width $\bar {w}$ defined in (2.4) and the Lenoach tip asymptote $w_{\tilde {m}}$ , we find that $\bar {w}=({\unicode{x03C0} }/{4})w_{\tilde {m}}=({3\unicode{x03C0} }/{8})^{1/4}\bar {w}_{\tilde {m}}(x)$ , which amounts to a $4\,\%$ difference between these two estimates for the average asymptote $\bar {w}$ .

Toughness asymptote: the toughness asymptote follows from the propagation condition (2.20)

(3.17) \begin{equation} \bar {w}_{\bar {k}}(x)=\frac {\bar {K}}{\bar {E}} (\ell -x)^{1/2} .\end{equation}

Note that if we apply the substitutions (3.14) to the plane strain tip asymptote $w_{k'}(x)= ({K^{\prime }}/{E^{\prime }})(\ell -x)^{1/2}$ , we obtain (3.17). Alternatively, if we combine the expression for the average width $\bar {w}$ defined in (2.4) and the plane strain tip asymptote $w_{k'}(x)$ , we recover the toughness asymptote (3.17), i.e. $\bar {w}=({\unicode{x03C0} }/{4}) w_{k'}(x)= \bar {w}_{\bar {k}}(x)$ .

3.3.2. Arrest asymptotes ( $0\;{\leqslant}\; K\lt \bar {K},\ V=0$ )

The arrest phase is characterised by $V=0$ while the stress intensity factor $K$ steadily decreases from $\bar {K}$ to $0$ as fluid continues to leak from the fracture into the porous medium. During propagation $\hat {g}$ is singular, however, after the fracture has come to rest $\hat {g}$ is no longer singular, and, in fact, becomes more spatially uniform in the tip region as time progresses. In this case, it is appropriate to assume that $\hat {g}(\hat {s},t)\sim \hat {g}_0(t)$ .

Arrest- $k$ -asymptote ( $K\gt 0,V = 0$ ): during arrest ( $V=0$ ) so the near-tip asymptote is thus given by LEFM

(3.18) \begin{equation} \bar {w}_k=\frac {K}{\bar {E}}(\ell -x)^{1/2}, \end{equation}

where $K$ is a decreasing function of time. We refer to $\bar {w}_k$ as the $k$ -asymptote.

Arrest-recession transition- $g$ -asymptote ( $K=0,V=0$ ): at the end of the arrest phase, both $K$ and $V$ vanish. Although the passage from arrest to recession is ephemeral, the transition is accompanied by a switch in the near-field asymptote. To construct the intermediate asymptote that applies at this transition point, we proceed as above and assume a power law for the aperture, replace $\hat {g}(\hat {s},t)$ in (3.11) by $\hat {g}_{0}(t),$ and match the third and fourth terms in (3.11) to find $\lambda =({3}/{4})$ and, solving for $\hat {A}$ , we obtain

(3.19) \begin{equation} \bar {w}_{g}=\beta _{\bar {g}}\left (\frac {\bar {{\mu} }\hat {g}_{0}}{\bar {E}}\right )^{1/4}(\ell -x)^{3/4} ,\end{equation}

where $\beta _{\bar {g}}= ({32}/{3}) ^{1/4}\simeq 1.807$ and we refer to $\bar {w}_{g}$ as the $\bar {g}$ -asymptote. We note that making the substitutions (3.14) for the $g$ -asymptote derived in (Peirce & Detournay Reference Peirce and Detournay2022b ) we recover (3.19).

3.3.3. The recession asymptote- $r$ -asymptote ( $K=0,V\lt 0$ )

If $V\lt 0$ , the $g$ -asymptote $\lambda =3/4$ is no longer admissible as the second term in (3.11) would become infinite and violate the dominant balance. Thus the only possibility is that $\lambda =1,$ and matching the second and fourth terms in (3.12), we obtain the $r$ -asymptote

(3.20) \begin{equation} \bar {w}=-\hat {g}_0(t)\frac {\ell }{\dot {\ell }}\hat {s}\quad \text{ or }\quad \bar {w}_{r}(x)=\hat {g}_0(t)\frac {(\ell -x)}{| V| } .\end{equation}

We observe that the first and third terms in (3.12) match at the next order.

3.4. The connection problem to establish multiscale solutions

Vertex and multiscale tip asymptotic solutions (Garagash et al. Reference Garagash, Detournay and Adachi2011; Dontsov & Peirce Reference Dontsov and Peirce2015b ; Detournay Reference Detournay2016) have proven to be extremely useful in the development of numerical algorithms to model propagating fractures (Adachi & Detournay Reference Adachi and Detournay2008; Peirce & Detournay Reference Peirce and Detournay2008; Lecampion et al. 2012; Peirce Reference Peirce2015, Reference Peirce2016; Dontsov & Peirce Reference Dontsov and Peirce2017), and, more recently, for the development of models for receding HFs in which fluid loss to the porous elastic medium causes the fracture to deflate and close (Peirce Reference Peirce2022; Peirce & Detournay Reference Peirce and Detournay2022b , Reference Peirce and Detournayc ). These algorithms can achieve solutions with a high degree of precision, able to account for the vertex and multiscale behaviour at the finest length scale, while using a relatively coarse mesh. The fundamental idea for propagating and receding HFs alike (Peirce & Detournay Reference Peirce and Detournay2008; Peirce Reference Peirce2015; Dontsov & Peirce Reference Dontsov and Peirce2017; Peirce Reference Peirce2022; Peirce & Detournay Reference Peirce and Detournay2022c ) is to use a trial value for the fracture aperture $\hat {w}$ , sampled at computational points within and adjacent to the fracture front, along with the applicable tip asymptote to estimate the fracture front location (or equivalently the local front velocity) at the current step. The front position for propagating or receding fractures is then adjusted iteratively and the aperture updated until both are consistent with the applicable asymptote. The efficacy of this approach was clearly demonstrated in a collaborative study that evaluated a number of numerical algorithms (Lecampion et al. 2012) and benchmarked their relative accuracy against an analytic solution and compared their relative efficiency.

The non-singular approach introduced in (Dontsov & Peirce Reference Dontsov and Peirce2015b ) provides an effective way to determine these multiscale asymptotes, namely: the viscous-toughness-leak-off multi-scale asymptotes for propagating fractures. (Peirce & Detournay Reference Peirce and Detournay2022b ) used the non-singular approach to determine the $k{-}g$ edge solution to capture the deflation of an arrested HF while the stress intensity factor decays to zero and the $r{-}g$ edge solution to capture the acceleration of the receding fracture from the arrested state $\dot {\ell }= 0$ till the $r$ -asymptote fully occupies the tip region.

The starting point for the non-singular approach (Dontsov & Peirce Reference Dontsov and Peirce2015b ) is the governing integral equation over the tip region that relates the fluid pressure to the fracture aperture. For the PKN fracture this integral equation is obtained by integrating (2.6) by parts, transforming to the tip coordinate $\hat {x}$ , zooming into the tip region, and making use of the leading-order expansion for $G$ , to yield

(3.21) \begin{equation} p(\hat {x},t) \sim -\frac {\bar {E}}{2\unicode{x03C0} }\int _{0}^{\infty }\frac {{\textrm d}\bar {w}(\hat {x}^{\prime })}{{\textrm d}\hat {x}^{\prime }}\frac {{\textrm d}\hat {x}^{\prime }}{\hat {x}^{\prime }-\hat {x}} .\end{equation}

The corresponding lubrication equation for the PKN fracture that expresses fluid conservation in the tip region is

(3.22) \begin{equation} \frac {\bar {w}^{3}}{\bar {{\mu} }}\frac {{\textrm d}p}{{\textrm d}\hat {x}}=V\bar {w}+\int \limits _{0}^{\hat {x}}\hat {g}(\hat {s}){\textrm d}\hat {s} ,\end{equation}

where the boundary conditions (2.16) have been used and

(3.23) \begin{equation} \hat {g}(\hat {x})=\left \{\!\!\begin{array}{l@{\quad}l} \dfrac {C^{\prime }V^{1/2}}{\hat {x}^{1/2}} & \text{for }V\gt 0 \\[8pt] \hat {g}_{0} & \text{for }V\;{\leqslant}\; 0.\end{array}\right . \end{equation}

We observe that the substitution (3.13) reduces (3.21) to the plane strain equation (A33) derived in Appendix A, while (3.22) is rendered formally identical to the corresponding lubrication equation for propagating and deflating fractures analysed in Dontsov & Peirce, (Reference Dontsov and Peirce2015b ) and Peirce & Detournay, (Reference Peirce and Detournay2022b ), respectively. The consequence of this formal equivalence is profound as the plane strain multiscale asymptotes can be used directly provided the appropriate transformations are made.

4.1. Propagation, shut-in and arrest

We consider a HF that propagates under the continuous injection of a viscous fluid until the shut-in time $t_{s}=$ $1500$ $\textrm{s}$ , after which injection is discontinued and the injection conduit is capped. The field-like parameters are: $E = 9.5$ $\textrm{GPa}$ , $\nu = 0.2$ , ${\mu} = 0.1$ $\textrm{Pa} \,\textrm{s}$ , $C = 8.2431\times 10^{-6}$ $\textrm{ms}^{-1/2},$ $K_{IC}= 1.0$ $\textrm{MPa}\,\textrm{m}^{1/2},$ $Q_{0} = 10^{-3}$ $\textrm{m}^{3}\,\textrm{s}^{-1}$ and $H= 20$ $\textrm{m}$ .

This parameter set corresponds to the dimensionless shut-in time ${\omega }\sim 2.0\times 10^{-3}$ and arrest regime parameter $\bar {\phi }^V\sim 2.49\times 10^{2}$ . The dimensionless shut-in time is defined to be ${\omega } = {t_{s}}/{\bar {t}_{m\tilde {m}}}$ , where $t_{s}$ is the actual shut-in time and $\bar {t}_{m\tilde {m}}$ is the time scale at which propagation transitions from being viscosity dominated to being dominated by leak-off. The arrest regime parameter is defined to be the ratio of two time scales $\bar {\phi }^{V}=({\bar {t}_{mk}^{V}})/({\bar {t}_{m\tilde {m}}^{V})}$ , where $\bar {t}_{mk}^{V}$ is the time (assuming a fixed injected volume, which is appropriate post shut-in) at which propagation transitions from viscosity to toughness dominated propagation, while $\bar {t}_{m\tilde {m}}^{V}$ again represents the transition time between viscosity dominated and leak-off dominated propagation – but now assuming a fixed injected volume. The $V$ superscript indicates fixed volume propagation.

If the time of arrest $t_{a}\sim \bar {t}_{mk}^{V}\ll \bar {t}_{m\tilde {m}}^{V}$ then $\phi ^{V}\ll 1$ , and the fracture is propagating in the toughness regime at the time of arrest but still contains a significant amount of fluid that needs to leak-off before recession can start. Conversely, if $t_{a}\sim \bar {t}_{m\tilde {m}}^{V}\ll t_{mk}^{V}$ then $\bar {\phi }^{V}\gg 1$ , and the fracture is propagating in the leak-off regime at the time of arrest and has lost sufficient fluid to preclude significant further propagation of the fracture. In this case recession will start almost immediately after a very short arrest period. See § 5 for more detail.

Figure 2. (a) The pre- and post-shut-in ILSA-PL3D fracture footprints (dashed red) and the corresponding IMMA-PKN fracture front positions (solid black), (b) the pre-shut-in fracture apertures that correspond to the fracture footprints in the top figure, (c) the post-shut-in fracture apertures that correspond to the fracture footprints in the top figure.

In the top plot in figure 2, we compare the first quadrant projections of the pre- and post-shut-in fracture footprints obtained using the ILSA-PL3D and IMMA-PKN algorithms. In the middle plot we compare the corresponding pre-shut-in fracture apertures, and in the bottom plot we compare the post-shut-in fracture apertures that correspond to the remaining fracture footprints shown in the top figure. As the fracture approaches arrest, the density of the footprints increases and the successive fracture apertures become much closer together. The PKN and PL3D solutions show very close agreement. The slight discrepancy for the first fracture footprint shown is to be expected as the initially radial planar fracture has just touched the upper boundary at which point the aspect ratio is close to unity so one would not expect the PKN solution to be particularly close to the PL3D solution.

Figure 3. (a) the ILSA-PL3D (dashed red) and the IMMA-PKN (solid black) aspect ratios as a function of time, (b) the efficiencies $\eta$ as functions of time, (c) the fracture apertures at the injection point as functions of time, (d) the fluid pressures at the injection point as functions of time.

In figure 3, we compare time evolution plots of the ILSA-PL3D (dashed red) and IMMA-PKN (solid black) solutions for the: aspect ratio $2 \ell /H$ , the fracture aperture at the injection point, the fluid pressure at the injection point and the efficiency, which is defined to be

(4.1) \begin{equation} \eta = 2 \int _{0}^{\ell (t)}\bar {w}(x,t)\,\textrm{d}x /V_f, \end{equation}

The shut-in time is indicated on the aspect ratio plot by the solid black square , which corresponds to the points at which the curvatures in the other plots in this figure change abruptly. As was the case with the solutions provided in figure 2, the ILSA-PL3D and IMMA-PKN solutions show close agreement. The only noticeable discrepancy is in the fluid pressure at the injection point near the start of injection, which is again due to the transition from a radial fracture (for the ILSA-PL3D algorithm) to a rectangular fracture during the initial stages of growth of the fracture.

The close agreement between the two sets of results, even in this environment in which multiple physical processes compete to determine the evolution of the fracture boundary, demonstrates the utility of the PKN asymptotes established by the analysis in § 3. Having established the efficacy of the PKN asymptotes in determining the post-shut-in dynamics of the HF propagation and the clear identification of the time of arrest, the IMMA-PKN algorithm is used to confirm the scaling analysis performed in § 5.

4.2. Propagation, shut-in, arrest and recession

We now provide IMMA-PKN and ILSA-PL3D solutions for a HF that evolves within a rectangular region of constant height $H$ and growing length $2\ell$ . In the case of the ILSA-PL3D algorithm the HF is constrained to propagate between two sedimentary layers by introducing a tenfold jump in the fracture toughness $K_{Ic}$ across layer interfaces located at $y=\pm H/2$ . Naturally, the rectangular footprint is a fundamental assumption of the PKN model. After shut-in, when the injection of fluid is discontinued, the fracture continues to propagate but gradually slows until it arrests when the stress intensity factor at the leading front of the fracture drops below the fracture toughness. This deceleration is not only due to the spreading of the fracture as the pressure gradients reduce but also due to the leak-off of fluid to the porous elastic medium. After arrest, the fracture continues to deflate due to fluid leak-off until the stress intensity factor is reduced to zero. Thereafter, with further loss of fluid due to leak off to the porous medium, the fracture tip starts to recede as the open parts of the fracture resume contact. This recession process is modelled by the ILSA-PL3D algorithm by the introduction of a minimum width constraint upon closure, while the IMMA-PKN algorithm uses the recession asymptote.

Rather than a selection of a set of parameters that might be encountered in the field, the following parameter set was chosen in order be able to accurately model the rectangular fracture evolution and recession process using the ILSA-PL3D algorithm: $E = 0.0096$ $\textrm{GPa}$ , $\nu = 0.2$ , ${\mu} = 300$ $\textrm{Pa} \,\textrm{s}$ , $C = 1.9758 \times 10^{-3}$ $\textrm{ms}^{-1/2},$ $K_{IC}= 0.01$ $\textrm{MPa}\,\textrm{m}^{1/2},$ $t_{s}=$ $10^{8}$ $\textrm{s},$ $Q_{0} = 1$ $\textrm{m}^{3}\textrm{s}\,^{-1}$ and $H= 216.8614$ $\textrm{m}$ . This choice of elastic moduli and fracture height $H$ leads to large fracture apertures so the minimum width is set to $w_{\textit{min}}=10^{-8}$ to avoid convergence issues that a smaller width constraint would entail due to round-off errors. Containment to a rectangular footprint is also ensured by the appropriate choice of layered fracture toughness (Mori et al. Reference Mori, Peruzzo, Garagash and Lecampion2024). This parameter set corresponds to the dimensionless shut-in time $\omega \sim 5.9\times 10^{-3}$ and arrest regime parameter $\bar {\phi }^V\sim 1.9571\times 10^{3}$ (see § 5 for the definition and interpretation of these parameters).

In the top plot in figure 4, we compare the first quadrant projections of the ILSA-PL3D (dashed red) and IMMA-PKN (solid black) pre-arrest footprints. Because both $x$ and $y$ axes are scaled to the factor $H/2$ , the sequence of footprints represent the evolving aspect ratio $2\ell /H$ . The aspect ratio at shut-in is approximately $2\ell _s/H\sim 11.44$ , while after shut-in the fracture continues to propagate until arrest when the aspect ratio is approximately $2\ell _a/H\sim 13.74$ . In the middle plot in figure 4, we compare the receding fracture fronts of the ILSA-PL3D (dashed blue) and IMMA-PKN (solid black) algorithms. In contrast to the nearly vertical fronts of the propagating ILSA-PL3D fracture footprints shown in the top plot, the receding ILSA-PL3D fracture fonts become more curved as recession progresses. These receding ILSA-PL3D fracture footprints eventually separate from the original containment rectangle when $2\ell /H\sim 8.5$ and ultimately assume a more elliptic shape. While the IMMA-PKN fracture fronts initially follow the maximum points of the corresponding fracture fronts for the ILSA-PL3D model, they ultimately approximate an averaged front position such that the fracture volume matches that of the fully planar model. The red vertical line in the middle plot when $2\ell /H\sim 4.5$ represents the last recession front of the PKN model before collapse. In the numerical IMMA-PKN simulations, an efficiency threshold $\eta \lt 0.5\,\%$ was used as an objective criterion to identify the collapse time. It is typically not possible for the IMMA-PKN model to remain stable beyond the efficiency threshold.

Figure 4. (a) pre- and post-shut-in IMMA-PKN (solid black) and ILSA-PL3D (dashed red) fracture footprints before arrest. (b) receding IMMA-PKN (solid black) and ILSA-PL3D (dashed blue) fracture footprints up till the point of collapse of the PKN solution, (c) ILSA-PL3D (solid blue) fracture footprints subsequent to the PKN collapse.

In the bottom plot in figure 4, we render the ILSA-PL3D fracture footprints (solid blue) as well as the last recession front of the PKN model, which is the same red vertical line shown in the middle figure. Since the width constraint logic used in the ILSA-PL3D model has been carefully calibrated, we assume that it provides an accurate capture of the dynamics of the receding fracture. Given that the PKN model is constrained to a rectangular geometry, and cannot faithfully capture the significant topological changes of the ILSA-PL3D model once it separates from the containment layer boundaries, it is to be expected that the two models deviate beyond this point. Indeed, the PKN front location tries to capture an average position, whose location is set to match the PKN fracture volume to that of the fully curved fracture ILSA-PL3D model. Since the PKN model significantly over-estimates the fracture footprint, to maintain volume balance it is forced to under-estimate the fracture aperture, which, in turn, leads to an earlier estimate of the collapse time. In contrast, the PL3D model has the freedom to model the non-rectangular post-separation fracture footprints and is, therefore, able to capture the dynamics of the HF much closer to the point of collapse.

We see that these PL3D recession fronts, which persist beyond the collapse of the PKN solution, approach an approximately elliptic shape. For the last five such elliptical fracture footprints, we also represent by solid red curves the ellipses whose semi-major axes correspond to the fracture lengths and whose semi-minor axes correspond to the fracture heights. We observe that these red ellipses are almost indistinguishable from the PL3D fracture footprints, which are represented by the solid blue lines.

In the top plot in figure 5, we compare the pre-shut-in fracture apertures for the IMMA-PKN (solid black) and ILSA-PL3D (dashed red) solutions. In the middle plot, we compare the post-shut-in fracture apertures representing the propagation up to the point of arrest. In the bottom plot, the dashed magenta curve represents the ILSA-PL3D solution during the brief period during which the fracture continues to deflate while arrested, while the corresponding IMMA-PKN solution is represented by a solid black curve. The dashed blue curves represent the ILSA-PL3D fracture apertures while the corresponding IMMA-PKN fracture apertures are represented by the solid black curves. The last of these receding fracture aperture comparisons is sampled just before the collapse of the IMMA-PKN model. We observe that the IMMA-PKN and ILSA-PL3D fracture apertures represented in figure 5 show close agreement.

Figure 5. (a) pre-shut-in IMMA-PKN (solid black) and ILSA-PL3D (dashed red) scaled fracture apertures $w/w_s$ plotted as functions of the scaled distance $2 x/H$ from the injection point. Here $w_s=w(0,t_s)$ is the aperture at the well-bore at the time of shut-in $t_s$ . (b) post-shut-in scaled fracture apertures. (c) scaled receding fracture apertures for the IMMA-PKN (solid black) and ILSA-PL3D (dashed blue) algorithms sampled up to the point of collapse of the PKN model.

In the top left plot in figure 6, we compare the scaled fracture dimensions for the IMMA-PKN and ILSA-PL3D solutions as functions of the scaled time $t/t_s$ . The horizontal dimension is determined by the scaled fracture length $\ell /\ell _s$ , which is represented by the solid black curve for the PKN model and the dashed black curve for the PL3D model. On the $\ell /\ell _s$ vs $t/t_s$ plot, the shut-in time is represented by the solid black square , the arrest point is represented by the solid black circle , while the point of initiation of recession is represented by the solid black triangle $\blacktriangledown$ . The vertical dimension is determined by the scaled fracture height $h/H$ , which is represented by the solid red curve for the PKN model and the dashed red curve for the PL3D model. The separation time $t/t_s\sim 4.5$ at which the dashed red curve deviates from the solid red curve, represents the time at which the PL3D fracture footprint first loses contact with the confining layers $y=\pm H/2$ . The fracture efficiency $\eta$ as a function of time is plotted in the top right figure. The scaled fracture aperture at the injection point $w/w_s$ is plotted as a function of time in the bottom left plot, while the scaled fluid pressure $p/p_s$ at the injection point is plotted as a function of time in the bottom right plot.

Figure 6. The IMMA-PKN solution is represented by (solid) curves and the ILSA-PL3D solution by (dashed) curves. Here $\ell _s$ , $w_s$ , $p_s$ , are the fracture half-length, and aperture and pressure at the injection point all sampled at the time of shut-in $t_s$ . Moreover, $h$ represents the vertical dimension of the fracture at the injection point: for the PKN model $h=H$ , while for the PL3D model $h\;{\leqslant}\; H$ . Before the PL3D fracture reaches the constraining layers $h\lt H$ briefly, followed by $h=H$ , and finally the PL3D fracture retreats from the constraining layers and $h\lt H$ once again (designated by the red dashes). (a) scaled fracture dimensions, (b) fracture efficiency $\eta$ , (c) scaled fracture aperture at the injection point $w/w_s$ , (d) scaled fluid pressure at the injection point $p/p_s$ , all plotted as functions of the scaled time $t/t_s$ .

5.1. Scaling for a PKN fracture subject to a constant injection rate $Q_{0}$

Following Detournay (Reference Detournay2004, Reference Detournay2016) we introduce a length scale ${l}_{\ast }$ , time scale $t_{\ast }$ , and characteristic aperture ${w}_{\ast }$ and net pressure ${p}_{\ast }$ such that $x={l}_{\ast }\xi$ , $t=t_{\ast }\tau$ , $\bar {w}={w}_{\ast }\varOmega$ , and $p={p}_{\ast }\varPi$ into the governing equations (2.6) and (2.13), enabling us to identify the following five dimensionless groups:

(5.1) \begin{equation} \mathcal{G}_{s}=\frac {{l}_{\ast }H{w}_{\ast }}{Q_{0}t_{\ast }},\text{ }\mathcal{G}_{m}=\frac {{l}_{\ast }\bar {{\mu} }Q_{0}}{{w}_{\ast }^{3}{p}_{\ast }H},\text{ }\mathcal{G}_{c}=\frac {C^{\prime }{l}_{\ast }H}{t_{\ast }^{1/2}Q_{0}},\text{ }\mathcal{G}_{e}=\frac {\bar {E}{w}_{\ast }}{H{p}_{\ast }},\text{ and }\mathcal{G}_{k}=\frac {\bar {K}}{\sqrt {H}{p}_{\ast }}. \end{equation}

The storage-viscosity $m$ -scaling can be identified by requiring $\mathcal{G}_{s}=\mathcal{G}_{m}=\mathcal{G}_{e}=1$ , from which it follows that the length ${l}_{m}$ , aperture ${w}_{m}$ and pressure ${p}_{m}$ scales are respectively given by

(5.2) \begin{equation} {l}_{m}=\left ( \frac {\bar {E}Q_{0}^{3}t^{4}}{\bar {{\mu} }H^{4}}\right ) ^{1/5}\quad {w}_{m}=\left ( \frac {\bar {{\mu} }Q_{0}^{2}t}{\bar {E}H}\right ) ^{1/5}\quad \text{and }\quad {p}_{m}=\left ( \frac {\bar {E}^{4}\bar {{\mu} }Q_{0}^{2}t}{H^6}\right ) ^{1/5}, \end{equation}

while the dimensionless toughness and leak-off coefficients become

(5.3) \begin{equation} \mathcal{G}_{k}:=\mathcal{K}_{m}(t)=\left ( \frac {\bar {t}_{km}}{t}\right ) ^{1/5}\quad \text{ and }\quad \mathcal{G}_{c}:=\mathcal{C}_{m}(t)=\left ( \frac {t}{\bar {t}_{m\tilde{m}}}\right ) ^{3/10} ,\end{equation}

where $\bar {t}_{km}= ({\bar {K}^{5}H^{7/2}})/({\bar {E}^{4}\bar {{\mu} }Q_{0}^{2}})$ and $\bar {t}_{m\tilde {m}}= ({\bar {{\mu} }Q_{0}^{2}}/{\bar {E}C^{\prime 5}H}) ^{2/3}$ are the $k{-}m$ and $m-\tilde {m}$ transition time scales, respectively.

We define the propagation mode parameter to be the ratio of these two characteristic time scales

(5.4) \begin{equation} \bar {\phi }=\frac {\bar {t}_{km}}{\bar {t}_{m\tilde {m}}}=\left ( \frac {\bar {K}^{6}H^{5}C^{\prime 4}}{\bar {{\mu} }^{2}\bar {E}^{4}Q_{0}^{4}}\right ) ^{5/6} .\end{equation}

The leak-off ( $\tilde {m}$ -scaling) can be obtained by requiring $\mathcal{G}_{c}=1$ instead of $\mathcal{G}_{s}$

(5.5) \begin{equation} {l}_{\tilde {m}}=\frac {Q_{0}t^{1/2}}{HC^{\prime }}\quad {w}_{\tilde {m}}=\left ( \frac {\bar {{\mu} }Q_{0}^{2}}{\bar {E}HC^{\prime }}\right ) ^{1/4}t^{1/8}\quad \text{ and }\quad {p}_{\tilde {m}}=\left ( \frac {\bar {E}^{3}\bar {{\mu} }Q_{0}^{2}}{H^{5}C^{\prime }}\right ) ^{1/4}t^{1/8}, \end{equation}

while the dimensionless toughness and storage coefficients become

(5.6) \begin{equation} \mathcal{G}_{k}:=\mathcal{K}_{\tilde {m}}(t)=\left ( \frac {\bar {t}_{k\tilde {m}}}{t}\right ) ^{1/8}\quad \text{ and }\quad \mathcal{G}_{s}:=\mathcal{S}_{\tilde {m}}(t)=\left ( \frac {\bar {t}_{m\tilde{m}}}{t}\right ) ^{3/8}, \end{equation}

where $\bar {t}_{k\tilde {m}}= ({\bar {K}^{8}H^{6}C^{\prime 2}})/({\bar {E}^{6}\bar {{\mu} }^{2}Q_{0}^{4}}).$ The toughness ( $k$ -scaling) can be obtained by requiring $\mathcal{G}_{k}=1$ instead of $\mathcal{G}_{m}$ , and the toughness-leak-off ( $\tilde {k}$ -scaling) can be obtained by requiring $\mathcal{G}_{k}=1$ and $\mathcal{G}_{c}=1$ see (Dontsov Reference Dontsov2022).

5.2. Characteristic power laws for arrest

In this paper we consider the injection of fluid at a constant flux $Q_{0}$ followed by shut-in at a certain time $t_{s}$ . Propagation of the HF, initiated at the beginning of the injection phase, may continue after shut-in depending on the regime of propagation. A propagating fracture in a permeable medium will ultimately come to rest either due to excessive leak-off if $\bar {K}=0$ or because the stress intensity factor $K$ has dropped below the fracture toughness when $\bar {K}\gt 0$ . There is an arrest period during which $K$ decreases as the fracture continues to lose fluid to the porous medium so that $\dot {w}\lt 0$ , until transition to the recession asymptote is initiated when $K=0$ and the fracture starts to recede. However, if the HF is already in the leak-off regime at the time of shut-in, then recession can be expected to start almost immediately after arrest.

The appropriate scaling for the dynamics with a fixed injected volume $V_{0}$ can be obtained Mori & Lecampion (Reference Mori and Lecampion2021) directly from those of a fracture driven by a constant flux $Q_{0}$ given in (5.2) by making the simple substitution $Q_{0}=V_{0}/t$ . In this case the length ${l}_{m}^{V}(t)$ and aperture ${w}_{m}^{V}(t)$ scales are given by

(5.7) \begin{equation} {l}_{m}^{V}(t)=\left ( \frac {\bar {E}V_{0}^{3}t}{\bar {{\mu} }H^{4}}\right ) ^{1/5}\quad \text{ and }\quad {w}_{m}^{V}(t)=\left ( \frac {\bar {{\mu} }V_{0}^{2}}{\bar {E}Ht}\right ) ^{1/5}, \end{equation}

while the dimensionless toughness $\mathcal{K}_{m}^{V}(t)$ and leak-off $\mathcal{C}_{m}^{V}(t)$ parameters become

(5.8) \begin{equation} \mathcal{K}_{m}^{V}(t)=\bigg( \frac {t}{\bar {t}_{mk}^{V}}\bigg) ^{1/5}\quad \text{ and }\quad \mathcal{C}_{m}^{V}(t)=\bigg( \frac {t}{\bar {t}_{m\tilde {m}}^{V}}\bigg) ^{7/10}, \end{equation}

where $\bar {t}_{mk}^{V}=({\bar {E}^{4}\bar {{\mu} }V_{0}^{2}})/({\bar {K}^{5}H^{7/2}})$ and $\bar {t}^V_{m\tilde {m}}= ( ({\bar {{\mu} }^{2}V_{0}^{4}})/({\bar {E}^{2}C^{\prime 10}H^{2}} )) ^{1/7}$ are the $m{-}k$ and $m-\tilde {m}$ transition time scales associated with the post-shut-in dynamics of a PKN HF, respectively.

In order to characterise the modes of arrest, we define the following arrest regime parameter $\bar {\phi }^{V}$ for a fixed injected volume PKN HF to be:

(5.9) \begin{equation} \bar {\phi }^{V}=\frac {\bar {t}_{mk}^{V}}{\bar {t}_{m\tilde {m}}^{V}}=\left ( \frac {\bar {{\mu} }^{2}\bar {E}^{12}C^{\prime 4}V_{0}^{4}}{\bar {K}^{14}H^{9}}\right ) ^{\frac {5}{14}}. \end{equation}

Once shut-in has occurred the two transition times defined in (5.8) identify two different ways the fracture behaves after the time of arrest $t_{a}$ , which, if $\bar {K}\gt 0$ , is characterised by the stress intensity factor dropping below the critical fracture toughness and, as a result, the velocity of the fracture going to zero. If $t_{a}\sim \bar {t}_{mk}^{V}\ll \bar {t}_{m\tilde {m}}^{V}$ , so that $\bar {\phi }^{V}\ll 1$ , then at the time of arrest, the fracture is propagating in the toughness regime and still contains a significant amount of fluid that needs to leak-off before recession can start. In this case, there will be a significant period during which the fracture deflates while it is in a state of arrest after which recession will begin. Conversely, if $t_{a}\sim \bar {t}_{m\tilde {m}}^{V}\ll t_{mk}^{V}$ , so that $\bar {\phi }^{V}\gg 1$ , then at the time of arrest the fracture is propagating in the leak-off regime and has lost sufficient fluid to preclude significant further propagation of the fracture. In this case, recession can be expected to start almost immediately and there will be a very short arrest period. We observe that the parameter $\bar {\phi }^{V}$ defined in (5.9) has no meaning in the zero toughness case since $\bar {t}_{mk}^{V}=\infty$ .

In the analysis that follows we consider a given shut-in time $t_{s}$ , which we define with respect to the constant injection rate storage-leak-off transition time $\bar {t}_{m\tilde {m}}$ in terms of the parameter $\omega$ defined by

(5.10) \begin{equation} {\omega } =\frac {t_{s}}{\bar {t}_{m\tilde {m}}}. \end{equation}

At shut-in, the injected volume is $V_{s}=Q_{0}t_{s}$ , the dimensionless leak-off coefficient $\mathcal{C}_{m}(t_{s}):=\mathcal{C}_{s}=\omega ^{3/10}$ , while the dimensionless toughness $\mathcal{K}_{m}(t_{s}):=\mathcal{K}_{s}=\left(\bar{\phi}\right)^{1/5}\omega ^{-1/5}$ . The fixed injected volume transition times, at which $\mathcal{K}_{m}^{V}(t)=1$ and $\mathcal{C}_{m}^{V}(t)=1$ , can be expressed in terms of $\mathcal{K}_{s}$ and $\mathcal{C}_{s}$ as follows:

(5.11) \begin{equation} \bar {t}_{mk}^{V}:=t_{s}\mathcal{K}_{s}^{-5}\quad \text{ and }\quad \bar {t}_{m\tilde {m}}^{V}:=t_{s}\mathcal{C}_{s}^{-10/7}. \end{equation}

Now, making use of the transition times defined in (5.3) and (5.8) as well as the definition (5.9) of the arrest regime parameter $\bar {\phi }^{V},$ the following relationship can be established between the fixed injected volume transition times $\bar {t}_{m\tilde {m}}^{V}$ and $\bar {t}_{mk}^{V}$ and the dimensionless shut-in parameter $\omega$ :

(5.12) \begin{equation} \bar {t}_{m\tilde {m}}^{V}=t_{s}\omega ^{-\frac {3}{7}}\quad \text{ and }\quad \bar {t}_{mk}^{V}=\bar {\phi }^{V}t_{s}\omega ^{-\frac {3}{7}} .\end{equation}

Naturally, because of the definition of the regime parameter $\bar {\phi }^{V}$ in (5.9), both the transition times $\bar {t}_{m\tilde {m}}^{V}$ and $\bar {t}_{mk}^{V}$ have the same power law dependence on $\omega$ . If $\omega \ll 1,$ then, at the time of shut-in $t_{s},$ the fracture is still propagating in the viscous-storage regime well away from the leak-off regime, so the arrest time $t_{a}$ is determined by the transition times $\bar {t}_{mk}^{V}$ if $\bar {K}\gt 0$ or by $\bar {t}_{m\tilde {m}}^{V}$ in the zero toughness case. The corresponding asymptotic dependence of the arrest length and arrest aperture can be obtained by substituting these transition times into (5.7). If $\omega \gg 1,$ then, at the time of shut-in $t_{s},$ the fracture is well into the leak-off regime so arrest occurs almost immediately, i.e. $t_{a}\sim t_{s}.$ Defining ${l}_{s}:=$ ${l}_{m}(t_{s})$ and ${w}_{s}:={w}_{m}(t_{s})$ , and combining (5.12) $_{b}$ with (5.7), we obtain

(5.13) \begin{align} t_{a}/t_{s} &\sim \left \{\!\! \begin{array}{c} \bar {\phi }^{V}\omega ^{-\frac {3}{7}} \text{ if }\omega \ll 1 \\[6pt] 1 \text{ if }\omega \gg 1\end{array}\right .\!\! ,\nonumber \\[-11pt] \end{align}

(5.14) \begin{align} {l}_{a}/{l}_{s} &\sim \left \{\!\! \begin{array}{c} ( \bar {\phi }^{V}) ^{1/5}\omega ^{-3/35} \text{ if }\omega \ll 1 \\[6pt] 1 \text{ if }\omega \gg 1\end{array}\right .\!\! ,\nonumber \\[-11pt] \end{align}

(5.15) \begin{align} {w}_{a}/{w}_{s} &\sim \left \{\!\! \begin{array}{c} ( \phi ^{V}) ^{-1/5}\omega ^{3/35} \text{ if }\omega \ll 1 \\[6pt] 1 \text{ if }\omega \gg 1\end{array}\right .\!\! .\nonumber \\[-11pt] \end{align}

For the zero toughness case $\bar {K}=0$ the arrest time $t_{a}$ will be determined by the transition time $\bar {t}_{m\tilde {m}}^{V},$ which has the same power law as $\bar {t}_{mk}^{V},$ without the factor involving the regime parameter $\bar {\phi }^{V}.$ Thus the power law asymptotes for the zero toughness case are identical to those in (5.13)–(5.15) without the factor involving the regime parameter $\bar {\phi }^{V}.$

The IMMA-PKN results for the non-zero toughness case $K^{\prime }\gt 0$ are presented in figure 7. The arrest time to shut-in time (a) and arrest length to shut-in length (b) ratios are plotted as functions of the dimensionless shut-in time $\omega$ for different values of the regime parameter $\bar {\phi }^V$ . Here the dashed red line represents the power law $A\omega ^{\alpha }$ for the particular case $\bar {\phi }^{V}=1957.1$ , with $A$ and $\alpha$ computed by a log linear regression of the first few data points. The regression yields $\alpha =-0.42$ for $t_a/t_s$ , which is consistent with the power law (5.13), and $\alpha =-0.0842$ for $\ell _a/\ell _s$ , again consistent with the power law (5.14). In both figures (a) and (b) the $\bullet$ symbol is used to indicate the curve corresponding to the regime parameter $\bar {\phi }^{V} =100$ , and the curve with the $\blacktriangle$ symbol corresponds to the regime parameter $\bar {\phi }^V=1957.1$ .

Figure 7. The solid black lines indicate the numerical solutions for (a) the arrest time to shut-in time ratios and (b) the arrest length to shut-in length ratios, both plotted as functions of $\omega$ for the values of the regime parameter $\bar {\phi }^{V} \in \{100\ ({\bullet }), 1957.1 ({\blacktriangle }) \}$ . The red symbols in each plot correspond to the parameter set $(\omega ,\bar {\phi }^{V})=(0.0059,1957.1)$ used in the PyFrac calibration runs in § 4.2. The dashed red lines represent linear regressions of the data set $({\blacktriangle })$ corresponding to $\bar {\phi }^{V}= 1957.1$ , assuming $t_a/t_s$ and $\ell _a/\ell _s$ are power laws of the form $A\omega ^{\alpha }$ .

Table 1. Scaled length $\ell (t)/\ell _s$ vs $t/t_s$ plotted for different $\bar {\phi }^V$ and $\omega$ values. In each of these plots: propagation under injection is represented by the black portion of the curve terminating with the symbol; the red portion of the curve, between the and the symbols, represents post shut-in propagation; the magenta portion of the curve, between the and the symbols, represent post shut-in propagation; the black portion of the curve, starting with the $\blacktriangledown$ symbol, represents recession.

5.3. Solution landscape represented by the evolving fracture length $\ell$

In table 1, the ratio of the fracture length scaled to the shut-in length $\ell (t)/\ell _s$ is plotted as a function of the scaled time $t/t_s$ for a range of different $\bar {\phi }^V$ and $\omega$ values. The value of the arrest regime parameter $\bar {\phi }^V$ is constant for each column of the table while the value of the dimensionless shut-in parameter $\omega$ is constant for each row. Fracture lengths are plotted for $\bar {\phi }^V= 10^{2}$ and $ 10^{3.29}\approx 1957.1$ and $\omega = 10^{-8}$ and $10^{-2}$ . The initially black parts of the curve represent the dynamics driven by injection up to the shut-in time $t_s$ , which is designated by a red asterisk . The subsequent red part of the curve, between the red shut-in point and the magenta arrest point , represents the length increase as the fracture continues to propagate between the time fluid injection ceases and the point of arrest. The magenta part of the curve, between the arrest point and the recession initiation point $\blacktriangledown$ , represents the period during which the fracture continues to lose fluid while at arrest, so that $V=0$ and $ {\partial w}/{\partial t}\lt 0$ . The subsequent black part of the curve, starting with $\blacktriangledown$ , represents the decreasing fracture length $\ell$ while the fracture recedes. Scanning down the first column $\bar {\phi }^V= 10^{2}$ , we observe that the collapse time occurs a factor of 1150 times later than the shut-in time for the case $\omega = 10^{-8}$ , which is reduced to 4.4 times for the case $\omega = 10^{-2}$ . Indeed, as $\omega$ increases, the post shut-in propagation time, arrest period and recession time till collapse all decrease.

6.1. Evidence of self-similarity

In order to motivate the analysis performed in this section we re-plot the data presented in § 4.2 in a form that reveals the self-similar nature of the solution as the fracture approaches the time of collapse $t_c$ . We recall the dimensionless parameters that characterise this set of data are given by the pair $(\omega ,\bar {\phi }^{V})=(0.0059,1957.1)$ .

Figure 8. In each of these plots the dashed red lines represent the lines obtained a linear regression on the first few points of each of the curves. (a) Log–log plot of the scaled fracture length $\ell /\ell _r$ for the PKN (solid black) and PL3D (solid blue) solutions and the scaled fracture height $h/H$ (solid magenta) vs the scaled reverse time $(t_c-t)/t_r$ . (b) The solid blue curve represents the PL3D scaled fracture aperture at the injection point $w/w_r$ vs the scaled reverse time $(t_c-t)/t_r$ .

In figure 8(a), we plot the scaled fracture length $\ell /\ell _r$ for the PKN (solid black) and PL3D (solid blue) solutions and the scaled fracture height $h/H$ for the PL3D solution (solid magenta) vs the scaled reverse time $(t_c-t)/t_r$ . Here $\ell _r$ and $t_r$ are the fracture length and time at which the fracture starts to recede, respectively. The linear trend in the asymptotic limit $t \rightarrow t_c$ indicates that $\ell$ behaves like a power law in this limit, i.e. $\ell \sim A (t_c-t)^\gamma$ . A linear regression on the first few points of each of these curves yields (represented by the dashed red lines) that for the PKN case $\gamma \approx 0.506$ and for the PL3D case $\gamma \approx 0.502$ . The exponent for the magenta curve representing the PL3D fracture height is $0.42$ , which is close to $0.5$ , indicating that the successive elliptic footprints are close to self-similar as the fracture recedes. In figure 8(b), we plot the PL3D scaled fracture aperture $w/w_r$ at the injection point (blue curve) as a function of the scaled reverse time $(t_c-t)/t_r$ . A plot of the PKN aperture is similar, except that the collapse time of the PKN fracture occurs earlier (and at a larger aperture value) than that of the PL3D fracture. To keep the figure uncluttered this plot is not provided. Due to the PKN model being constrained to a rectangular geometry, it is not able to capture the recession process once the closure boundary retreats from the containment layers. We observe that as $t \rightarrow t_c$ the fracture aperture tends a linear function of time, so that $w\sim \mathtt {g} \times (t_c-t)$ . A linear regression on the first few points of the PL3D curve shown in figure 8(b) yields the gradient $\mathtt {g}\approx 1.73 \times 10^{-7}\, \textrm{ms}^{-1}$ , while from a similar a linear regression of the corresponding PKN curve we obtain the estimate $\mathtt {g}\approx 2.54 \times 10^{-7}\, \textrm{ms}^{-1}$ .

Figure 9. (a) Scaled PL3D fracture aperture $w(x,t)/w(0,t)$ vs the scaled distance from the injection point $x/\ell (t)$ . (b) Scaled PKN fracture aperture $w(x,t)/w(0,t)$ vs the scaled distance from the injection point $x/\ell (t)$ . The dashed red curve represents the universal aperture to which the scaled numerical solutions tend as $t \rightarrow t_c$ .

In figure 9, for a sequence of sample times, starting with the recession initiation time $t_r$ , we plot the dynamically scaled fracture aperture $w(x,t)/w(0,t)$ against the scaled distance from the injection point $x/\ell (t)$ for: (a) the PL3D model (blue), and (b) the PKN model (black). We observe that in both cases the fracture apertures approach a universal function of the scaled distance $x/\ell (t)$ , which is represented by the dashed red curve. Since these dynamic scale factors are precisely the fracture length $\ell (t)$ for the abscissae and the injection aperture $w(0,t)$ for the ordinates, both of which exhibit power law behaviour according to figure 8, there is strong evidence of the existence of a similarity solution for both the PL3D model and the PKN model.

6.2. Reverse-time equations and the similarity ansatz

The evidence presented in § 6.1 indicates that the both the PL3D and PKN models for receding rectangular fractures ultimately tend to a self-similar structure. Though the PL3D equations are much more complex and not amenable to analysis, it is possible to analyse the PKN equations to investigate the nature of this self-similar behaviour, which is what we explore in the remainder of this section. Given the focus on the behaviour of the solution close to the closure time $t_c$ , to simplify the analysis we rewrite the governing equations in terms of the reverse time $t^{\prime }=t_{c}-t$ . In this case, $\dot {\ell }(t^{\prime })\gt 0$ for recession, the elasticity equation (3.1) remains unchanged, while the signs are changed in the last two terms in (3.2). As the HF approaches collapse, the current time $t$ is assumed to be sufficiently more advanced than the initiation times active in the collapsing fracture, i.e. $t\gg t_{0}(s,t)$ , that the leak-off term has the following asymptotic behaviour:

(6.1) \begin{equation} \hat {g}(\hat {s},t)=\frac {C^{\prime }}{\sqrt {t-t_{0}(\ell (1-\hat {s}))}}\overset {t\gg t_{0}}{\sim }\frac {C^{\prime }}{\sqrt {t}} .\end{equation}

Thus for sufficiently advanced times $t$ , the leak-off function $\hat {g}(\hat {s},t)$ is slowly varying in time and approximately independent of $\hat {s}$ , so that the leak-off function $\hat {g}$ may be considered to be approximately constant. Thus, in what follows, $\hat {g}$ is replaced by the constant $\hat {g}_{0}$ . The lubrication equation (3.2), expressed in terms of the reverse time $t^{\prime }$ , becomes

(6.2) \begin{equation} \frac {\partial \hat {w}}{\partial t^{\prime }}+(1-\hat {s})\frac {\dot {\ell }}{\ell }\frac {\partial \hat {w}}{\partial \hat {s}}=-\frac {1}{\bar {{\mu} }\ell ^{2}}\frac {\partial }{\partial \hat {s}}\left ( \hat {w}^{3}\frac {\partial \hat {p}}{\partial \hat {s}}\right ) +\hat {g}_{0} .\end{equation}

We observe that these reverse-time equations are equivalent to those of an inflating fracture subjected to a constant source distributed throughout its growing length.

Motivated by the evidence presented in figures 8 and 9, and the accompanying observations made in § 6.1, we look for a similarity solution to this ‘growing’ HF driven by an influx of fluid from a constant distributed source in terms of $s=x/\ell (t^{\prime })$ and by assuming a solution of the form

(6.3) \begin{equation} \bar {w}(s,t^{\prime })=t^{\prime \alpha }W(s)\quad p(s,t^{\prime })=t^{\prime \beta }P(s)\quad \ell (t^{\prime })=\varLambda t^{\prime \gamma } .\end{equation}

Note that for the reverse time, $t^{\prime }=0$ represents the time of closure.

6.3. Fracture volume

Define the average fracture aperture $\mathcal{W}$ to be

(6.4) \begin{equation} \mathcal{W}=\int \limits _{0}^{1}\bar {w}(s,t^{\prime }){\textrm d}s=t^{\prime \alpha }\int \limits _{0}^{1}W(s){\textrm d}s=t^{\prime \alpha }\overline {W}. \end{equation}

The crack volume $V_{c}$ can be expressed in terms of the similarity variables as follows:

(6.5) \begin{equation} V_{c}(t)=\mathcal{W}\ell =\varLambda t^{\prime {\alpha +\gamma }}\overline {W} .\end{equation}

Since the rate of change of volume should match the rate of efflux of fluid, it follows that

(6.6) \begin{equation} \dot {V_{c}}=\varLambda ( \alpha +\gamma ) t^{\prime {\alpha +\gamma -1}}\overline {W}=g_{0}\ell =g_{0}\varLambda t^{\prime \gamma } .\end{equation}

Matching powers we obtain the following:

(6.7) \begin{equation} \alpha =1\quad \text{ and }\quad \overline {W}=\frac {g_{0}}{1+\gamma }\text{ } .\end{equation}

6.4. Tip asymptotics and Taylor expansion

Now $\bar {w}(s,t^{\prime })=\hat {w}(\hat {s},t^{\prime })=t^{\prime }W(s)=t^{\prime }W(1-\hat {s})$ and, motivated by the linear asymptote in (3.20), we assume a Taylor expansion for $W$ about $s=1$ in powers of $\hat {s}$ of the form

(6.8) \begin{equation} \hat {w}(\hat {s},t^{\prime })=t^{\prime }\sum \limits _{n=1}^{\infty }\frac {(-1)^{n}w_{n}}{n!}\hat {s}^{n} ,\end{equation}

where $w_{n}= ({d^{n}W}/{{\textrm d}s^{n}})|_{s =1}$ and $w_0=W(1)=0$ .

Taking time and space derivatives of (6.8), the left side of (6.2) can be written in the form

(6.9) \begin{equation} \frac {\partial \hat {w}}{\partial t^{\prime }}+(1-\hat {s})\frac {\dot {\ell }}{\ell }\frac {\partial \hat {w}}{\partial \hat {s}}=-\gamma w_{1}+\sum \limits _{n=1}^{\infty }(-1)^{n}\frac {( 1-n\gamma ) w_{n}-\gamma w_{n+1}}{n!}\hat {s}^{n} .\end{equation}

Now because the asymptotic behaviour of the non-local kernel function $G((\hat{s}^{\prime }-\hat{s})/\varepsilon )$ changes depending on the distance between the sending point $\hat {s}^{\prime }$ and the receiving point $\hat {s},$ in order to determine the asymptotic behaviour of the pressure $\hat {p},$ it is necessary to split the action of integral operator into an inner region in which the receiving point is close to the tip $\hat {s}\lt \varepsilon$ , and an outer region for which the receiving point is remote from the tip $\hat {s}\gt \varepsilon .$

Inner region: using (3.10) to determine the leading behaviour for the pressure $\hat {p}$ from the action of the Cauchy operator on the leading term in (6.8), it follows that

(6.10) \begin{equation} \hat {p}\overset {\hat {s}\ll \varepsilon }{\sim }-\frac {\bar {E}t^{\prime {1-\gamma }}}{2\unicode{x03C0} \varLambda }w_{1}\ln \hat {s}. \end{equation}

Comparing (6.3) and (6.10) it follows that the time exponent of the pressure $\beta$ is given by

(6.11) \begin{equation} \beta =1-\gamma .\end{equation}

Using the expansion (6.10) for the pressure, the leading behaviour of the flux gradient in the tip region $\hat {s}\lt \varepsilon$ can be shown to be of the form

(6.12) \begin{equation} \frac {1}{\ell ^{2}}\frac {\partial }{\partial \hat {s}}\left ( \frac {\hat {w}^{3}}{\bar {{\mu} }}\frac {\partial \hat {p}}{\partial \hat {s}}\right ) \sim \frac {\bar {E}t^{\prime {4-3\gamma }}}{\unicode{x03C0} \bar {{\mu} }\varLambda ^{3}}w_{1}^{4}\hat {s} .\end{equation}

Outer region: away from the tip $\hat {s}\gt \varepsilon ,$ the leading behaviour for the pressure $\hat {p}$ is given by

(6.13) \begin{equation} \hat {p}(\hat {s},t)\sim \frac {\bar {E}}{H}\hat {w}(\hat {s},t) ,\end{equation}

so the leading behaviour of the flux gradient away from the tip can be shown to be of the form

(6.14) \begin{equation} \frac {1}{\ell ^{2}}\frac {\partial }{\partial \hat {s}}\left ( \frac {\hat {w}^{3}}{\bar {{\mu} }}\frac {\partial \hat {p}}{\partial \hat {s}}\right ) \sim \frac {\bar {E}t^{\prime {4-2\gamma }}}{\bar {{\mu} }\varLambda ^{2}H}3w_{1}^{4}\hat {s}^{2} .\end{equation}

Numerical evidence suggests that, as the fracture approaches closure it accelerates rather than slowing down, which implies that $\gamma \lt 1$ . Thus, in the small time limit $t^{\prime }\ll 1$ , it follows that $t^{\prime {4-3\gamma }}\ll 1$ and, a fortiori, $t^{\prime {4-2\gamma }}\ll 1,$ so we can neglect the flux term on the right of (6.2). The significance of the emerging subdominance of the pressure gradient terms in (6.2) in the limit $t^{\prime }\ll 1$ is profound as it signifies a decoupling of the elastic force balance equation from the fluid conservation equation, so that the remaining dominant terms in (6.2) essentially enforce a local kinematic condition between the fluid leak-off velocity and the rate of decrease of the aperture and length of the receding fracture.

Now matching the powers of $\hat {s}$ in (6.9) to the only non-zero term $\hat {g}_0$ on the right side of (6.2) (assuming that the flux terms are negligible on this time scale), we obtain the following value for $w_{1}$ and recursion for $w_{n},\ n \gt 1$ :

(6.15) \begin{equation} w_{1}=-\hat {g}_0/\gamma ,\text{ and }w_{n+1}=\frac {( 1-n\gamma ) }{\gamma }w_{n} .\end{equation}

From this recursion it follows that

(6.16) \begin{equation} w_{n}=-( 1-\gamma ) ( 1-2\gamma ) \ldots ( 1-( n-1) \gamma ) \hat {g}_0/\gamma ^{n} .\end{equation}

Substituting these values for $w_{n}$ into (6.8) we obtain the following expansion for $\hat {w}(\hat {s},t^{\prime })$ :

(6.17) \begin{equation} \hat {w}(\hat {s},t^{\prime })=\hat {g}_0t^{\prime }\left [ \frac {1}{\gamma }\hat {s}+\sum \limits _{n=2}^{\infty }\frac {(-1)^{n+1}( 1-\gamma ) ( 1-2\gamma ) \ldots ( 1-( n-1) \gamma) }{n!\gamma ^{n}}\hat {s}^{n}\right ] .\end{equation}

Note that the power series (6.17) provides a local solution centred on the fracture tip, which only converges for $\hat {s}\lt 1$ and will not yield a global solution valid at the centre of the fracture $\hat {s}=1$ . However, we observe from (6.17) that, for $\gamma \lt 1$ , there are a countable infinity of global solutions valid for $0\;{\leqslant}\; \hat {s}\;{\leqslant}\; 1$ , each corresponding to the reciprocals $\gamma =1/2,1/3,\ldots 1/N\ldots$ of the integers greater than or equal to $2$ , which terminate for finite $N$ to polynomial solutions. Indeed, for $\gamma =1/N,\text{ }N\;{\geqslant}\; 2$ the series solution (6.17) terminates after $N$ terms to yield the globally valid solution of degree $N$

(6.18) \begin{equation} \hat {w}(\hat {s},t^{\prime })=\hat {g}_0t^{\prime }(-1) ^{N}[ 1-( 1-\hat {s}) ^{N}] .\end{equation}

From (6.18) we observe that the solutions for $N$ odd are a-physical since, if this were the case $\hat {w}(\hat {x},t)\lt 0$ , so only solutions with even powers of $N$ are admissible. Furthermore, from (6.18) we also see that the spatial derivative for all these polynomial solutions satisfy the symmetry condition $({\partial \hat {w}}/{\partial \hat {s}})(1,t) = 0$ and have $N-1$ derivatives that vanish at the centre of the fracture $\hat {s}=1$ . The gradient at the fracture tip $({\partial \hat {w}}/{\partial \hat {s}})(0,t^{\prime }) = {\hat {g}_0t^{\prime } ( -1 ) ^{N}N}$ increases with $N$ so, of these even polynomials, the solution corresponding to $\gamma = 1/2$ is the last remaining admissible shape as the fracture approaches closure. Thus the polynomial associated with $\gamma = 1/2$ forms an attractor for the receding fracture solution, which we call the sunset solution.

6.5. The sunset solution

6.5.1. Closed form expression for the sunset solution

The similarity solution corresponding to $\gamma =1/2$ , when expressed as a function of $s$ , is given by

(6.19) \begin{equation} \bar {w}(s,t^{\prime })=\hat {g}_0t^{\prime }(1-s^{2}) ,\text{ } s=x/\ell ,\text{ }\ell =\varLambda t^{\prime {1/2}} .\end{equation}

It is interesting to note that slowly draining gravity currents through a layered porous medium also exhibit similarity solutions with a parabolic structure (Prichard et al. Reference Prichard, Woods and Hogg2001) .

The asymptotic behaviour of the pressure $p$ can be determined by substituting (6.19) into (2.6) and using the asymptotic behaviour of the kernel $G$ given in (2.9)

(6.20) \begin{equation} p_{\textit{PKN}}(s,t)=\left \{\!\! \begin{array}{l@{\quad}l} \dfrac {t^{\prime }\bar {E}\hat {g}_0}{H}( 1-s^{2}) & \text{ for } | 1-s|\gt \varepsilon \\[12pt] \dfrac {t^{\prime }\bar {E}}{2\unicode{x03C0} \varLambda t^{\prime 1/2}}\dfrac {\hat {g}_{0}}{| V| }\ln | 1-s| & \text{ for }| 1-s|\lt \varepsilon. \end{array} \right . \end{equation}

Since $\gamma =1/2$ for the sunset solution, it follows from (6.11) that the time exponent for the pressure $\beta =1-\gamma =1/2$ , which is consistent with (6.20) $_b$ , in which the similarity solution for the fracture half-length $\ell =\varLambda t^{\prime 1/2}$ has been used. In order to close the solution to the forward problem for a receding HF with initial conditions (2.15), the only remaining parameter to specify is the prefactor $\varLambda$ in the power law for $\ell$ defined in (6.3).

In figure 10, we plot the scaled fracture apertures $w/w_s$ against the scaled distance from the injection point $x/\ell _s$ for (a) the PL3D solution coloured blue and (b) the PKN solutions coloured black compared with the corresponding sunset solutions, which are designated by the dashed red curves. In both plots (a) and (b) the different symbols labelling each curve correspond to the times indicated by the corresponding symbols in figure 8(a).

Figure 10. Comparison of the receding (a) PL3D and (b) PKN solutions and the sunset solution for different sample times.

6.5.2. Estimation of the leak-off coefficient $C^{\prime }$

In the analysis to establish the sunset solution for plane strain and radial fractures (Peirce & Detournay Reference Peirce and Detournay2022a ), it was recognised that by determining the time derivative ${\partial \bar {w}}/{\partial t^{\prime }}$ at the injection point as the fracture approaches collapse, i.e. $t^{\prime }\rightarrow 0$ , we can obtain an estimate of $\hat {g}_0$ , which, according to (6.19), is given by the coefficient of the aperture profile that can be determined by linear regression of the aperture vs $t^{\prime }$ function. Thus as the fracture approaches collapse, $t\rightarrow t_{c}\gg t_{0}(\hat {x})$ and (6.1) yields the following estimate for the leak-off coefficient $C^{\prime }$ :

(6.21) \begin{equation} C^{\prime } \sim \hat {g}_0 \times \sqrt {t_{c}} .\end{equation}

Using (6.21) and the estimate of $\hat {g}_0$ obtained from the linear regression in figure 8(b), we obtain the estimate $C^{\prime }\approx 3.97 \times 10^{-3}$ $\textrm{ms}^{-1/2}$ , whereas the value actually used in the simulation was $C^{\prime }\approx 3.9516 \times 10^{-3}$ $\textrm{ms}^{-1/2}$ .