2.1 Governing equations

For such a planar axisymmetric pure opening mode fracture of radius $R(t)$ , quasi-static elasticity reduces to the following boundary integral equation between the fracture aperture $w(r,t)$ and the distribution of net loading $p(r,t)$ (Sneddon & Lowengrub Reference Sneddon and Lowengrub1969; Cleary & Wong Reference Cleary and Wong1985):

(2.1) $$\begin{eqnarray}\displaystyle w(r,t)=\frac{8R(t)}{\unicode[STIX]{x03C0}E^{\prime }}\int _{0}^{1}G(r/R(t),\unicode[STIX]{x1D709})p(\unicode[STIX]{x1D709}R(t),t)\,\text{d}\unicode[STIX]{x1D709}, & & \displaystyle\end{eqnarray}$$

where $E^{\prime }=E/(1-\unicode[STIX]{x1D708}^{2})$ is the plane-strain Young’s modulus, and $p(r,t)=p_{f}(r,t)-\unicode[STIX]{x1D70E}_{o}$ is the net pressure loading on the fracture faces, i.e. the fluid pressure in excess of the in situ compressive stress $\unicode[STIX]{x1D70E}_{o}$ normal to the fracture plane. The elastic kernel $G$ is recalled in appendix A for completeness. For such a penny-shaped fracture, the stress intensity factor can be expressed using the following weight function integral (Rice Reference Rice1972)

(2.2) $$\begin{eqnarray}\displaystyle K_{I}=\frac{2}{\sqrt{\unicode[STIX]{x03C0}R}}\int _{0}^{R}\frac{p(r,t)}{\sqrt{R^{2}-r^{2}}}r\,\text{d}r. & & \displaystyle\end{eqnarray}$$

The quasi-static fracture propagation condition stipulates that the mode I stress intensity factor $K_{I}$ is always equal to the material fracture toughness $K_{Ic}$ . In other words, $K_{I}=K_{Ic}$ in (2.2).

This penny-shaped fracture is driven by a point injection of a fluid at its centre. Under the assumption of an impermeable rock, for an incompressible fluid, the width-averaged fluid mass conservation in the fracture reduces to the following continuity equation in axisymmetry for $r\in [0,R(t)]$ :

(2.3) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}w}{\unicode[STIX]{x2202}t}+\frac{1}{r}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}(rwv)=0, & & \displaystyle\end{eqnarray}$$

where $v$ denotes the width-averaged radial fluid velocity and $q=wv$ is the corresponding fluid flux. The point-source injection condition at the fracture centre (neglecting the wellbore radius $a$ compared to the fracture size) is given by

(2.4) $$\begin{eqnarray}\displaystyle \lim _{r\rightarrow 0}2\unicode[STIX]{x03C0}r\,w(r,t)v(r,t)=Q_{o}. & & \displaystyle\end{eqnarray}$$

Correspondingly, at the fracture tip, the flow rate $q=wv$ as well as the fracture opening are zero (Detournay & Peirce Reference Detournay and Peirce2014),

(2.5a,b ) $$\begin{eqnarray}\displaystyle q(R(t),t)=0,\quad w(R(t),t)=0. & & \displaystyle\end{eqnarray}$$

The width-averaged fluid momentum equation for such a parallel plate axisymmetric flow is given by

(2.6) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70C}\left(\frac{\unicode[STIX]{x2202}v}{\unicode[STIX]{x2202}t}+v\frac{\unicode[STIX]{x2202}v}{\unicode[STIX]{x2202}r}\right)=-\frac{\unicode[STIX]{x2202}p_{f}}{\unicode[STIX]{x2202}r}-\frac{2}{w}\unicode[STIX]{x1D70F}_{w}, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D70C}$ is the fluid density and $\unicode[STIX]{x1D70F}_{w}$ is the average wall shear stress typically expressed as

(2.7) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70F}_{w}=f\left(Re,\frac{k_{R}}{w}\right)\times \frac{\unicode[STIX]{x1D70C}v|v|}{2}. & & \displaystyle\end{eqnarray}$$

The fluid velocity has to be understood as a local time average in the presence of turbulence. The averaged wall shear stress is a function of the Fanning friction factor $f$ which depends on the local Reynolds number $Re$ , the relative roughness of the flow geometry $k_{R}/w$ and the shape of the cross-section perpendicular to the flow direction. In our case, it is natural to define the Reynolds number as

(2.8) $$\begin{eqnarray}\displaystyle Re=\frac{\unicode[STIX]{x1D70C}wv}{\unicode[STIX]{x1D707}}. & & \displaystyle\end{eqnarray}$$

It is important to note that the Reynolds number at a given point varies inside the fracture as both the fracture width and fluid velocity vary in space and time. The Reynolds number (2.8) notably decreases from the fracture inlet to the fracture tip as the fracture width $w$ tends to zero near the tip. As a result, a laminar region is always expected to be present near the tip, albeit reducing to a boundary layer. For flow in a thin channel/fracture, it can be shown that the inertial terms on the left-hand side of the balance of momentum (2.6) scale as $R_{e}\times w/R$ (e.g. Szeri Reference Szeri2010) which remains typically small for a fracture as $w/R\ll 1$ . For a liquid, inertial terms are thus negligible compared to the viscous terms on the right-hand side (more precisely, they only play a role at extremely early times, see Garagash (Reference Garagash2006) and Zia & Lecampion (Reference Zia and Lecampion2017) for more discussion). In the following, we neglect inertial terms, such that the width-averaged balance of momentum reduces to the following relation between pressure gradient, fluid velocity and friction factor:

(2.9) $$\begin{eqnarray}\displaystyle v\times |v|=-\frac{w}{\unicode[STIX]{x1D70C}\times f\left(Re,{\displaystyle \frac{k_{R}}{w}}\right)}\frac{\unicode[STIX]{x2202}p_{f}}{\unicode[STIX]{x2202}r}. & & \displaystyle\end{eqnarray}$$

2.2 Friction factor with and without drag reducing additives

The evolution of the friction factor as a function of the Reynolds number and relative roughness $f(Re,k_{R}/w)$ is well known for circular pipes (see figure 2). A large number of phenomenological relations for the evolution of the friction factor have been proposed. All these phenomenological expressions are based on steady-state experimental data of pressure drop as a function of velocity. They are thus valid only when the time–space average of inertial terms is indeed negligible.

The transition to turbulence as well as the effect of pipe roughness on the wall shear stress can be clearly seen on Moody and Prandtl–Kármán type plot (figure 2). The effect of the addition of heavy molecular weight polymers has also been extensively studied due to its drastic effect on turbulence (see figure 2). The long polymer chain induces viscoelasticity, reducing turbulence significantly (see Virk (Reference Virk1975) and White & Mungal (Reference White and Mungal2008) and references therein). The effect saturates at relatively low polymer concentration: a maximum drag reduction asymptote (MDR) emerges in turbulent flow. Such a maximum drag reduction is therefore always obtained at a reasonable cost in industrial hydraulic fracturing operations.

2.2.1 Friction in fractures

In order to apply the phenomenological relations originally derived for pipe flow to the parallel plate type flow, Jones (Reference Jones1976) suggests obtaining the characteristic dimension $D_{eq}$ used to determine the Reynolds number for parallel plate flow by matching the laminar friction (which is given analytically by Poiseuille law for a parallel plate flow) with the expression for a circular pipe in the laminar flow regime $f_{laminar}=16\,Re_{D_{eq}}^{-1}$ . This ensures that the proper friction is recovered in the laminar regime for a parallel plate flow. The predictions of friction in the turbulent regime using $Re_{D_{eq}}$ in the expressions of friction for a circular pipe agrees very well with experimental data for rectangular ducts of different aspect ratios (Jones Reference Jones1976). In the case of parallel plate type flow, in the laminar regime, the mean velocity is given by Poiseuille law (see e.g. Lamb Reference Lamb1932) $v=-(w^{2}/12\unicode[STIX]{x1D707})(\unicode[STIX]{x2202}p/\unicode[STIX]{x2202}r)$ . Writing $f_{laminar}=16\times Re_{D_{eq}}^{-1}=16\unicode[STIX]{x1D707}/(\unicode[STIX]{x1D70C}D_{eq}v)$ in the laminar regime, the laminar equivalent pipe diameter is simply obtained by comparing (2.9) with the previous equation as $D_{eq}=(4/3)w$ . We therefore use $D_{eq}$ in the definition of the Reynolds number entering the expressions for the friction factor for circular pipes. In other words, following Jones (Reference Jones1976), we use the following definition for the equivalent laminar pipe Reynolds number $Re_{D_{eq}}$ :

(2.10) $$\begin{eqnarray}\displaystyle Re_{D_{eq}}={\textstyle \frac{4}{3}}Re, & & \displaystyle\end{eqnarray}$$

where $Re$ corresponds to the Reynolds number (2.8) defined with the mean fracture width as the characteristic flow dimension. Such a choice allows us to directly use the different friction models derived for circular pipe while ensuring that we recover the proper expression in the laminar flow regime. As can be seen in figure 2, the transition to turbulent flow starts around $Re_{D}\approx 2200$ in a pipe, which gives for parallel plate flow a critical number ${\mathcal{R}}_{c}\approx 2200\times 3/4\approx 1650$ – a value only slightly lower than the one obtained experimentally for a plane Poiseuille flow in a rectangular duct (Jones Reference Jones1976).

Note that we can also rewrite the lubrication law (2.9), following Dontsov (Reference Dontsov2016a )

(2.11) $$\begin{eqnarray}\displaystyle v=-\frac{w^{2}}{12\unicode[STIX]{x1D707}\times \tilde{f}\left(Re_{D_{eq}},{\displaystyle \frac{k_{R}}{w}}\right)}\frac{\unicode[STIX]{x2202}p_{f}}{\unicode[STIX]{x2202}r}, & & \displaystyle\end{eqnarray}$$

where $\tilde{f}(Re_{D_{eq}},k_{R}/w)=f((4/3)Re,k_{R}/w)/f_{laminar}$ .

2.2.2 Limiting turbulent regimes

The friction factor $f$ depends on both the local value of the Reynolds number and relative roughness of the fracture, as depicted in figure 2. Although empirical models exist to somehow capture the complete evolution of $f$ , it is interesting to focus on limiting regimes for friction: laminar, turbulent rough, turbulent smooth etc. Explicit expressions for the friction factor exist in these different limiting turbulent regimes, and can be subsequently used to derive limiting solutions for hydraulic fracture growth.

Turbulent rough regime

At very large Reynolds numbers ( $Re\gg 10^{5}$ ), friction for plain water is solely governed by the relative roughness $k_{R}/w$ , where $k_{R}$ is the roughness scale. In that rough turbulent limiting regime, the friction factor evolves according to the Gauckler–Manning–Strickler scaling (Tsai & Rice Reference Tsai and Rice2010; Zia & Lecampion Reference Zia and Lecampion2017)

(2.12a-c ) $$\begin{eqnarray}\displaystyle f=f_{R}\times \left(\frac{k_{R}}{w}\right)^{\unicode[STIX]{x1D6FC}_{R}},\quad \unicode[STIX]{x1D6FC}_{R}=1/3,\quad f_{R}=2^{1/3}\times 0.143/4, & & \displaystyle\end{eqnarray}$$

where the factor $f_{R}$ allows for a perfect match with the experimental data (Nikuradse Reference Nikuradse1950).

Turbulent smooth regime

The turbulent regime for plain water in smooth pipes (where roughness is negligible) is well captured by the Blasius (Reference Blasius1913) scaling (see figure 2)

(2.13a-c ) $$\begin{eqnarray}\displaystyle f=f_{B}Re_{D_{eq}}^{-n_{B}},\quad n_{B}=1/4,\quad f_{B}=0.316/4, & & \displaystyle\end{eqnarray}$$

which is an empirical formula with a response similar to the Prandtl–Kármán theoretical predictions for smooth pipes.

Turbulent maximum drag reduction asymptote

As already mentioned, the addition of a friction reducer drastically changes the turbulent behaviour (Virk Reference Virk1975). When friction reducers are added in a sufficient concentration, the maximum drag reduction asymptote is reached. In that so-called MDR limit, the turbulent flow – even in a rough pipe, as can be seen from figure 2 – can be well captured by the following expression for Reynolds number in the range $[10^{3}{-}3\times 10^{4}]$

(2.14a-c ) $$\begin{eqnarray}\displaystyle f=f_{mdr}Re_{D_{eq}}^{-n_{mdr}},\quad n_{mdr}\approx 7/10,\quad f_{mdr}\approx 1.78. & & \displaystyle\end{eqnarray}$$

We will use this expression for simplicity in the following. It is interesting to point out that the last two regimes have a similar power-law relation between the friction factor and the Reynolds number, i.e. $f=f_{o}Re^{-n}$ .

2.2.3 Generalised Poiseuille law for limiting regimes

Combining the definition for the equivalent Reynolds number (2.10) and the relation between the width-averaged velocity and pressure gradient (2.9), it is possible to rewrite the general equation (2.9) for any limiting flow regime previously discussed as

(2.15) $$\begin{eqnarray}\displaystyle v\times |v|^{m-1}=-\frac{w^{n+1}}{c}\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}r}, & & \displaystyle\end{eqnarray}$$

where both the exponents $m$ and $n$ as well as the constant $c$ depend on the limiting flow regimes. For the laminar case, we simply have $c=12\unicode[STIX]{x1D707}$ , $m=1$ and $n=1$ . In the limiting turbulent flow regimes, these parameters take the following values:

1. (i) Turbulent rough – Gauckler–Manning–Strickler scaling $f=f_{R}\times (k_{R}/w)^{1/3}$ (2.12):

   (2.16a-c ) $$\begin{eqnarray}\displaystyle c=f_{R}k_{R}^{1/3}\unicode[STIX]{x1D70C},\quad n=1/3,\quad m=2. & & \displaystyle\end{eqnarray}$$

2. (ii) Turbulent regimes where $f=f_{o}Re^{-n}$ :

   (2.17a-c ) $$\begin{eqnarray}\displaystyle c=\underbrace{(4/3)^{-n}f_{o}}_{f_{o}^{\prime }}\unicode[STIX]{x1D707}^{n}\unicode[STIX]{x1D70C}^{1-n},\quad m=2-n,\quad f_{o}^{\prime }=(4/3)^{-n}f_{o}. & & \displaystyle\end{eqnarray}$$

The case of a smooth turbulent regime (Blasius scaling) as well as the case of the maximum drag reduction asymptote both fall under this form where $f=f_{o}Re^{-n}$ . Table 1 summarises the parameters of these three turbulent regimes in the generalised lubrication equation (2.15).

Table 1.

Generalised lubrication parameters of the different flow regimes. We use the notation $f_{o}^{\prime }=(4/3)^{-n}f_{o}$ for short; $f_{o}\approx 1.78$ for the MDR approximation, and $f_{o}\approx 0.316/4$ for turbulent smooth.

We therefore can use the general form of the averaged velocity/pressure gradient relation (2.15) to solve for the different limiting flow regimes in a similar way. Moreover, as first noted by Tsai & Rice (Reference Tsai and Rice2010) for the turbulent rough case, this general lubrication relation (2.15) is actually exactly similar in its mathematical form to the Poiseuille law for a power-law rheology (with $m=n$ ) for which $c$ and $n$ are the power-law consistency and fluid index respectively. One can therefore take advantage of the solution of hydraulic fracture problems obtained for a power-law rheology: notably, the zero toughness near-tip solution derived by Desroches et al. (Reference Desroches, Detournay, Lenoach, Papanastasiou, Pearson, Thiercelin and Cheng1994), as well as its extension to power-law fluid (Madyarova & Detournay Reference Madyarova and Detournay2004) for a penny-shaped hydraulic fracture propagating in a zero toughness medium, originally derived for a Newtonian fluid and laminar flow (Savitski & Detournay Reference Savitski and Detournay2002). We shall see that few adaptations are necessary to extend this zero toughness power-law laminar solution to the turbulent limiting regimes.

2.3 Inlet behaviour

The constant point-source injection (2.4) directly results in a $1/r$ singularity for the Reynolds number $Re=\unicode[STIX]{x1D70C}wv/\unicode[STIX]{x1D707}$ as it can be rewritten as $\lim _{r\rightarrow 0}2\unicode[STIX]{x03C0}rRe=\unicode[STIX]{x1D70C}Q_{o}/\unicode[STIX]{x1D707}$ . Obviously, in practice, injection always occurs from a finite wellbore (of radius  $a$ ) such that the inlet Reynolds number is limited to $Re_{max}=\unicode[STIX]{x1D70C}Q_{o}/(2\unicode[STIX]{x03C0}a\times \unicode[STIX]{x1D707})$ (see figure 1). Note that this also is the case for numerical evaluation, when using a finite grid to model the growth of a planar three-dimensional hydraulic fracture (Dontsov & Peirce Reference Dontsov and Peirce2017). Nevertheless, for simplicity we neglect the finiteness of the wellbore compared to the fracture radius in what follows. Combining the constant point-source injection boundary condition (2.4) with the generalised lubrication law (2.15), and recognising that the fracture width should be finite at the fracture inlet $r=0$ , we obtain the following inlet behaviour for the pressure gradient:

(2.18) $$\begin{eqnarray}\displaystyle \lim _{r\rightarrow 0}\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}r}\propto -\frac{1}{r^{m}}. & & \displaystyle\end{eqnarray}$$

The net pressure is thus singular at the inlet: $p\sim r^{1-m}$ . The singularity (associated with the point-source nature of the injection) is different for turbulent flow regimes compared to the laminar case ( $m=1$ ) where the pressure is log singular for a Newtonian fluid (Savitski & Detournay Reference Savitski and Detournay2002). Note that it is also different from a power-law fluid with fluid index $n$ under laminar conditions, where the pressure is not singular at the inlet as $p\sim r^{1-n}$ ( $0<n<1$ ).

Before further developing the solution of a finite radial hydraulic fracture, we now turn to examine the near-tip fracture behaviour.

3.1 Laminar Newtonian case

For a Newtonian fluid under laminar conditions ( $\tilde{f}=1$ ), the solution (see Garagash et al. Reference Garagash, Detournay and Adachi2011) evolves from the linear elastic fracture mechanics asymptote (subscript $k$ ) near the tip to a viscous asymptote in the far field (subscript $m$ ),

(3.4a-c ) $$\begin{eqnarray}\displaystyle w_{k}=\ell _{k}^{1/2}\hat{x}^{1/2},\quad p_{k}=0,\quad \ell _{k}=\left(\frac{K^{\prime }}{E^{\prime }}\right)^{2},\quad \hat{x}\ll 1, & & \displaystyle\end{eqnarray}$$

(3.5a-c ) $$\begin{eqnarray}\displaystyle w_{m}=\unicode[STIX]{x1D6FD}_{o}\ell _{m}^{1/3}\hat{x}^{2/3},\quad p_{m}=E^{\prime }\unicode[STIX]{x1D6FF}_{o}\frac{\ell _{m}^{1/3}}{\hat{x}^{1/3}},\quad \ell _{m}=V\frac{\unicode[STIX]{x1D707}^{\prime }}{E^{\prime }},\quad \hat{x}\gg 1, & & \displaystyle\end{eqnarray}$$

with $\unicode[STIX]{x1D6FD}_{o}=2^{1/3}\times 3^{5/6},\unicode[STIX]{x1D6FF}_{o}=(2\unicode[STIX]{x1D6FD}_{o}/12)\cot (\unicode[STIX]{x03C0}(2/3))$ . The viscous asymptote can be obtained by combining lubrication and elasticity assuming zero toughness – see e.g. Desroches et al. (Reference Desroches, Detournay, Lenoach, Papanastasiou, Pearson, Thiercelin and Cheng1994). A transition length scale $\ell _{mk}=K^{\prime 6}/(E^{\prime 4}V^{2}\unicode[STIX]{x1D707}^{\prime 2})$ can be defined as the location where these two width asymptotes intersect, i.e. $\ell _{k}^{1/2}\ell _{mk}^{1/2}=\ell _{m}^{1/3}\ell _{mk}^{2/3}$ (Garagash et al. Reference Garagash, Detournay and Adachi2011).

As discussed in Dontsov & Peirce (Reference Dontsov and Peirce2015), the complete solution for fracture width spanning the transition between the toughness and viscosity asymptotes can be well approximated (within a few per cent) as,

(3.6a-c ) $$\begin{eqnarray}\displaystyle \frac{w(\unicode[STIX]{x1D709}=\hat{x}/\ell _{mk})}{w_{mk}}=\unicode[STIX]{x1D709}^{1/2}\times (1+{\unicode[STIX]{x1D6FD}_{o}}^{3}\unicode[STIX]{x1D709}^{1/2})^{1/3},\quad w_{mk}=\frac{K^{\prime 4}}{E^{\prime 3}V\unicode[STIX]{x1D707}^{\prime }},\quad \ell _{mk}=\frac{K^{\prime 6}}{E^{\prime 4}V^{2}\unicode[STIX]{x1D707}^{\prime 2}},\quad \qquad & & \displaystyle\end{eqnarray}$$

where $w_{mk}$ and $\ell _{mk}$ are the characteristic width and length scales in the $mk$ edge scaling, spanning the transition between the near field $k$ and far-field $m$ asymptote (see also table 2 for a summary of all different transition scales).

Table 2.

Characteristic scales in the different possible edge scalings – note that $p_{\ast }=E^{\prime }w_{\ast }/\ell _{\ast }$ and $Re_{\ast }=\unicode[STIX]{x1D70C}w_{\ast }V/\unicode[STIX]{x1D707}$ .

3.2 Turbulent asymptote

The previous laminar solution is valid up to the abscissa $\hat{x}_{t}$ where the Reynolds number $Re=\unicode[STIX]{x1D70C}wV/\unicode[STIX]{x1D707}$ reaches the critical value for the onset of turbulent flow. Another asymptote – that we shall denote as the turbulent asymptote (subscript $t$ ) – can be obtained assuming that the turbulent flow is strictly in one of the limiting regimes previously discussed (Blasius, Gauckler–Manning–Strickler or maximum drag reduction case). The power-law-like lubrication relation (2.15) reads, in the moving tip coordinate system,

(3.7) $$\begin{eqnarray}\displaystyle V^{m}=\frac{w^{n+1}}{c}\frac{\text{d}p}{\text{d}\hat{x}}, & & \displaystyle\end{eqnarray}$$

such that the semi-infinite hydraulic fracture asymptote in the far field for large $\hat{x}$ /large $Re$ (where the effect of fracture toughness can be neglected) follows the zero toughness asymptote for a power-law fluid obtained by Desroches et al. (Reference Desroches, Detournay, Lenoach, Papanastasiou, Pearson, Thiercelin and Cheng1994), as pointed out by Tsai & Rice (Reference Tsai and Rice2010),

(3.8a-c ) $$\begin{eqnarray}\displaystyle w_{t}=\unicode[STIX]{x1D6FD}_{n}\ell _{t}^{n/(n+2)}\hat{x}^{2/(n+2)},\quad p_{t}=-\unicode[STIX]{x1D6FF}_{n}E^{\prime }\left(\frac{\ell _{t}}{\hat{x}}\right)^{n/(n+2)},\quad \ell _{t}=\left(\frac{cV^{m}}{E^{\prime }}\right)^{1/n}, & & \displaystyle\end{eqnarray}$$

with the prefactors (see Desroches et al. Reference Desroches, Detournay, Lenoach, Papanastasiou, Pearson, Thiercelin and Cheng1994)

(3.9a,b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FD}_{n}=\left(2\frac{(n+1)^{2}}{n}\tan \left[-\frac{2\unicode[STIX]{x03C0}}{n+2}\right]\right)^{1/(n+2)},\quad \unicode[STIX]{x1D6FF}_{n}=\frac{\unicode[STIX]{x1D6FD}_{n}(n+2)}{2}\cot \left(\frac{2\unicode[STIX]{x03C0}}{n+2}\right).\qquad & & \displaystyle\end{eqnarray}$$

This asymptote is notably valid for any type of limiting flow regime: turbulent rough (Gauckler–Manning–Strickler), turbulent smooth (Blasius) as well as the case with friction reducers at maximum drag reduction). It is worthwhile noting that the turbulent exponent $m$ does not enter into the functional form of the width and net pressure asymptote (but in the length scale $\ell _{t}$ ). Only the exponent $n$ appearing in the generalised lubrication relation (2.15) impacts this far-field asymptote.

3.3 Complete solution for $f\propto Re^{-\unicode[STIX]{x1D6FD}}$ in the turbulent regime for slickwater

It is possible to solve the complete semi-infinite fracture problem numerically for a given evolution of the friction factor $f$ – see Dontsov (Reference Dontsov2016a ) for the case of water without friction reducers. Here, focusing chiefly on the case of a friction reducer at maximum drag reduction, a simple piecewise power-law expression for the friction factor above and below a critical Reynolds number approximation appears appropriate (see figure 2)

(3.10) $$\begin{eqnarray}\displaystyle f=\left\{\begin{array}{@{}ll@{}}\begin{array}{@{}l@{}}12/Re,\\ f_{o}^{\prime }Re^{-\unicode[STIX]{x1D6FD}},\end{array}\quad & \begin{array}{@{}l@{}}Re<Re_{c},\\ Re>Re_{c},\end{array}\end{array}\right. & & \displaystyle\end{eqnarray}$$

with $Re=\unicode[STIX]{x1D70C}wV/\unicode[STIX]{x1D707}$ and $Re_{c}$ is chosen to ensure continuity of friction: $Re_{c}=(12/f_{o}^{\prime })^{1/(1-\unicode[STIX]{x1D6FD})}$ , i.e. $Re_{c}\approx 1132$ for the maximum drag reduction approximation. It is clear that a similar approximation for the case without a friction reducer results in a gross approximation of the transition region, which plays a significant role as discussed in Dontsov (Reference Dontsov2016a ). The method and qualitative evolution however remain the same. Using the scaling for the toughness to viscous asymptote transition in the laminar case (Garagash et al. Reference Garagash, Detournay and Adachi2011), we scale the coordinate, width and pressure as,

(3.11a-c ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D709}=\hat{x}/\ell _{mk},\quad \unicode[STIX]{x1D6FA}=w/w_{mk},\quad \unicode[STIX]{x1D6F1}=p/p_{mk}, & & \displaystyle\end{eqnarray}$$

using the characteristic scales corresponding to the laminar $mk$ transition previously mentioned (see also table 2)

(3.12a-c ) $$\begin{eqnarray}\displaystyle \ell _{mk}=\frac{K^{\prime 6}}{E^{\prime 4}V^{2}\unicode[STIX]{x1D707}^{\prime 2}},\quad w_{mk}=\frac{K^{\prime 4}}{E^{\prime 3}V\unicode[STIX]{x1D707}^{\prime }},\quad p_{mk}=\frac{E^{\prime 2}V\unicode[STIX]{x1D707}^{\prime }}{K^{\prime 2}}. & & \displaystyle\end{eqnarray}$$

The Reynolds number reduces to

(3.13a,b ) $$\begin{eqnarray}\displaystyle Re=\frac{\unicode[STIX]{x1D70C}wV}{\unicode[STIX]{x1D707}}=Re_{mk}\times \unicode[STIX]{x1D6FA},\quad Re_{mk}=12\frac{\unicode[STIX]{x1D70C}K^{\prime 4}}{E^{\prime 3}\unicode[STIX]{x1D707}^{\prime 2}}\left(=\frac{\unicode[STIX]{x1D70C}K^{\prime 4}}{12E^{\prime 3}\unicode[STIX]{x1D707}^{2}}\right). & & \displaystyle\end{eqnarray}$$

Here, $Re_{mk}$ is a characteristic Reynolds number which does not depend on the velocity of the fracture but solely on the fluid and solid material properties. The dependence on velocity lies in the transition length scale $\ell _{mk}$ . In such a scaling, elasticity, the fracture propagation condition and generalised lubrication flow read

(3.14) $$\begin{eqnarray}\displaystyle \hat{\unicode[STIX]{x1D6F1}}=\frac{1}{4\unicode[STIX]{x03C0}}\int _{0}^{\infty }\frac{\unicode[STIX]{x2202}\hat{\unicode[STIX]{x1D6FA}}(\unicode[STIX]{x1D702})}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}\frac{\text{d}\unicode[STIX]{x1D702}}{\unicode[STIX]{x1D709}-\unicode[STIX]{x1D702}}, & & \displaystyle\end{eqnarray}$$

(3.15a,b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FA}=\unicode[STIX]{x1D709}^{1/2},\quad \unicode[STIX]{x1D709}\ll 1, & & \displaystyle\end{eqnarray}$$

(3.16) $$\begin{eqnarray}\displaystyle 1=\frac{\hat{\unicode[STIX]{x1D6FA}}}{\tilde{f}(\hat{\unicode[STIX]{x1D6FA}})}\frac{\text{d}\hat{\unicode[STIX]{x1D6F1}}}{\text{d}\unicode[STIX]{x1D709}}, & & \displaystyle\end{eqnarray}$$

with

(3.17a,b ) $$\begin{eqnarray}\displaystyle \tilde{f}=1,\quad \hat{\unicode[STIX]{x1D6FA}}<Re_{c}/Re_{mk}, & & \displaystyle\end{eqnarray}$$

(3.18a,b ) $$\begin{eqnarray}\displaystyle \tilde{f}=\frac{f_{o}^{\prime }}{12}\times {Re_{mk}}^{1-n}\hat{\unicode[STIX]{x1D6FA}}^{1-n},\quad \hat{\unicode[STIX]{x1D6FA}}>Re_{c}/Re_{mk}. & & \displaystyle\end{eqnarray}$$

There is therefore a family of solutions parameterised by the dimensionless tip Reynolds number $Re_{mk}$ . We solve this system of nonlinear singular integro-differential equations numerically using a Gauss–Chebyshev quadrature mapping the interval $[0,\infty [$ to $[-1,1[$ and a simple finite difference approximation of the pressure gradient. The discretisation yields a nonlinear system of equations that can then be solved with standard root finding techniques. The method is strictly similar – pending the nonlinearity associated with $\tilde{f}(\hat{\unicode[STIX]{x1D6FA}})$ – to the one described in Moukhtari & Lecampion (Reference Moukhtari and Lecampion2018) and Viesca & Garagash (Reference Viesca and Garagash2018) and will not be repeated here.

Figure 3.

(a) Parametric diagram of the near-tip hydraulic fracture solution structure accounting for the transition to a turbulent flow regime ( $f\propto Re^{-n}$ ) with a few trajectories for different values of the tip Reynolds number $Re_{mk}$ ; (b) structure of the solution for small and large values of $Re_{mk}$ . For order $Re_{mk}\lesssim 100$ , the viscous laminar asymptote ( $m$ ) is clearly visible in between the far-field turbulent ( $t$ ) and near-tip linear elastic fracture mechanics asymptote ( $k$ ). For large $Re_{mk}$ , the laminar viscous asymptote disappears altogether: the solution transitions directly from $t$ to $k$ as one moves closer to the tip.

Figure 4 displays the results for different values of the tip Reynolds number $Re_{mk}$ . For small values of $Re_{mk}$ , the solution evolves from the toughness ( $k$ ) to the laminar viscous ( $m$ ) and ultimately the turbulent flow ( $t$ ) asymptotes as one moves away from the fracture tip. On the other hand, for large $Re_{mk}$ (above $10^{5}$ ), the viscous laminar region disappears and the solution transitions directly from the toughness ( $k$ ) to the turbulent ( $t$ ) flow asymptote, with the toughness region also shrinking. This is in line with the fact that for very large Reynolds numbers, the laminar region shrinks, as figure 3 shows. The transition length scale between the different limiting asymptotes can be obtained from arguments similar to the solely laminar case (i.e. equalising the width scales of the two different asymptotes to obtain the transition scale), thus defining different edge scalings (see table 2). From these transition length scales, and our numerical results, it is possible to estimate the distance from the tip $\hat{x}_{transition}$ where the flow transitions from laminar to turbulent: intuitively $\hat{x}_{transition}$ should scale with the length scale $\ell _{tm}$ governing the transition between the turbulent and laminar asymptotic regions. This is indeed the case, as can be seen in figure 5(a), where the location at which the local Reynolds number equals the critical value obtained from our numerical results, is compared with the transition length scale $\ell _{tm}$ as a function of the tip characteristic Reynolds number $Re_{mk}$ . One sees from figure 5(a) that the extent of the laminar layer is well approximated as

(3.19) $$\begin{eqnarray}\displaystyle \hat{x}_{transition}\approx 5000Re_{mk}^{-3/2}\ell _{mk}. & & \displaystyle\end{eqnarray}$$

From figure 5, we also observe that for values of $Re_{mk}$ greater than ${\sim}10^{3}$ , the turbulent to viscous laminar transition scale $\ell _{tm}$ becomes smaller than the turbulent to toughness transition scale $\ell _{tk}$ such that the laminar viscous asymptote ( $m$ ) disappears, as observed in the solution for width (figure 4). The domain of validity of the different tip asymptotes can also be well grasped from figure 5(b). It can be well observed that the laminar viscous asymptote disappears for $Re_{mk}>10^{3}$ while in the meantime the turbulent asymptotic region moves closer to the tip (and its location scales with $\ell _{tk}$ for large $Re_{mk}$ ). The fact that the extent of the laminar layer ( $\hat{x}_{transition}$ ) scales with $\ell _{tm}$ even for large Reynolds numbers is due to the fact that, even though the laminar viscous asymptote is not fully realised, laminar flow still influences the solution in the intermediate region. The transition from the turbulent to the toughness asymptote scales qualitatively with $\ell _{tk}$ for large $Re_{mk}$ . However $\ell _{tk}$ does not physically include the effect of laminar flow (and thus is not the proper scale for the laminar region).

Figure 4.

Near-tip solution for a slickwater hydraulic fracture accounting for turbulent flow at maximum drag reduction: numerical results (continuous lines) and $t$ , $m$ and $k$ asymptotes. Width as a function of tip distance (in the $mk$ scaling) for different values of the tip Reynolds number (inset: same results but scaled by the toughness near-tip asymptote $w_{k}$ ). The dots correspond to the boundary of the laminar flow region near the tip.

Figure 5.

(a) Evolution of the laminar region in the near-tip region as function of the tip Reynolds number. The laminar–turbulent transition distance $\hat{x}_{transition}$ is obtained from the numerical results as the distance (up to the discretisation) where $Re=Re_{c}$ (bullet points). The transition length scale $\ell _{tm}$ provides a qualitative estimate of the laminar region (continuous blue line), while a better approximation (with the same dependence on $Re_{mk}$ ) as the numerical results is obtained at $5000\times Re_{mk}^{-3/2}$ (dashed-dotted line). (b) Domains of validity of the different asymptotes ( $m,k,t$ ) as a function of normalised distance from the tip $\hat{x}/\ell _{mk}$ and $Re_{mk}$ . The boundaries of the domains correspond to the locations where the numerical solution is within a $5\,\%$ relative difference from the corresponding tip asymptote. The laminar–turbulent transition distance $\hat{x}_{transition}/\ell _{mk}$ is also plotted.

Such a simple estimation of the extent of the laminar layer (3.19) as well as the disappearance of the viscous asymptote for $Re_{mk}\gtrsim 10^{3}$ are extremely useful in practice to quickly determine the importance of taking into account the turbulent flow behaviour in a numerical simulation of hydraulic fracture growth. It is also worthwhile noting that the entire solution for the $tm$ and $tk$ edges of this near-tip solution can be very well approximated (up to 1.5 % of maximum relative error) by a simple first-order approximation similar to the one proposed for the laminar $mk$ edge by Dontsov & Peirce (Reference Dontsov and Peirce2015). These approximations are listed in appendix B for convenience.

This semi-infinite hydraulic fracture solution provides the near-tip width asymptote of the radial hydraulic fracture by replacing $V$ with the current radial fracture velocity $\text{d}R/\text{d}t$ (see discussions in Detournay & Peirce (Reference Detournay and Peirce2014), Dontsov & Kresse (Reference Dontsov and Kresse2018) and Lecampion, Bunger & Zhang (Reference Lecampion, Bunger and Zhang2018)).

4.1 Friction/viscous scalings

In the case where the energy dissipated in lubrication flow is larger than the energy spent in the creation of new surfaces, the dominance of the fluid friction dictates setting ${\mathcal{G}}_{f}$ to unity. We thus obtain the following estimate of the characteristic fracture length scale $L$ for the universal lubrication relation (2.15):

(4.13) $$\begin{eqnarray}\displaystyle L_{f}=\left(\frac{E^{\prime }Q_{o}^{n+2}t^{m+n+2}}{c}\right)^{1/(6+m+2n)}, & & \displaystyle\end{eqnarray}$$

where the subscript $f$ refers to such a fluid friction/viscous scaling. The dimensionless group ${\mathcal{G}}_{k}$ associated with fracture energy can be re-expressed in that scaling as a dimensionless toughness

(4.14) $$\begin{eqnarray}\displaystyle {\mathcal{K}}=K^{\prime }\times \left(\frac{t}{E^{\prime ^{(7+2m+4n)/(3m+n-2)}}Q_{o}^{5m/(3m+n-2)-1}c^{5/(3m+n-2)}}\right)^{(3m+n-2)/(2(6+m+2n))}, & & \displaystyle\end{eqnarray}$$

while the characteristic Reynolds number (4.9) in such a scaling is

(4.15) $$\begin{eqnarray}\displaystyle {\mathcal{R}}=\frac{\unicode[STIX]{x1D70C}Q_{o}}{\unicode[STIX]{x1D707}L_{f}}=\left(\frac{\unicode[STIX]{x1D70C}}{\unicode[STIX]{x1D707}}\frac{Q_{o}^{(4+m+n)/(2+m+n)}}{t}\left(\frac{c}{E^{\prime }}\right)^{1/(2+m+n)}\right)^{(2+m+n)/(6+m+2n)}. & & \displaystyle\end{eqnarray}$$

It is worthwhile noting that the dimensionless toughness increases with time. It can be expressed via the introduction of a characteristic time scale $t_{f\rightarrow k}$

(4.16a,b ) $$\begin{eqnarray}\displaystyle {\mathcal{K}}=\left(\frac{t}{t_{f\rightarrow k}}\right)^{3m+n-2/(2(6+m+2n))},\quad t_{f\rightarrow k}=\frac{E^{\prime (7+2m+4n)/(3m+n-2)}c^{5/(3m+n-2)}Q_{o}^{(2m+2-n)/(3m+n-2)}}{K^{\prime (2(6+m+2n))/(3m+n-2)}}.\hspace{-5.0pt} & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

Similarly, the characteristic Reynolds number decreases with time, indicating that turbulent flow dominates at early times, and becomes negligible when ${\mathcal{R}}<{\mathcal{R}}_{c}$ . This allows us to define a transition time scale for the transition towards the laminar regime $t_{t\rightarrow l}$ , such that the characteristic Reynolds number can be written as

(4.17a,b ) $$\begin{eqnarray}\displaystyle {\mathcal{R}}={\mathcal{R}}_{c}\left(\frac{t}{t_{t\rightarrow l}}\right)^{-((2+m+n)/(6+m+2n))},\quad t_{t\rightarrow l}=\frac{c^{1/(2+m+n)}Q_{o}^{1+(2/(2+m+n))}\unicode[STIX]{x1D70C}^{(6+m+2n)/(2+m+n)}}{E^{\prime 1/(2+m+n)}\unicode[STIX]{x1D707}^{2+(2-m/2+m+n)}{\mathcal{R}}_{c}^{2+(2-m/(2+m+n))}}. & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

Table 3.

Characteristic length scale $L_{f}$ and dimensionless toughness in the different fluid friction dominated scalings. The corresponding width scales as $W=\unicode[STIX]{x1D716}L_{f}=Q_{o}t/L_{f}^{2}$ , while the characteristic pressure is $P=E^{\prime }Q_{o}t/L_{f}^{3}$ . Note that the complete friction factor depends in the general case on both Reynolds number and relative roughness. The Reynolds number scales with the characteristic value ${\mathcal{R}}=\unicode[STIX]{x1D70C}Q_{o}/(\unicode[STIX]{x1D707}\,L_{f}(t))$ , while the relative roughness scales as $k_{R}/W$ . For the turbulent smooth and MDR regimes the corresponding expressions for $f_{o}$ are listed in (2.13) and (2.14) respectively.

Such a friction scaling is valid for the different flow regimes (different $m,n$ and  $c$ ) where the lubrication velocity–pressure gradient relation can be expressed in the general form (2.15). The exact expressions for the turbulent rough, turbulent smooth as well as the maximum drag reduction regimes can be found in table 3. We obviously recover the viscosity scaling in the laminar case. It is interesting to note that the power laws of time for the characteristic length are $1/2$ , $16/33$ and $40/87$ respectively for the rough, smooth and MDR regimes – all very close to $1/2$ . As a result, the fracture width, which scales as $W=Q_{o}t/L_{f}^{2}$ , will either be constant (rough regimes) or evolves slowly as $t^{1/33}$ and $t^{7/87}$ .

It is important to keep in mind that these scalings are based on the underlying assumption that the flow regime is the same over the whole fracture length. Changes in the flow regimes will occur within the fracture in relation to the local value of the Reynolds number which depends on the characteristic Reynolds number ${\mathcal{R}}$ which scales as $1/L_{f}$ , as well as the local dimensionless fracture width and velocity (see (4.8)). Therefore, a transition between flow regimes will necessarily occur within the growing fracture as time progresses and a laminar region, albeit small, always exist near the fracture tip. Nevertheless, at early times, one can build a limiting self-similar solution for the growth of a radial hydraulic fracture for each of the flow regimes: laminar, rough turbulent, smooth turbulent or at maximum drag reduction (in the case with friction reducing additives).

4.2 Toughness scaling

The other limiting scaling corresponds to the case where most of the energy is dissipated in the solid, creating new fracture surfaces. Setting ${\mathcal{G}}_{k}=1$ , one recovers the well-known toughness dominated scaling (Savitski & Detournay Reference Savitski and Detournay2002), where the characteristic length scale is independent of the fluid properties

(4.18) $$\begin{eqnarray}\displaystyle L_{k}=\left(\frac{E^{\prime }Q_{o}t}{K^{\prime }}\right)^{2/5} & & \displaystyle\end{eqnarray}$$

and the corresponding dimensionless ‘viscosity’ ${\mathcal{G}}_{f}$ becomes

(4.19) $$\begin{eqnarray}\displaystyle {\mathcal{M}}=c\times \left(\frac{Q_{o}^{(2+2m-n)/(3m+n-2)}E^{\prime (7+2m+4n)/(3m+n-2)}}{K^{\prime ((2(6+m+2n))/(2m+n-2))}t}\right)^{5/(3m+n-2)}. & & \displaystyle\end{eqnarray}$$

There is of course a correspondence with the friction/lubrication scalings

(4.20a,b ) $$\begin{eqnarray}\displaystyle {\mathcal{K}}={\mathcal{M}}^{-5/(2(6+m+2n))},\quad \frac{L_{k}}{L_{f}}={\mathcal{K}}^{-2/5}. & & \displaystyle\end{eqnarray}$$

4.3 Structure of the solution

In the case of solely laminar flow, the hydraulic fracture growth solution evolves from the viscosity dominated regime of propagation at small times to the toughness dominated regime at large times (Savitski & Detournay Reference Savitski and Detournay2002). The time scale of this transition $t_{m\rightarrow k}$ is the one appearing in the expression for the dimensionless toughness ${\mathcal{K}}$ (see table 3). Such a transition can be physically grasped as the energy spent in creating new surfaces increases as $2\unicode[STIX]{x03C0}R\,\text{d}R$ for such an axisymmetric fracture and thus necessarily dominates at large times. In the laminar case, it was shown that the viscosity dominated regime is valid up to ${\mathcal{K}}=1$ while the toughness dominated regime starts for dimensionless toughness larger than $3.5$ (Savitski & Detournay Reference Savitski and Detournay2002).

When accounting for turbulent regimes, the solution structure is similar at large times, as here again the dimensionless toughness ${\mathcal{K}}$ increases with time as $\propto t^{\unicode[STIX]{x1D6FC}}$ (with $\unicode[STIX]{x1D6FC}>0$ ), as can be seen in table 3 for all the different limiting flow regimes. However, other time scales appear associated with the transition between the different flow regimes as functions of the Reynolds number and characteristic roughness. We first discuss the case of plain water without friction reducing additives, where both rough and smooth turbulent regimes can occur.

4.3.1 Water (without friction reducers)

As discussed previously, the characteristic Reynolds number evolves as $1/L$ (4.9), such that turbulent regimes will dominate at the early times of fracture growth. Moreover, at early times, fracture roughness is bound to dominate as the fracture width will be small and comparable to the roughness scale. The following picture therefore emerges for the case of plain water without friction reducers. At early times, the propagation will be dominated by the rough turbulent regime. As time progresses, one may transition to the turbulent smooth regime as the fracture widens compared to the roughness scales, then to the laminar viscous dominated regime and ultimately to the toughness dominated regime. However, it may well be possible that the HF growth directly transitions from the turbulent rough to the toughness dominated regime. In order to grasp the transition between the different flow regimes, one can easily determine the time at which the characteristic fracture scale $L_{f}$ is equal for both limiting regimes. In doing so, one can define three time scales corresponding respectively to (i) the transition between the rough and smooth turbulent regimes ( $t_{r\rightarrow s}$ ), (ii) the turbulent rough and laminar regimes ( $t_{r\rightarrow m}$ ) and (iii) the turbulent smooth and laminar regimes ( $t_{s\rightarrow m}$ ). All these time scales are listed in table 4. The transition time scales from the turbulent to the laminar regime could also have been obtained by arguing that the laminar regime is reached when the characteristic Reynolds number ${\mathcal{R}}$ falls below the critical laminar–turbulent transition value ${\mathcal{R}}_{c}$ (similar results are obtained up to a constant factor).

Table 4.

Transition time scales between regimes of propagation governed by the different limiting flow regimes.

As a result, four different situations may be encountered during the growth of a radial fracture, as depicted in figure 6, depending on the relative value of the transition time scales between flow regimes as well as the transition time scales from a given flow regime to the toughness dominated regime (table 3).

Figure 6.

Propagation diagram for a radial hydraulic fracture: water without friction reducers (a) and slickwater with friction reducer at maximum drag reduction (b).

For the sake of illustration, let us consider the following values of the parameters (corresponding to water injected at a very high rate in an ‘average’ sedimentary rock):

(4.21a-c ) $$\begin{eqnarray}\displaystyle E^{\prime }=20~\text{GPa},\quad Q_{o}=0.079~\text{m}^{3}~\text{s}^{-1},\quad \unicode[STIX]{x1D70C}=1000~\text{kg}~\text{m}^{-3}, & & \displaystyle\end{eqnarray}$$

(4.21d-f ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D707}=1~\text{cP},\quad K_{Ic}=0.5~\text{MPa}\,\sqrt{\text{m}},\quad k_{R}=500~\unicode[STIX]{x03BC}\text{m}, & & \displaystyle\end{eqnarray}$$

for which we obtain the following transition time scales between the different limiting regimes (all times in seconds):

(4.22) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}l@{}}t_{r\rightarrow s}=97\,681,\quad t_{r\rightarrow m}=844,\quad t_{s\rightarrow m}=142,\\ t_{r\rightarrow k}=5049,\quad t_{s\rightarrow k}=2975,\quad t_{m\rightarrow k}=47\,245.\end{array}\right\} & & \displaystyle\end{eqnarray}$$

For that particular set of parameters, the propagation will evolve from the roughness turbulent dominated regime, toward the laminar viscous dominated regime to ultimately the toughness dominate regime at large times: along trajectory (2) in figure 6. The smooth turbulent regime will not appear as $t_{r\rightarrow s}\gg t_{r\rightarrow m}$ . Moreover the transition to the laminar viscous dominated regime occurs faster than the transition from the rough to the toughness dominated regime $t_{r\rightarrow m}\ll t_{r\rightarrow k}$ .

Note that other types of trajectories are possible as functions of the problem parameters. For a larger toughness, a transition directly from the rough turbulent to the toughness regime can occur (trajectory (4) in figure 6), or for a smaller value of the roughness scale $k_{R}$ , the smooth turbulent regimes may appear (trajectories (1) and (3) in figure 6). The diagram and limiting trajectories depicted in figure 6 in combination with the different transition time scales allows us to directly quantify the regimes encountered during the growth of a radial hydraulic fracture without doing any numerical simulations.

4.3.2 Slickwater (water with friction reducers)

In the case where friction reducing agents are added in sufficient concentration to water, the turbulent flow is governed by the drag reduction asymptote and roughness appears to have little effect here for Reynolds numbers lower than $3\times 10^{4}$ (see figure 2). From the expected ranges of Reynolds number in applications (always lower than $10^{5}$ ), the situation is therefore simpler than for the previously discussed water case. Beside the transition towards the toughness dominated regime at large times, there exists here only one additional time scale governing the transition from the turbulent/MDR regime toward the laminar viscous dominated one. Similarly, the expression for such a time scale $t_{mdr\rightarrow m}$ can be easily obtained by solving for the time at which the characteristic fracture scales in the MDR and laminar regimes are equal – see table 4.

The diagram picturing the possible trajectory from one regime to the next for a radial HF growing with friction reducing additives is depicted in figure 6(b). For the same set of parameters used previously, we obtain in the slickwater case the following transition time scales (time in seconds)

(4.23) $$\begin{eqnarray}\displaystyle t_{mdr\rightarrow m}=83,\quad t_{mdr\rightarrow k}=9288,\quad t_{m\rightarrow k}=47\,245. & & \displaystyle\end{eqnarray}$$

For that particular case, the HF quickly transitions to the viscous laminar dominated regime ( $t_{mdr\rightarrow m}\ll t_{mdr\rightarrow k}$ ) and then moves toward the toughness dominated regime at large times. For a larger value of fracture toughness (2.5 instead of 0.5), we obtain the opposite case $t_{mdr\rightarrow m}(=30~\text{s})>t_{mdr\rightarrow k}(=30~\text{s})$ such that the HF will move directly from the MDR dominated to the toughness dominated regime.

5.1 Method of solution

The method of solution is similar to the one described in Savitski & Detournay (Reference Savitski and Detournay2002). It consists in seeking a solution ${\mathcal{F}}(\bar{\unicode[STIX]{x1D6FA}},\unicode[STIX]{x1D6F1},\unicode[STIX]{x1D6FE})$ as ${\mathcal{F}}={\mathcal{F}}^{\ast }+B{\mathcal{F}}^{\ast \ast }$ , where ${\mathcal{F}}^{\ast \ast }$ denotes a particular solution satisfying the inlet behaviour associated with the boundary conditions (see § 2.3), while ${\mathcal{F}}^{\ast }$ is a general solution. We write the net pressure and reduced opening profiles as the following series:

(5.9) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6F1}(\unicode[STIX]{x1D709})=\mathop{\sum }_{i=1}^{\infty }A_{i}\unicode[STIX]{x1D6F1}_{i}^{\ast }(\unicode[STIX]{x1D709})+B\unicode[STIX]{x1D6F1}^{\ast \ast }(\unicode[STIX]{x1D709}), & \displaystyle\end{eqnarray}$$

(5.10) $$\begin{eqnarray}\displaystyle & \displaystyle \bar{\unicode[STIX]{x1D6FA}}(\unicode[STIX]{x1D709})=\mathop{\sum }_{i=1}^{\infty }C_{i}\bar{\unicode[STIX]{x1D6FA}}_{i}^{\ast }(\unicode[STIX]{x1D709})+B\bar{\unicode[STIX]{x1D6FA}}^{\ast \ast }(\unicode[STIX]{x1D709}), & \displaystyle\end{eqnarray}$$

where $A_{i},C_{i}$ and $B$ are unknowns coefficients. The general basis functions $\unicode[STIX]{x1D6F1}_{i}^{\ast }$ and $\bar{\unicode[STIX]{x1D6FA}}_{i}^{\ast }$ are taken to be a pairwise solution of the elasticity equation (5.2) and to also satisfy the zero toughness condition (5.4) as well as the asymptotic behaviour at the fracture tip (3.8) – see Savitski & Detournay (Reference Savitski and Detournay2002) for more details. The particular function $\unicode[STIX]{x1D6F1}^{\ast \ast }$ is taken to satisfy the inlet pressure singularity as well as the zero toughness condition, and the function $\bar{\unicode[STIX]{x1D6FA}}^{\ast \ast }$ is the corresponding solution of the elastic equation. After the proper choice of general and particular solutions, we will discretise the continuity equation (5.7) at a number of collocation points and solve for the unknown coefficients of the truncated series (5.9)–(5.10).

5.1.1 Particular inlet solution

As discussed in sub-section (2.3), under the assumption of a finite width at the fracture inlet, from the inlet flux boundary conditions, the gradient of net pressure at the inlet behaves as

(5.11) $$\begin{eqnarray}\displaystyle \frac{\text{d}\unicode[STIX]{x1D6F1}}{\text{d}\unicode[STIX]{x1D709}}\propto -\frac{1}{\unicode[STIX]{x1D709}^{m}} & & \displaystyle\end{eqnarray}$$

such that after integration, for $m\neq 1$ , we obtain the following particular solution for the net pressure which has the proper behaviour for $\unicode[STIX]{x1D709}\rightarrow 0$ :

(5.12) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6F1}^{\ast \ast }=\frac{1}{m-1}\frac{1}{\unicode[STIX]{x1D709}^{m-1}}-b. & & \displaystyle\end{eqnarray}$$

The constant $b$ is determined such that this particular solution satisfies the condition of zero toughness propagation (5.4),

(5.13) $$\begin{eqnarray}\displaystyle b=-\frac{\sqrt{\unicode[STIX]{x03C0}}}{4}\frac{\unicode[STIX]{x1D6E4}\left({\displaystyle \frac{1-m}{2}}\right)}{\unicode[STIX]{x1D6E4}(2-m/2)}. & & \displaystyle\end{eqnarray}$$

The fracture width $\bar{\unicode[STIX]{x1D6FA}}^{\ast \ast }$ corresponding to the previous particular net pressure profile can be obtained from the elasticity equation (5.2). After integration, one obtains

(5.14) $$\begin{eqnarray}\displaystyle \bar{\unicode[STIX]{x1D6FA}}^{\ast \ast }(\unicode[STIX]{x1D709}) & = & \displaystyle \frac{2}{\sqrt{\unicode[STIX]{x03C0}}}\frac{\unicode[STIX]{x1D6E4}\left({\displaystyle \frac{1-m}{2}}\right)}{\unicode[STIX]{x1D6E4}\left({\displaystyle \frac{4-m}{2}}\right)}\left(\sqrt{1-\unicode[STIX]{x1D709}^{2}}-\frac{\sqrt{\unicode[STIX]{x03C0}}}{2}\frac{\unicode[STIX]{x1D6E4}\left({\displaystyle \frac{m-2}{2}}\right)}{\unicode[STIX]{x1D6E4}\left({\displaystyle \frac{m-1}{2}}\right)}\unicode[STIX]{x1D709}^{2-m}\right.\nonumber\\ \displaystyle & & \displaystyle -\,\left.\frac{\text{}_{2}F_{1}\left({\displaystyle \frac{1}{2}},{\displaystyle \frac{m-2}{2}},{\displaystyle \frac{1}{2}},\unicode[STIX]{x1D709}^{2}\right)}{2-m}\right)\quad \text{for }m\neq 2,\end{eqnarray}$$

(5.15) $$\begin{eqnarray}\displaystyle \bar{\unicode[STIX]{x1D6FA}}^{\ast \ast }(\unicode[STIX]{x1D709}) & = & \displaystyle -4(\sqrt{1-\unicode[STIX]{x1D709}^{2}}+\log \unicode[STIX]{x1D709}-\log (1+\sqrt{1-\unicode[STIX]{x1D709}^{2}}))\quad \text{for }m=2,\end{eqnarray}$$

where $\unicode[STIX]{x1D6E4}$ denotes the Euler gamma and $_{2}F_{1}$ the hyper-geometric function. Such particular solutions ${\mathcal{F}}^{\ast \ast }(\unicode[STIX]{x1D6F1}^{\ast \ast },\bar{\unicode[STIX]{x1D6FA}}^{\ast \ast })$ satisfy the inlet behaviour of the net pressure, the elasticity equation as well as the zero toughness condition.

The inlet logarithmic singularity for the fracture width obtained for the rough turbulent regimes ( $m=2$ ) violates the original assumption of a finite inlet width. Zolfaghari & Bunger (Reference Zolfaghari and Bunger2018a ) have discussed the use of a more appropriate particular inlet solution for that particular case. In what follows, we nevertheless use the previous functional form in the series for simplicity. This has little impact on the obtained results, as we shall see.

5.1.2 General solution

Following Savitski & Detournay (Reference Savitski and Detournay2002), Madyarova & Detournay (Reference Madyarova and Detournay2004), we choose the following general basis for the net pressure and corresponding opening general solution

(5.16) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6F1}_{i}^{\ast }(\unicode[STIX]{x1D709})=\underbrace{-\frac{(1-\unicode[STIX]{x1D709})^{-n/(2+n)}}{\sqrt{h_{i-1}\left({\displaystyle \frac{4}{2+n}},2\right)}}J_{i-1}\left(\frac{4}{2+n},2,\unicode[STIX]{x1D709}\right)}_{\bar{\unicode[STIX]{x1D6F1}}^{\ast }}+\,\unicode[STIX]{x1D714}_{i}, & \displaystyle\end{eqnarray}$$

(5.17) $$\begin{eqnarray}\displaystyle & \displaystyle \bar{\unicode[STIX]{x1D6FA}}_{i}^{\ast }(\unicode[STIX]{x1D709})=\frac{(1-\unicode[STIX]{x1D709})^{2/(2+n)}}{\sqrt{h_{i-1}\left({\displaystyle \frac{2(4+n)}{2+n}},2\right)}}J_{i-1}\left(\frac{2(4+n)}{2+n},2,\unicode[STIX]{x1D709}\right), & \displaystyle\end{eqnarray}$$

where $J$ is the Jacobi polynomial of order $i$ and $h_{i}$ is the norm of the corresponding Jacobi polynomial. Both functions $\unicode[STIX]{x1D6F1}^{\ast }{-}\bar{\unicode[STIX]{x1D6FA}}^{\ast }$ provide an orthogonal basis. They satisfy the elasticity equation (5.2) as well as the zero toughness tip asymptotes (3.8), i.e.

(5.18) $$\begin{eqnarray}\displaystyle & \displaystyle \bar{\unicode[STIX]{x1D6FA}}_{i}^{\ast }\propto (1-\unicode[STIX]{x1D709})^{2/(2+n)},\quad 1-\unicode[STIX]{x1D709}\ll 1, & \displaystyle\end{eqnarray}$$

(5.19) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6F1}_{i}^{\ast }\propto (1-\unicode[STIX]{x1D709})^{-n/(2+n)},\quad 1-\unicode[STIX]{x1D709}\ll 1. & \displaystyle\end{eqnarray}$$

The constants $\unicode[STIX]{x1D714}_{i}$ ensure that the zero toughness condition (5.4) is satisfied and are thus simply obtained as,

(5.20) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D714}_{i}=-\int _{0}^{1}\frac{\bar{\unicode[STIX]{x1D6F1}}_{i}^{\ast }(\unicode[STIX]{x1D709})}{\sqrt{1-\unicode[STIX]{x1D709}^{2}}}\unicode[STIX]{x1D709}\,\text{d}\unicode[STIX]{x1D709}. & & \displaystyle\end{eqnarray}$$

Truncating the series, using $n_{\unicode[STIX]{x1D6F1}}$ terms for the net pressure and $n_{\unicode[STIX]{x1D6FA}}=n_{\unicode[STIX]{x1D6F1}}+1$ for the opening, one obtains the following relations between the $C$ and $A$ coefficients via the elastic relation:

(5.21) $$\begin{eqnarray}\displaystyle & \displaystyle C_{i}^{(n_{\unicode[STIX]{x1D6FA}},n_{\unicode[STIX]{x1D6F1}})}=\mathop{\sum }_{j=1}^{n_{\unicode[STIX]{x1D6F1}}}L_{ij}A_{j}^{(n_{\unicode[STIX]{x1D6FA}},n_{\unicode[STIX]{x1D6F1}})}, & \displaystyle\end{eqnarray}$$

(5.22) $$\begin{eqnarray}\displaystyle & \displaystyle L_{ij}=\frac{8}{\unicode[STIX]{x03C0}}\int _{0}^{1}\bar{\unicode[STIX]{x1D6FA}}_{i}(\unicode[STIX]{x1D702})\left[\int _{0}^{1}G(\unicode[STIX]{x1D702},\unicode[STIX]{x1D709})\unicode[STIX]{x1D6F1}_{j}^{\ast }(\unicode[STIX]{x1D709})\unicode[STIX]{x1D709}\,\text{d}\unicode[STIX]{x1D709}\right]\unicode[STIX]{x1D702}\,\text{d}\unicode[STIX]{x1D702}. & \displaystyle\end{eqnarray}$$

The coefficients $L_{ij}$ can be pre-computed (see Savitski & Detournay (Reference Savitski and Detournay2002) for more details).

5.1.3 Discretised lubrication equation

With the previously discussed general and particular solutions in hand, we collocate the lubrication equation (5.7) at $N$ distinct spatial points $\unicode[STIX]{x1D709}_{j}$ in the interval $[0,1]$ and construct the following error cost function to be minimised in terms of the coefficients $\boldsymbol{A}=(A_{1},\ldots ,A_{i},\ldots )$ and $B$ :

(5.23) $$\begin{eqnarray}\displaystyle \text{Err}(\boldsymbol{A},B)=\mathop{\sum }_{j=1}^{N}\left\Vert \frac{\text{rhs}(\boldsymbol{A},B,\unicode[STIX]{x1D709}_{j})}{\text{lhs}(\boldsymbol{A},B,\unicode[STIX]{x1D709}_{j})}-1\right\Vert ^{2}, & & \displaystyle\end{eqnarray}$$

where the ‘rhs’ and ‘lhs’ stand for the right-hand side and left-hand side of the lubrication equation (5.7). The details of their expressions are listed in appendix C. For a given order of truncation $n_{\unicode[STIX]{x1D6F1}}$ , we use the built in NMininize function of Mathematica to minimise this cost function and obtain the corresponding series coefficients.

5.3 Approximate solutions

Simple approximate solutions for these limiting zero toughness propagation regimes can also be devised following an approximate approach similar to the one described by Dontsov (Reference Dontsov2016b ) for the laminar case. The underlying idea is to recognise that the near-tip asymptote governs most of the fracture width, such that the fracture width profile of the finite radial hydraulic fracture is approximated as

(5.24) $$\begin{eqnarray}\displaystyle w(r)=\left(\frac{R+r}{2R}\right)^{\unicode[STIX]{x1D706}}w_{a}(R-r), & & \displaystyle\end{eqnarray}$$

where $w_{a}$ corresponds to the near-tip solution which depends on material parameters and the fracture velocity; $\unicode[STIX]{x1D706}$ is a factor to stitch the asymptote solution to a finite fracture profile – we will use $\unicode[STIX]{x1D706}=1/2$ in the following for simplicity. Assuming a limiting turbulent propagation regime (turbulent rough, smooth or at maximum drag reduction), the growth of the fracture is self-similar, e.g. $R=\unicode[STIX]{x1D6FE}L_{f}(t)$ , while the near-tip asymptote $w_{a}$ is given by the turbulent tip asymptote (3.8). The global volume balance then provides a simple equation to solve for the dimensionless fracture length $\unicode[STIX]{x1D6FE}$ . In the turbulent lubrication scaling, one obtains (using $\text{d}R/\text{d}t=\unicode[STIX]{x1D6FE}\text{d}L_{f}/\text{d}t$ as the fracture velocity in the expression of the tip asymptote (3.8)),

(5.25) $$\begin{eqnarray}\displaystyle 2\unicode[STIX]{x03C0}\left(\frac{2+m+n}{6+2m+n}\right)^{m/(2+n)}\unicode[STIX]{x1D6FE}^{(6+m+2n)/(2+n)}\int _{0}^{1}\left(\frac{1+\unicode[STIX]{x1D709}}{2}\right)^{\unicode[STIX]{x1D706}}\times \unicode[STIX]{x1D6FD}_{n}\times (1-\unicode[STIX]{x1D709})^{2/(n+2)}\unicode[STIX]{x1D709}\,\text{d}\unicode[STIX]{x1D709}=1, & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

which can be easily solved for $\unicode[STIX]{x1D6FE}$ for a given flow regime ( $m$ and $n$ ), and a choice of $\unicode[STIX]{x1D706}$ .

Following Dontsov (Reference Dontsov2016b ), using $\unicode[STIX]{x1D706}=1/2$ , we obtain the following approximate solutions for the different turbulent flow regimes – expressed here in dimensionless form (see table 3 for the corresponding scaling):

1. (i) turbulent rough

   (5.26a,b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FE}\approx 0.8473,\quad \bar{\unicode[STIX]{x1D6FA}}(\unicode[STIX]{x1D709})=\unicode[STIX]{x1D6FA}(\unicode[STIX]{x1D709})/\unicode[STIX]{x1D6FE}\approx 1.1466(1-\unicode[STIX]{x1D709})^{6/7}(1+\unicode[STIX]{x1D709})^{1/2}, & & \displaystyle\end{eqnarray}$$

2. (ii) turbulent smooth

   (5.27a,b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FE}\approx 0.8417,\quad \bar{\unicode[STIX]{x1D6FA}}(\unicode[STIX]{x1D709})=\unicode[STIX]{x1D6FA}(\unicode[STIX]{x1D709})/\unicode[STIX]{x1D6FE}\approx 1.1868(1-\unicode[STIX]{x1D709})^{8/9}(1+\unicode[STIX]{x1D709})^{1/2}, & & \displaystyle\end{eqnarray}$$

3. (iii) turbulent with friction reducers at maximum drag reduction

   (5.28a,b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FE}\approx 0.7534,\quad \bar{\unicode[STIX]{x1D6FA}}(\unicode[STIX]{x1D709})=\unicode[STIX]{x1D6FA}(\unicode[STIX]{x1D709})/\unicode[STIX]{x1D6FE}\approx 1.4369(1-\unicode[STIX]{x1D709})^{0.7407}(1+\unicode[STIX]{x1D709})^{1/2}. & & \displaystyle\end{eqnarray}$$

One can obtain the corresponding net pressure profiles via the elastic equation. It is interesting to note that the dimensionless lengths obtained are in very good agreement with the more rigorous solution previously obtained with relative differences of $1.15$ , $1.8$ and 0.6 % for the rough, smooth and MDR cases respectively. However, the width profiles are not as precise as these approximations do not capture the details of the fracture inlet behaviour, as can be observed in figure 7. These approximate solutions nevertheless provide a simple estimate of fracture length and opening that can be easily used in practice.