2 Mathematical model

We consider the problem of a semi-infinite fluid-driven fracture propagating with a constant velocity $V$ in a permeable elastic rock under condition of plane strain. The analysis is carried out on the basis of the following assumptions: (i) fracture propagates according to the linear elastic fracture mechanics (LEFM, see Rice Reference Rice and Liebowitz1968); (ii) the flow of incompressible power-law fluid in the crack can be modelled by the lubrication theory (Batchelor Reference Batchelor1967; Ben-Naceur Reference Ben-Naceur, Economies and Nolte1989); (iii) the fluid loss velocity is given by Carter’s leak-off model (Carter Reference Carter, Howard and Fast1957); (iv) the fluid front coincides with the fracture tip, i.e. there is no fluid lag (Detournay & Garagash Reference Detournay and Garagash2003).

In addition to the fracture tip velocity, $V$ , the material parameters that characterize the problem are

(2.1a-d ) $$\begin{eqnarray}E^{\prime }=\frac{E}{1-\unicode[STIX]{x1D708}^{2}},~~K^{\prime }=4\sqrt{\frac{2}{\unicode[STIX]{x03C0}}}K_{Ic},~~C^{\prime }=2C_{l},~~M^{\prime }={\displaystyle \frac{2^{n+1}(2n+1)^{n}}{n^{n}}}K,\end{eqnarray}$$

where $E$ denotes the Young’s modulus, $\unicode[STIX]{x1D708}$ is the Poisson’s ratio, $K_{Ic}$ is the rock fracture toughness, $C_{l}$ is Carter’s leak-off coefficient, $K$ is the consistency index and $n$ is the power-law exponent. Note that the power-law rheology corresponds to the situation when the shear stress $\unicode[STIX]{x1D70F}$ and the shear rate $\dot{\unicode[STIX]{x1D6FE}}$ for the fluid are related though the relation $\unicode[STIX]{x1D70F}=K\dot{\unicode[STIX]{x1D6FE}}^{n}$

To address the problem, we introduce the coordinate $x$ , such that the origin is moving with the constant velocity $V$ and is always located at the fracture tip. In this case, the positive $x$ values correspond to the distance from a point inside the fracture to the tip, as shown in figure 1(a). Within the given moving coordinate system, the fracture width is denoted by $w(x)$ , the net pressure (which is equal to the difference between the fluid pressure $p_{f}$ and the far-field stress $\unicode[STIX]{x1D70E}_{0}$ ) by $p(x)$ , the fluid flux by $q(x)$ and fluid loss velocity by $g(x)$ . The system of governing equations that describes the crack propagation is briefly outlined below.

Figure 1. (a) Schematics of a semi-infinite hydraulic fracture model. (b) The triangular parametric space, see Garagash et al. (Reference Garagash, Detournay and Adachi2011) for a similar representation for the case of Newtonian fluids. The solid lines inside the triangle schematically indicate trajectories of the solution in the parametric space for different values of leak-off.

Elasticity equation. The elastic relation between $p(x)$ and $w(x)$ is expressed by the singular integral equation (Hills et al. Reference Hills, Kelly, Dai and Korsunsky1996)

(2.2) $$\begin{eqnarray}p(x)=\frac{E^{\prime }}{4\unicode[STIX]{x03C0}}\int _{0}^{\infty }\frac{\unicode[STIX]{x2202}w}{\unicode[STIX]{x2202}s}\,\frac{\text{d}s}{x-s}.\end{eqnarray}$$

Poiseuille law. According to the lubrication theory, the equation that governs the flow of a power-law fluid within the fracture is given by (Ben-Naceur Reference Ben-Naceur, Economies and Nolte1989)

(2.3) $$\begin{eqnarray}q|q|^{n-1}=\frac{w^{2n+1}}{M^{\prime }}\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}x}.\end{eqnarray}$$

Leak-off law. The fluid loss velocity $g(x)$ for a steadily propagating crack, which accounts for both fracture faces, is given by (Carter Reference Carter, Howard and Fast1957)

(2.4) $$\begin{eqnarray}g(x)=C^{\prime }\sqrt{\frac{V}{x}.}\end{eqnarray}$$

Fluid mass balance. By assuming an incompressible fracturing fluid, the local fluid mass balance, that is integrated over the crack width, can be expressed by the continuity equation

(2.5) $$\begin{eqnarray}V\frac{\unicode[STIX]{x2202}w}{\unicode[STIX]{x2202}x}-\frac{\unicode[STIX]{x2202}q}{\unicode[STIX]{x2202}x}+g=0.\end{eqnarray}$$

Integration of this expression with respect to $x$ and the use of (2.4) yields the global volume balance

(2.6) $$\begin{eqnarray}q=Vw+2C^{\prime }V^{1/2}x^{1/2}.\end{eqnarray}$$

The combination of (2.3) and (2.6) leads to the Reynolds lubrication equation

(2.7) $$\begin{eqnarray}(Vw+2C^{\prime }V^{1/2}x^{1/2})^{n}=\frac{w^{2n+1}}{M^{\prime }}\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}x},\end{eqnarray}$$

which accounts for the fact that the pressure gradient is always positive for the problem under consideration.

Propagation criterion. The propagation condition imposes the asymptotic form of $w$ at the tip, which is given by the LEFM solution as (Rice Reference Rice and Liebowitz1968)

(2.8) $$\begin{eqnarray}w=\frac{K^{\prime }}{E^{\prime }}x^{1/2},\quad x\rightarrow 0.\end{eqnarray}$$

The propagation condition is complemented by zero flux condition $q=0$ at the fracture tip $x=0$ .

3 Non-singular formulation of the problem

To develop a computationally efficient numerical solution, we utilize the non-singular formulation, which was developed in Dontsov & Peirce (Reference Dontsov and Peirce2015b ) for Newtonian fluids. To construct the non-singular formulation, the elasticity equation (2.2) can be inverted (see Garagash & Detournay Reference Garagash and Detournay2000) and integrated by parts to yield

(3.1) $$\begin{eqnarray}w={\displaystyle \frac{K^{\prime }}{E^{\prime }}}x^{1/2}-{\displaystyle \frac{4}{\unicode[STIX]{x03C0}E^{\prime }}}\int _{0}^{\infty }F(x,s){\displaystyle \frac{\text{d}p}{\text{d}s}}\,\text{d}s,\quad F(x,s)=(s-x)\ln \left|{\displaystyle \frac{x^{1/2}+s^{1/2}}{x^{1/2}-s^{1/2}}}\right|-2x^{1/2}s^{1/2},\end{eqnarray}$$

where $E^{\prime }$ is the plane strain elastic modulus and $K^{\prime }$ is the scaled mode I fracture toughness of the material, as defined in (2.1).

By following the approach in Dontsov & Peirce (Reference Dontsov and Peirce2015b ), the pressure gradient is eliminated by combining the lubrication equation (2.7) and the elasticity equation (3.1), in which case the result can be rewritten in terms of a non-singular equation as

(3.2a,b ) $$\begin{eqnarray}\displaystyle \tilde{w}(\tilde{x})=1+{\displaystyle \frac{8}{\unicode[STIX]{x03C0}}}\int _{0}^{\infty }G\left(\frac{\tilde{s}}{\tilde{x}}\right){\displaystyle \frac{\tilde{s}^{1-n}}{\tilde{w}(\tilde{s})^{1+n}}}\left(1+\frac{\unicode[STIX]{x1D712}}{\tilde{w}(\tilde{s})}\right)^{n}\,\text{d}\tilde{s},\quad G(t)={\displaystyle \frac{1-t^{2}}{t}}\log \left|{\displaystyle \frac{1+t}{1-t}}\right|+2, & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

where the scaled quantities are defined as follows

(3.3) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\tilde{w}={\displaystyle \frac{E^{\prime }w}{K^{\prime }x^{1/2}}},\quad \unicode[STIX]{x1D712}={\displaystyle \frac{2C^{\prime }E^{\prime }}{V^{1/2}K^{\prime }}},\\ \tilde{x}=(x/l_{km})^{1/2},\quad \tilde{s}=(s/l_{km})^{1/2},\quad l_{km}=\left({\displaystyle \frac{{K^{\prime }}^{n+2}}{{E^{\prime }}^{1+n}M^{\prime }V^{n}}}\right)^{2/(2-n)}.\end{array}\right\} & & \displaystyle\end{eqnarray}$$

The kernel in the integral equation (3.2) is non-singular ( $0<G(t)\leqslant 4$ ) and is identically the same as for the Newtonian fluid case (Dontsov & Peirce Reference Dontsov and Peirce2015b ). The scaled fracture opening $\tilde{w}$ is always greater than unity, as can be clearly seen from (3.2). As a result, the only singular term in (3.2) is $\tilde{s}^{1-n}$ . However, this is a weak integrable singularity, which can be addressed by combining this term with the differential. It is important to note that the power-law exponent arises in both the governing equation (3.2) and in the scaling (3.3), and all the relations reduce to that for Newtonian fluids for $n=1$ .

4 Vertex solutions

It is instructive to outline the vertex solutions written in terms of the scaling (3.3) and their relation to (3.2). The toughness solution $\tilde{w}_{k}$ corresponds to the situation of negligible viscosity and leak-off, which can be obtained by setting $M^{\prime }\rightarrow 0$ , $C^{\prime }=0$ , so that the whole integral term in (3.2) can be neglected. The leak-off limit of the solution $\tilde{w}_{\tilde{m}}$ occurs for small toughness and large leak-off, i.e. $K^{\prime }=0$ , $C^{\prime }\rightarrow \infty$ . The latter situation can be realized in (3.2) by neglecting the unity in front of the integral and letting $\unicode[STIX]{x1D712}\gg 1$ . Finally, the viscous limit of the solution $\tilde{w}_{m}$ can be obtained by neglecting both the toughness and leak-off ( $K^{\prime }=0$ , $C^{\prime }=0$ ), in which case $\unicode[STIX]{x1D712}=0$ and the unity in front of the integral term in (3.2) should be neglected. To determine such solutions for the leak-off and viscosity limits, one needs to consider solution in the form $\tilde{w}=\unicode[STIX]{x1D6FD}\tilde{x}^{\unicode[STIX]{x1D6FF}}$ and use the relation

(4.1) $$\begin{eqnarray}{\mathcal{G}}(\unicode[STIX]{x1D6FC})={\displaystyle \frac{8}{\unicode[STIX]{x03C0}}}\int _{0}^{\infty }{\displaystyle \frac{G(t)}{t^{\unicode[STIX]{x1D6FC}}}}\,\text{d}t={\displaystyle \frac{16}{\unicode[STIX]{x1D6FC}(2-\unicode[STIX]{x1D6FC})}}\tan \left({\displaystyle \frac{\unicode[STIX]{x03C0}}{2}}\unicode[STIX]{x1D6FC}\right),\quad -1<\unicode[STIX]{x1D6FC}<1,\end{eqnarray}$$

to obtain

(4.2a-c ) $$\begin{eqnarray}\tilde{w}_{k}=1,\quad \tilde{w}_{\tilde{m}}=\unicode[STIX]{x1D6FD}_{\tilde{m}}\unicode[STIX]{x1D712}^{n/(2(1+n))}\tilde{x}^{(2-n)/(2(1+n))},\quad \tilde{w}_{m}=\unicode[STIX]{x1D6FD}_{m}\tilde{x}^{(2-n)/(2+n)},\end{eqnarray}$$

where the constants $\unicode[STIX]{x1D6FD}_{\tilde{m}}$ and $\unicode[STIX]{x1D6FD}_{m}$ are given respectively by

(4.3a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D6FD}_{\tilde{m}}=\left[{\displaystyle \frac{64(1+n)^{2}}{3n(4+n)}}\tan \left(\frac{3\unicode[STIX]{x03C0}n}{4(1+n)}\right)\right]^{1/(2(1+n))},\quad \unicode[STIX]{x1D6FD}_{m}=\left[{\displaystyle \frac{2(2+n)^{2}}{n}}\tan \left(\frac{\unicode[STIX]{x03C0}n}{2+n}\right)\right]^{1/(2+n)}.\end{eqnarray}$$

Note that $\unicode[STIX]{x1D6FD}_{m}=\unicode[STIX]{x1D6FD}_{\tilde{m}}=\sqrt{4\unicode[STIX]{x03C0}}$ for $n=0$ and the unscaled expressions for (4.2) are

(4.4) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle w_{k}={\displaystyle \frac{K^{\prime }}{E^{\prime }}}x^{1/2},\quad w_{\tilde{m}}=\unicode[STIX]{x1D6FD}_{\tilde{m}}\left[{\displaystyle \frac{(2C^{\prime })^{n}V^{n/2}M^{\prime }}{E^{\prime }}}\right]^{1/(2(1+n))}x^{(4+n)/4(1+n)},\\ \displaystyle \quad w_{m}=\unicode[STIX]{x1D6FD}_{m}\left({\displaystyle \frac{V^{n}M^{\prime }}{E^{\prime }}}\right)^{1/(2+n)}x^{2/(2+n)}.\end{array}\right\} & & \displaystyle\end{eqnarray}$$

It is interesting to observe that for $w_{m}=w_{\tilde{m}}\propto x$ for $n\rightarrow 0$ . At the same time, $w_{m}\propto x^{1/2}$ and $w_{\tilde{m}}\propto x^{1/2}$ for $n\rightarrow 2$ . Note that for the latter case $\unicode[STIX]{x1D6FD}_{m}\rightarrow \infty$ and $\unicode[STIX]{x1D6FD}_{\tilde{m}}\rightarrow \infty$ , which is not physical. This happens because the viscous dissipation integrated over the tip region becomes singular, therefore there must be a fluid lag to prevent the singularity. Situations with $n\approx 2$ are beyond the scope of this study and are left for future investigations, especially since $n\leqslant 1$ in typical practical applications. Note that the limiting solutions (4.4) were obtained earlier in Desroches et al. (Reference Desroches, Detournay, Lenoach, Papanastasiou, Pearson, Thiercelin and Cheng1994), Lenoach (Reference Lenoach1995).

5 Solution inside $m\tilde{m}k$ parametric space

Figure 1(b) shows the triangular parametric space $m\tilde{m}k$ , which is a schematic representation of the parametric space of the solutions, see e.g. Garagash et al. (Reference Garagash, Detournay and Adachi2011) for a similar representation for Newtonian fluids. Each vertex corresponds to one of the limiting solutions (4.3), and the global solution transitions from one vertex to another with increasing distance from the fracture tip. In order to obtain the solution inside the parametric space, this section presents the numerical solution of the integral equation (3.2). As for the Newtonian fluid case (Dontsov & Peirce Reference Dontsov and Peirce2015b ), the non-singular formulation allows us to use standard numerical techniques to obtain accurate numerical solution. In particular, the integral in (3.2) is discretized using Simpson’s rule and the resulting system of nonlinear algebraic equations is solved iteratively using Newton’s method. Further details of the numerical scheme are omitted for brevity.

Figure 2. Numerical solution for the normalized fracture width versus $\tilde{x}$ and $\unicode[STIX]{x1D712}$ for $n=\{0,0.5,1,1.5\}$ , (a–d) respectively. Red, blue and green regions respectively indicate the applicability of the toughness ( $k$ ), viscous ( $m$ ) and leak-off ( $\tilde{m}$ ) vertex solutions defined in (4.2a-c ).

To illustrate the behaviour of the solution for different values of $n$ , figure 2 plots the variation of the normalized fracture width $\tilde{w}$ (calculated numerically from (3.2)) versus $\tilde{x}$ and $\unicode[STIX]{x1D712}$ for $n=\{0,0.5,1,1.5\}$ . Coloured zones indicate applicability zones of the vertex solutions (4.2), where the red coloured region corresponds to the toughness solution (indicated by the $k$ symbol), the blue coloured region stands for the viscous solution (indicated by the $m$ symbol) while the green coloured region shows the leak-off solution (indicated by the $\tilde{m}$ symbol). The applicability regions are defined as the zones in which the relative difference between the numerical solution and the corresponding vertex solution does not exceed 1 %, i.e.

(5.1) $$\begin{eqnarray}\displaystyle {\displaystyle \frac{\tilde{w}-\tilde{w}_{i}}{\tilde{w}}}=\unicode[STIX]{x1D716},\quad i=k,m,\tilde{m}, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D716}=0.01$ , $\tilde{w}$ is the numerical solution and $\tilde{w}_{i}$ is one of the vertex solutions (4.2a-c ). Black solid lines in figure 2 indicate the boundaries of the vertex solution applicability zones. There is no leak-off limit for $n=0$ , which is consistent with the governing equation (3.2), since it does not depend on $\unicode[STIX]{x1D712}$ in this case. Note that $\tilde{w}_{m}\equiv \tilde{w}_{\tilde{m}}$ for $n=0$ , in which case the $m\tilde{m}k$ triangular parametric space collapses to the $mk$ edge. The leak-off zone appears for $n>0$ and increases in size with respect to $\tilde{x}$ and moves towards smaller values of $\tilde{x}$ with the rise of $n$ . It is also interesting to observe the widening of the transition regions, which occurs in the normalized coordinate space $(\tilde{x},\unicode[STIX]{x1D712})$ . Also note that the limiting solutions (4.2) correspond to planes on a logarithmic scale in figure 2 so that the red, the blue and the green regions are planar. Therefore, geometrically, the solution transitions gradually from one plane to another and the transitions occur in the vicinity of plane intersections.

To illustrate spatial variation of the solution, figure 3(a) plots the normalized fracture width $\tilde{w}$ versus spatial coordinate $\tilde{x}$ for $\unicode[STIX]{x1D712}=0$ and $n=\{0,0.25,0.5,0.75,1,1.25,1.5,175\}$ . Figure 3(b) shows similar solution but for larger leak-off $\unicode[STIX]{x1D712}=10^{2}$ . Dashed red, green and blue lines indicate respectively the toughness, leak-off and viscous vertex solutions (4.2). Note that to reduce overlapping of multiple curves, all the vertex solutions are plotted at fixed spatial locations rather than their respective applicability regions. One can observe from the figure 3 that solutions transition from the toughness limiting solution to the viscous limiting solution for $\unicode[STIX]{x1D712}=0$ , and from toughness, to leak-off, and then to viscous solutions for $\unicode[STIX]{x1D712}=10^{2}$ . Transition regions become wider and the solution becomes ‘flatter’ for $n$ approaching $2$ , which is consistent with the results shown in figure 2.

Figure 3. Numerical solution for the normalized fracture width versus $\tilde{x}$ for $n=\{0,0.25,0.5,0.75,1,1.25,1.5,175\}$ and $\unicode[STIX]{x1D712}=0$ (a), $\unicode[STIX]{x1D712}=10^{2}$ (b). Red, blue and green dashed lines respectively indicate the toughness ( $k$ ), viscous ( $m$ ) and leak-off ( $\tilde{m}$ ) vertex solutions defined in (4.2).

6 Transition regions

6.1 Scaling of the transitions

There are three relevant transitions that determine the asymptotic bounds: the $km$ transition, the $k\tilde{m}$ transition and the $\tilde{m}m$ transition. The $km$ transition describes the transition of the solution from the toughness ( $k$ ) solution to the viscous ( $m$ ) solution. Similarly, the $k\tilde{m}$ and $\tilde{m}m$ transitions capture the evolution from $k$ to $\tilde{m}$ and $\tilde{m}$ to $m$ vertex solutions. Here the notation is such that the first letter corresponds to the vertex solution that applies closer to the fracture tip and the second letter represents the vertex solution that works further away from the fracture tip.

The $km$ transition occurs in the limit of no leak-off and is governed by

(6.1) $$\begin{eqnarray}\tilde{w}(\tilde{x})=1+{\displaystyle \frac{8}{\unicode[STIX]{x03C0}}}\int _{0}^{\infty }G\left(\frac{\tilde{s}}{\tilde{x}}\right){\displaystyle \frac{\tilde{s}^{1-n}}{\tilde{w}(\tilde{s})^{1+n}}}\,\text{d}\tilde{s},\end{eqnarray}$$

which is obtained from (3.2) by setting $\unicode[STIX]{x1D712}=0$ . Equation (6.1) has no independent parameters, in which case the $km$ transition is determined solely by the spatial coordinate

(6.2) $$\begin{eqnarray}\tilde{x}_{km}=\tilde{x}=\left[{\displaystyle \frac{{E^{\prime }}^{1+n}V^{n}M^{\prime }}{{K^{\prime }}^{n+2}}}\right]^{1/(2-n)}x^{1/2}=\left({\displaystyle \frac{x}{l_{km}}}\right)^{1/2},\end{eqnarray}$$

see figure 2(d).

The $k\tilde{m}$ transition occurs for large values of leak-off, i.e. $\unicode[STIX]{x1D712}\gg 1$ , in which case (3.2) can be reduced to

(6.3a-c ) $$\begin{eqnarray}\tilde{w}(\tilde{x}_{k\tilde{m}})=1+{\displaystyle \frac{8}{\unicode[STIX]{x03C0}}}\int _{0}^{\infty }G\left(\frac{\tilde{s}_{k\tilde{m}}}{\tilde{x}_{k\tilde{m}}}\right){\displaystyle \frac{\tilde{s}_{k\tilde{m}}^{1-n}}{\tilde{w}(\tilde{s}_{k\tilde{m}})^{1+2n}}}\,\text{d}\tilde{s}_{k\tilde{m}},\quad \tilde{x}_{k\tilde{m}}=\unicode[STIX]{x1D712}^{n/(2-n)}\tilde{x},\quad \tilde{s}_{k\tilde{m}}=\unicode[STIX]{x1D712}^{n/(2-n)}\tilde{s},\end{eqnarray}$$

where an additional scaling is performed to eliminate the explicit dependence on $\unicode[STIX]{x1D712}$ . As a result, the parameter that describes the $k\tilde{m}$ transition is

(6.4) $$\begin{eqnarray}\tilde{x}_{k\tilde{m}}=\unicode[STIX]{x1D712}^{n/(2-n)}\tilde{x}=\left[{\displaystyle \frac{(2C^{\prime })^{n}{E^{\prime }}^{1+2n}V^{n/2}M^{\prime }}{{K^{\prime }}^{2n+2}}}\right]^{1/(2-n)}x^{1/2}=\left({\displaystyle \frac{x}{l_{k\tilde{m}}}}\right)^{1/2},\end{eqnarray}$$

which defines the transition length scale $l_{k\tilde{m}}$ . Figure 2(d) indicates lines in the parametric space that correspond to fixed values of the parameter $\tilde{x}_{k\tilde{m}}$ .

To reduce the governing equation (3.2) to the $\tilde{m}m$ transition, it is necessary to disregard the contribution of the fracture toughness to the solution. This can be achieved by neglecting the constant term (which is equal to one) that appears outside of the integral, in which case equation (3.2) can be written as

(6.5) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle \tilde{w}_{\tilde{m}m}(\tilde{x}_{\tilde{m}m})={\displaystyle \frac{8}{\unicode[STIX]{x03C0}}}\int _{0}^{\infty }G\left(\frac{\tilde{s}_{\tilde{m}m}}{\tilde{x}_{\tilde{m}m}}\right){\displaystyle \frac{\tilde{s}_{\tilde{m}m}^{1-n}}{\tilde{w}_{\tilde{m}m}(\tilde{s}_{\tilde{m}m})^{1+n}}}\left(1+\frac{1}{\tilde{w}_{\tilde{m}m}(\tilde{s}_{\tilde{m}m})}\right)^{n}\,\text{d}\tilde{s}_{\tilde{m}m},\\ \displaystyle \tilde{w}_{\tilde{m}m}=\unicode[STIX]{x1D712}^{-1}\tilde{w},\quad \tilde{x}_{m\tilde{m}}=\unicode[STIX]{x1D712}^{-((2+n)/(2-n))}\tilde{x},\quad \tilde{s}_{m\tilde{m}}=\unicode[STIX]{x1D712}^{-((2+n)/(2-n))}\tilde{s},\end{array}\right\} & & \displaystyle\end{eqnarray}$$

where an additional scaling is used for both the width and the spatial coordinate to remove the explicit dependence on the parameter $\unicode[STIX]{x1D712}$ . The parameter that describes the $\tilde{m}m$ transition is then

(6.6) $$\begin{eqnarray}\tilde{x}_{m\tilde{m}}=\unicode[STIX]{x1D712}^{-((2+n)/(2-n))}\tilde{x}=\left[{\displaystyle \frac{V^{1+3n/2}M^{\prime }}{(2C^{\prime })^{2+n}E^{\prime }}}\right]^{1/(2-n)}x^{1/2}=\left({\displaystyle \frac{x}{l_{m\tilde{m}}}}\right)^{1/2},\end{eqnarray}$$

which can also be used as the definition of the length scale $l_{m\tilde{m}}$ . Lines that correspond to constant values of the parameter $\tilde{x}_{m\tilde{m}}$ are indicated in figure 2 $(d)$ .

6.2 Numerical calculation of the bounds

Each transition region can be characterized by two numbers which indicate the beginning (denoted by $\tilde{x}^{(1)}$ ) and the end (denoted by $\tilde{x}^{(2)}$ ) of the transition. In this study, these numbers are determined by the requirement that the difference between the numerical solution and the corresponding vertex solution is equal to 1 % (5.1), which is consistent with the black lines shown in figure 2. Figure 4 shows variation of the numerically computed boundaries of all three transition zones $km$ , $k\tilde{m}$ , and $m\tilde{m}$ versus $n$ . Note that the abscissa axis is $(2-n)\log _{10}\tilde{x}$ , which accounts the singular behaviour at $n=2$ . To have an ability to reconstruct the boundaries of applicability of the vertex solutions (such as shown in figure 2) for any $n$ , the data points are fitted to the second degree polynomial (indicated by dashed lines in figure 4)

(6.7) $$\begin{eqnarray}\log _{10}\tilde{x}=(2-n)^{-1}(p_{1}n^{2}+p_{2}n+p_{3}),\end{eqnarray}$$

where the multiplier $2-n$ captures the singular behaviour near $n=2$ . The coefficients $p_{1}$ , $p_{2}$ , and $p_{3}$ for each transition are summarized in the table 1. Equation (6.7), the data in table 1, and the scaling relations (6.2), (6.4) and (6.6) can be used to reconstruct zones of applicability of the vertex solution inside the whole parametric space (the black lines shown in figure 2(d)) for any value of $n$ .

Figure 4. Boundaries of the $km$ , $k\tilde{m}$ and $m\tilde{m}$ transitions versus $n$ . Quantities with the superscripts ‘ $(1)$ ’ and ‘ $(2)$ ’ correspond to the beginning and end of the transitions. Dashed lines represent the curves that are calculated based on the fitting (6.7).

Figure 5. Validity regions of the vertex solutions that are calculated using asymptotic expansions (6.8) and (6.9) (blue solid, dashed and dotted lines) for different values of the power-law exponent $n$ . The black dashed lines indicate the bounds calculated using the equation and values in table 1.

Table 1. Coefficients of the second degree polynomial (6.7) that is used to approximate boundaries of the transitions.

6.3 Calculation of the bounds from asymptotic expansions

It is also possible to calculate bounds of applicability of the vertex solutions using an alternative technique that involves asymptotic expansions. The procedure is similar to that for Newtonian fluids, see (Garagash et al. Reference Garagash, Detournay and Adachi2011), but requires overcoming complexities that are associated with the power-law rheology. As a result, only a brief outline of the procedure is presented, while cumbersome mathematical calculations and lengthy expressions for the higher-order terms are omitted.

The vertex solutions (4.2) represent the leading behaviour of the solution near the corresponding vertices. Asymptotic analysis of the governing equations (2.2), (2.7) and (2.8) or simply (3.2) can be used to calculate the respective second terms in the expansions, i.e.

(6.8) $$\begin{eqnarray}\tilde{w}(\tilde{x},\unicode[STIX]{x1D712})\approx \tilde{w}_{i}(\tilde{x},\unicode[STIX]{x1D712})+\tilde{w}_{i,2}(\tilde{x},\unicode[STIX]{x1D712}),\quad i=k,m,\tilde{m},\end{eqnarray}$$

where $\tilde{w}_{i}$ represents the vertex solution and $\tilde{w}_{i,2}$ is the second term in the expansion. The bounds for each vertex are then defined similarly to (5.1) via the relations

(6.9) $$\begin{eqnarray}{\displaystyle \frac{\tilde{w}_{i,2}(\tilde{x},\unicode[STIX]{x1D712})}{\tilde{w}_{i}(\tilde{x},\unicode[STIX]{x1D712})}}=\unicode[STIX]{x1D716},\quad i=k,m,\tilde{m},\end{eqnarray}$$

where $\unicode[STIX]{x1D716}=0.01$ is taken in the calculations. Figure 5 plots zones of applicability of the vertex solutions that are calculated using the asymptotic expansion (blue solid, dashed and dotted lines) for different values of $n$ , as well as the zones of applicability that are evaluated using (6.7) and the results in table 1 for $n=1$ (back dashed lines). As expected, both approaches give nearly identical results.

Despite the approach with the asymptotic expansion is more rigorous, it is less practical. The reason for this lies in the fact that some of the coefficients in the expansion can only be computed numerically. This was also the case for Newtonian fluids, see the related discussion in Garagash et al. (Reference Garagash, Detournay and Adachi2011). The unknown coefficients can be computed once for Newtonian fluids, while for power-law fluids they become functions of $n$ . These functions can be computed numerically and then approximated or interpolated. In this situation, however, the method loses its primary advantage of being independent of the numerical solution, so that one may use the direct numerical computation of the applicability bounds, as presented in section 6.2.

7 Approximate solution

Despite the fact that the numerical solution for the given problem of a semi-infinite hydraulic fracture can always be obtained, there are situations when rapid evaluation of the solution is necessary. As shown for the Newtonian fluid case, such a rapid solution can be used as a propagation condition in a hydraulic fracturing simulator for a planar or multiple planar fractures (Dontsov & Peirce Reference Dontsov and Peirce2016, Reference Dontsov and Peirce2017), and to construct rapid solutions for plane strain and radial hydraulic fracture geometries (Dontsov Reference Dontsov2016a , Reference Dontsov2017).

To develop an approximate solution for a semi-infinite hydraulic fracture driven by a power-law fluid, assume that $\tilde{w}\propto \tilde{x}^{\tilde{\unicode[STIX]{x1D6FF}}_{1}}$ and $(\tilde{w}+\unicode[STIX]{x1D712})\propto \tilde{x}^{\tilde{\unicode[STIX]{x1D6FF}}_{2}}$ , in which case the integral equation (3.2) can be differentiated and simplified to

(7.1) $$\begin{eqnarray}{\displaystyle \frac{\text{d}\tilde{w}}{\text{d}\tilde{x}}}=[C_{1}(\tilde{\unicode[STIX]{x1D6FF}}_{1})-C_{3}(\tilde{\unicode[STIX]{x1D6FF}}_{1})]{\displaystyle \frac{\tilde{x}^{1-n}}{\tilde{w}(\tilde{x})^{1+n}}}+C_{2}(\tilde{\unicode[STIX]{x1D6FF}}_{1},\tilde{\unicode[STIX]{x1D6FF}}_{2}){\displaystyle \frac{\tilde{x}^{1-n}}{\tilde{w}(\tilde{x})^{1+n}}}\left(1+\frac{\unicode[STIX]{x1D712}}{\tilde{w}(\tilde{x})}\right)^{n},\quad \tilde{w}(0)=1,\end{eqnarray}$$

where $\tilde{\unicode[STIX]{x1D6FF}}_{1}=(\tilde{x}/\tilde{w})\,\text{d}\tilde{w}/\text{d}\tilde{x}$ and $\tilde{\unicode[STIX]{x1D6FF}}_{2}=(\tilde{x}/(\tilde{w}+\unicode[STIX]{x1D712}))\,\text{d}\tilde{w}/\text{d}\tilde{x}$ . Here the coefficients $C_{1}$ , $C_{2}$ and $C_{3}$ are calculated as

(7.2a ) $$\begin{eqnarray}\displaystyle & \displaystyle C_{1}(\tilde{\unicode[STIX]{x1D6FF}}_{1})=C_{3}(\tilde{\unicode[STIX]{x1D6FF}}_{1})={\mathcal{C}}(2-n-(1+n)\tilde{\unicode[STIX]{x1D6FF}}_{1}), & \displaystyle\end{eqnarray}$$

(7.2b ) $$\begin{eqnarray}\displaystyle & \displaystyle \quad C_{2}(\tilde{\unicode[STIX]{x1D6FF}}_{1},\tilde{\unicode[STIX]{x1D6FF}}_{2})={\mathcal{C}}(2-n-(1+2n)\tilde{\unicode[STIX]{x1D6FF}}_{1}+n\tilde{\unicode[STIX]{x1D6FF}}_{2}), & \displaystyle\end{eqnarray}$$

(7.3) $$\begin{eqnarray}{\mathcal{C}}(\unicode[STIX]{x1D6FE})=-{\displaystyle \frac{8}{\unicode[STIX]{x03C0}}}\int _{0}^{\infty }G^{\prime }(t)t^{\unicode[STIX]{x1D6FE}}\,\text{d}t={\displaystyle \frac{16\unicode[STIX]{x1D6FE}}{1-\unicode[STIX]{x1D6FE}^{2}}}\tan \left({\displaystyle \frac{\unicode[STIX]{x03C0}}{2}}(1-\unicode[STIX]{x1D6FE})\right).\end{eqnarray}$$

The multiplier $C_{1}$ captures the ‘viscous’ term ( $C_{2}\equiv C_{3}$ for $\unicode[STIX]{x1D712}=0$ since $\tilde{\unicode[STIX]{x1D6FF}}_{2}=\tilde{\unicode[STIX]{x1D6FF}}_{1}$ ), while $C_{2}$ and $C_{3}$ account for the ‘leak-off’ term. If one uses accurate expressions for $\tilde{\unicode[STIX]{x1D6FF}}_{1}$ and $\tilde{\unicode[STIX]{x1D6FF}}_{2}$ , then $C_{1}\equiv C_{3}$ , which greatly simplifies the governing equation (7.1). At the same time, for the purpose of constructing the approximate solution, it is beneficial to use constant values for the parameters $C_{1}$ , $C_{2}$ and $C_{3}$ . To ensure that the solution of ordinary differential equation (7.1) precisely captures the limiting solutions (4.2), one can use the following constant approximations for the parameters $C_{1}$ , $C_{2}$ and $C_{3}$

(7.4a-d ) $$\begin{eqnarray}C_{1}={\mathcal{C}}(\unicode[STIX]{x1D6FF}_{m})=\unicode[STIX]{x1D6FF}_{m}\unicode[STIX]{x1D6FD}_{m}^{2+n},\quad C_{2}=C_{3}={\mathcal{C}}(\unicode[STIX]{x1D6FF}_{\tilde{m}})=\unicode[STIX]{x1D6FF}_{\tilde{m}}\unicode[STIX]{x1D6FD}_{\tilde{m}}^{2(1+n)},\quad \unicode[STIX]{x1D6FF}_{m}={\displaystyle \frac{2-n}{2+n}},\quad \unicode[STIX]{x1D6FF}_{\tilde{m}}={\displaystyle \frac{2-n}{2(1+n)}},\end{eqnarray}$$

where $\unicode[STIX]{x1D6FD}_{m}$ and $\unicode[STIX]{x1D6FD}_{\tilde{m}}$ are defined in (4.3a,b ). To develop a more accurate approximation, it is useful to consider a more general case, in which the whole differential equation (7.1) (with constant approximations for $C_{1}$ , $C_{2}$ and $C_{3}$ defined in (7.4)) is multiplied by $\tilde{\unicode[STIX]{x1D6FF}}_{1}^{\unicode[STIX]{x1D703}}$ (where $\unicode[STIX]{x1D703}$ is a small parameter), which leads to

(7.5) $$\begin{eqnarray}\displaystyle \left({\displaystyle \frac{\tilde{x}}{\tilde{w}}}{\displaystyle \frac{\text{d}\tilde{w}}{\text{d}\tilde{x}}}\right)^{\unicode[STIX]{x1D703}}{\displaystyle \frac{\text{d}\tilde{w}}{\text{d}\tilde{x}}}=[\unicode[STIX]{x1D6FF}_{m}^{1+\unicode[STIX]{x1D703}}\unicode[STIX]{x1D6FD}_{m}^{2+n}-\unicode[STIX]{x1D6FF}_{\tilde{m}}^{1+\unicode[STIX]{x1D703}}\unicode[STIX]{x1D6FD}_{\tilde{m}}^{2(1+n)}]{\displaystyle \frac{\tilde{x}^{1-n}}{\tilde{w}(\tilde{x})^{1+n}}}+\unicode[STIX]{x1D6FF}_{\tilde{m}}^{1+\unicode[STIX]{x1D703}}\unicode[STIX]{x1D6FD}_{\tilde{m}}^{2(1+n)}{\displaystyle \frac{\tilde{x}^{1-n}}{\tilde{w}(\tilde{x})^{1+n}}}\left(1+\frac{\unicode[STIX]{x1D712}}{\tilde{w}(\tilde{x})}\right)^{n}. & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

Here the definition of $\tilde{\unicode[STIX]{x1D6FF}}_{1}$ is used for the left-hand side of the equation, while $\tilde{\unicode[STIX]{x1D6FF}}_{1}=\unicode[STIX]{x1D6FF}_{m}$ for the viscous term and $\tilde{\unicode[STIX]{x1D6FF}}_{1}=\unicode[STIX]{x1D6FF}_{\tilde{m}}$ for the leak-off terms on the right-hand side of the equation, which ensures that the modified differential equation still captures the vertex solutions precisely for any value of the parameter $\unicode[STIX]{x1D703}$ . With the use of the initial condition $\tilde{w}(0)=1$ , the differential equation (7.5) can be integrated to obtain

(7.6) $$\begin{eqnarray}\displaystyle & & \displaystyle \int _{1}^{\tilde{w}}\left\{\unicode[STIX]{x1D6FF}_{m}^{1+\unicode[STIX]{x1D703}}\unicode[STIX]{x1D6FD}_{m}^{2+n}+\unicode[STIX]{x1D6FF}_{\tilde{m}}^{1+\unicode[STIX]{x1D703}}\unicode[STIX]{x1D6FD}_{\tilde{m}}^{2(1+n)}\left[\left(1+\frac{\unicode[STIX]{x1D712}}{\bar{w}}\right)^{n}-1\right]\right\}^{-(1/(1+\unicode[STIX]{x1D703}))}\bar{w}^{(1+n-\unicode[STIX]{x1D703})/(1+\unicode[STIX]{x1D703})}\,\text{d}\bar{w}\nonumber\\ \displaystyle & & \displaystyle \quad ={\displaystyle \frac{1+\unicode[STIX]{x1D703}}{2-n}}\tilde{x}^{(2-n)/(1+\unicode[STIX]{x1D703})}.\end{eqnarray}$$

The integral in (7.6) can be calculated analytically for $\unicode[STIX]{x1D712}=0$ and $\unicode[STIX]{x1D712}\gg 1$ . This result can be used to construct the approximation, which reduces (7.6) to

(7.7) $$\begin{eqnarray}\tilde{x}^{(2-n)/(1+\unicode[STIX]{x1D703})}=\tilde{\unicode[STIX]{x1D6FF}}_{1}(\tilde{w},\unicode[STIX]{x1D712})\,\tilde{w}^{(2+n)/(1+\unicode[STIX]{x1D703})}\left\{\unicode[STIX]{x1D6FF}_{m}^{1+\unicode[STIX]{x1D703}}\unicode[STIX]{x1D6FD}_{m}^{2+n}+\unicode[STIX]{x1D6FF}_{\tilde{m}}^{1+\unicode[STIX]{x1D703}}\unicode[STIX]{x1D6FD}_{\tilde{m}}^{2(1+n)}\left[\left(1+{\displaystyle \frac{\unicode[STIX]{x1D712}}{\tilde{w}}}\right)^{n}-1\right]\!\right\}^{-(1/(1+\unicode[STIX]{x1D703}))},\end{eqnarray}$$

where the parameter $\tilde{\unicode[STIX]{x1D6FF}}_{1}$ is determined from the relation $\tilde{\unicode[STIX]{x1D6FF}}_{1}=(\tilde{x}/\tilde{w})d\tilde{w}/d\tilde{x}$ (since $\tilde{w}\propto \tilde{x}^{\tilde{\unicode[STIX]{x1D6FF}}_{1}}$ ) as

(7.8) $$\begin{eqnarray}\displaystyle \tilde{\unicode[STIX]{x1D6FF}}_{1}(\tilde{w},\unicode[STIX]{x1D712}\!) & = & \displaystyle {\displaystyle \frac{\unicode[STIX]{x1D6FF}_{\tilde{m}}\unicode[STIX]{x1D6FF}_{m}V_{\tilde{m}}V_{m}\left[\!\unicode[STIX]{x1D6FD}_{m}^{(2+n)/n}V_{\tilde{m}}^{(1+\unicode[STIX]{x1D703})/n}+{\displaystyle \frac{\unicode[STIX]{x1D712}}{\tilde{w}}}\unicode[STIX]{x1D6FD}_{\tilde{m}}^{(2(1+n))/n}V_{m}^{(1+\unicode[STIX]{x1D703})/n}\!\right]}{\unicode[STIX]{x1D6FF}_{\tilde{m}}V_{\tilde{m}}\unicode[STIX]{x1D6FD}_{m}^{(2+n)/n}V_{\tilde{m}}^{(1+\unicode[STIX]{x1D703})/n}+\unicode[STIX]{x1D6FF}_{m}V_{m}{\displaystyle \frac{\unicode[STIX]{x1D712}}{\tilde{w}}}\unicode[STIX]{x1D6FD}_{\tilde{m}}^{(2(1+n))/n}V_{m}^{(1+\unicode[STIX]{x1D703})/n}}},\quad \nonumber\\ \displaystyle & & \displaystyle \left.\begin{array}{@{}l@{}}V_{m}=1-\tilde{w}^{-((2+n)/(1+\unicode[STIX]{x1D703}))},\\ V_{\tilde{m}}=1-\tilde{w}^{-((2+2n)/(1+\unicode[STIX]{x1D703}))}.\end{array}\right\}\end{eqnarray}$$

Figure 6. (a) Maximum error of the approximate solution versus parameter $\unicode[STIX]{x1D703}$ for the case $n=0$ (results do not depend on $\unicode[STIX]{x1D712}$ for $n=0$ ). (b) Optimal values of $\unicode[STIX]{x1D703}$ versus $n$ for different values of the leak-off parameter $\unicode[STIX]{x1D712}$ .

To find the optimal values of $\unicode[STIX]{x1D703}$ that produce the most accurate approximation, the maximum error between the numerical solution (for given values of $n$ and $\unicode[STIX]{x1D712}$ ) and the approximation is first plotted against $\unicode[STIX]{x1D703}$ , see figure 6(a). Such plots always have a local minimum, which corresponds to the optimal value of $\unicode[STIX]{x1D703}=\unicode[STIX]{x1D703}_{opt}$ . Figure 6(b) shows the calculated optimal values of $\unicode[STIX]{x1D703}_{opt}$ versus $n$ for $\unicode[STIX]{x1D712}=\{0,1,10^{2},10^{3}\}$ . Clearly, the variation of $\unicode[STIX]{x1D703}_{opt}$ versus $\unicode[STIX]{x1D712}$ is small. As a result, it is possible to use $\unicode[STIX]{x1D712}=0$ data to determine $\unicode[STIX]{x1D703}_{opt}(n)$ . By fitting the $\unicode[STIX]{x1D712}=0$ data to a second degree polynomial, the expression for the parameter $\unicode[STIX]{x1D703}$ becomes

(7.9) $$\begin{eqnarray}\unicode[STIX]{x1D703}=\unicode[STIX]{x1D703}_{opt}(n)=0.0452\,n^{2}-0.1780\,n+0.1753.\end{eqnarray}$$

Finally, the approximate solution $\tilde{w}(\tilde{x},\unicode[STIX]{x1D712},n)$ for any $0\leqslant n\leqslant 1.9$ is implicitly determined by equations (7.7) and (7.8), where the values of $\unicode[STIX]{x1D6FD}_{m}$ , $\unicode[STIX]{x1D6FD}_{\tilde{m}}$ , $\unicode[STIX]{x1D6FF}_{m}$ , and $\unicode[STIX]{x1D6FF}_{\tilde{m}}$ are given in (4.3) and (7.4), while the parameter $\unicode[STIX]{x1D703}$ is specified in (7.9).

It is interesting to observe that the equations (7.7) and (7.8) reduce to much simpler expressions for the $km$ and $k\tilde{m}$ transitions. In particular, there is no leak-off for the $km$ transition (i.e. $\unicode[STIX]{x1D712}=0$ ), so that $\tilde{\unicode[STIX]{x1D6FF}}_{1}(\tilde{w},\unicode[STIX]{x1D712})=\unicode[STIX]{x1D6FF}_{m}V_{m}$ and the solution becomes

(7.10) $$\begin{eqnarray}\tilde{w}_{km}=\left[1+(\unicode[STIX]{x1D6FD}_{m}^{2+n}\tilde{x}^{2-n})^{1/(1+\unicode[STIX]{x1D703})}\right]^{(1+\unicode[STIX]{x1D703})/(2+n)},\quad \unicode[STIX]{x1D712}=0.\end{eqnarray}$$

For the $k\tilde{m}$ transition, on the other hand, there is large leak-off (i.e. $\unicode[STIX]{x1D712}\gg 1$ ), in which case $\tilde{\unicode[STIX]{x1D6FF}}_{1}(\tilde{w},\unicode[STIX]{x1D712})=\unicode[STIX]{x1D6FF}_{\tilde{m}}V_{\tilde{m}}$ and the solution is

(7.11) $$\begin{eqnarray}\tilde{w}_{k\tilde{m}}=\left[1+(\unicode[STIX]{x1D6FD}_{\tilde{m}}^{2(1+n)}\unicode[STIX]{x1D712}^{n}\tilde{x}^{2-n})^{1/(1+\unicode[STIX]{x1D703})}\right]^{(1+\unicode[STIX]{x1D703})/(2(1+n))},\quad \unicode[STIX]{x1D712}\gg 1.\end{eqnarray}$$

As can be seen from (7.10) and (7.11), the solution evolves from the $k$ vertex to either $m$ or $\tilde{m}$ vertex as the scaled distance $\tilde{x}$ increases, see (4.2) for the vertex solutions. The unscaled solutions can be easily obtained by using the scaling relations in (3.3).

Figure 7. Relative difference between the numerical solution and the approximate solution for the normalized fracture width versus $\tilde{x}$ and $\unicode[STIX]{x1D712}$ for $n=\{0,0.5,1,1.5\}$ , (a–d) respectively.

Figure 8. Maximum relative error of the approximate solution versus $n$ .

To quantify the error of the approximate solution, figure 7 shows the relative error between the numerical solution of (3.2) and the approximation versus $\tilde{x}$ and $\unicode[STIX]{x1D712}$ for $n=\{0,0.5,1,1.5\}$ . It can be clearly seen that the error is concentrated in the transition regions, and reaches its maximum at the intersection of the three transition regions. Figure 8 shows variation of the maximum error versus $n$ , where the maximum error is computed over the range of $\tilde{x}$ and $\unicode[STIX]{x1D712}$ that captures all the transitions. The largest error of approximately 1 % is observed for $n$ in the range between $0.1$ and $0.25$ . The error then decreases for larger values of $n$ and for $n=0$ . These results demonstrate that the developed approximate solution is able to accurately and uniformly approximate the solution for a semi-infinite hydraulic fracture with leak-off that is driven by a power-law fluid. The approximation is given implicitly by a relatively simple expression (7.7), which permits one to obtain the solution rapidly. This is especially important for hydraulic fracturing simulators that utilize the tip asymptotic solution for the propagation condition, see e.g. Peirce & Detournay (Reference Peirce and Detournay2008), Peirce (Reference Peirce2015), Dontsov & Peirce (Reference Dontsov and Peirce2017).

8 Assumptions and validity of the model

The mathematical model for the semi-infinite hydraulic fracture, outlined in § 2, relies on a series of assumptions that lead to a relatively simple problem formulation. The aim of this section is to revisit some of these assumptions and to analyse regions of validity of the developed solution. Readers should be warned that some typical parameters are used to justify the assumptions. At the same time, there can be situations, in which the model assumptions will not hold. Therefore, readers are strongly encouraged to double check validity of the model for their particular problem parameters before using the results of this paper.

8.1 Fracture process zone

One of the assumptions of the model arises in using the linear elastic fracture mechanics (LEFM) for the fracture propagation criterion. By doing this, it is implicitly assumed that the size of the plastic zone around the crack tip is much smaller than the distances at which the solution is computed, i.e. $x\gg x_{p}$ . The size of the plastic zone can be estimated as

(8.1) $$\begin{eqnarray}x_{p}={\displaystyle \frac{K_{Ic}^{2}}{T^{2}}},\end{eqnarray}$$

where $K_{Ic}$ is the fracture toughness and $T$ is the tensile strength of the formation. For instance, for fracture toughness $K_{Ic}=1$  MPa m $^{1/2}$ and tensile strength $T=20$  MPa, the size of the process zone becomes of the order of a few millimetres. Typical tip element size, however, is of the order of one metre or a fraction of a metre. Therefore, the use of LEFM is justified for these parameters.

8.2 Fluid lag

The model assumes no fluid lag, which is a cavity between the fluid front and the crack front. To estimate the length of the lag size, we need to have $p=-\unicode[STIX]{x1D70E}_{0}$ , where $\unicode[STIX]{x1D70E}_{0}$ is the compressive stress at the location of the fracture (Detournay & Garagash Reference Detournay and Garagash2003). The net pressure can be estimated by integrating the lubrication equation (2.7), so that the fluid lag $\unicode[STIX]{x1D706}$ can be estimated form the equation

(8.2) $$\begin{eqnarray}p(\unicode[STIX]{x1D706})=-\int _{\unicode[STIX]{x1D706}}^{\infty }(Vw+2C^{\prime }V^{1/2}x^{1/2})^{n}\frac{M^{\prime }}{w^{2n+1}}\,\text{d}x=-\unicode[STIX]{x1D70E}_{0},\end{eqnarray}$$

where $p(\infty )=0$ is used in calculations, which is valid for $n>0$ . Equation (8.2) can be rewritten in the dimensionless form using (3.3) as

(8.3) $$\begin{eqnarray}\int _{\tilde{\unicode[STIX]{x1D706}}}^{\infty }\left(1+{\displaystyle \frac{\unicode[STIX]{x1D712}}{\tilde{w}}}\right)^{n}\frac{1}{\tilde{w}^{n+1}\tilde{x}^{n}}\,\text{d}\tilde{x}=\tilde{\unicode[STIX]{x1D70E}}_{0},\end{eqnarray}$$

where the normalized lag and stress are given respectively by

(8.4a,b ) $$\begin{eqnarray}\tilde{\unicode[STIX]{x1D706}}=\left({\displaystyle \frac{\unicode[STIX]{x1D706}}{l_{km}}}\right)^{1/2},\quad \tilde{\unicode[STIX]{x1D70E}}_{0}={\displaystyle \frac{\unicode[STIX]{x1D70E}_{0}}{2}}\left({\displaystyle \frac{{K^{\prime }}^{2n}}{M^{\prime }V^{n}{E^{\prime }}^{n+1}}}\right)^{1/(2-n)}={\displaystyle \frac{\unicode[STIX]{x1D70E}_{0}l_{km}^{1/2}}{2K^{\prime }}}.\end{eqnarray}$$

Vertex solutions can be substituted into (8.3) to estimate the size of the fluid lag. The toughness limiting solution should not be influenced by the presence of a fluid lag, therefore we need to consider only the viscous and leak-off vertex solutions. By substituting the viscous and the leak-off limiting solutions into (8.3), the lag size can be estimated as

(8.5a,b ) $$\begin{eqnarray}\tilde{\unicode[STIX]{x1D706}}_{m}=\left[\frac{2+n}{2n\unicode[STIX]{x1D6FD}_{m}^{n+1}\tilde{\unicode[STIX]{x1D70E}}_{0}}\right]^{(2+n)/2n},\quad \tilde{\unicode[STIX]{x1D706}}_{\tilde{m}}=\left[\frac{2(1+n)\unicode[STIX]{x1D712}^{n/(2(1+n))}}{3n\unicode[STIX]{x1D6FD}_{\tilde{m}}^{2n+1}\tilde{\unicode[STIX]{x1D70E}}_{0}}\right]^{(2(1+n))/3n},\end{eqnarray}$$

where $\unicode[STIX]{x1D712}=0$ and $\unicode[STIX]{x1D712}\gg 1$ are used to simplify the integrand respectively for the viscous and leak-off solutions.

It is instructive to compare the estimates (8.5) with calculations in Garagash & Detournay (Reference Garagash and Detournay2000), where similar problem for Newtonian fluid ( $n=1$ ) and no leak-off ( $\unicode[STIX]{x1D712}=0$ ) is considered together with fluid lag. Solution in Garagash & Detournay (Reference Garagash and Detournay2000) is given in terms of the dimensionless lag $\unicode[STIX]{x1D6EC}=8\tilde{\unicode[STIX]{x1D706}}^{2}\tilde{\unicode[STIX]{x1D70E}}_{0}^{3}$ and the dimensionless toughness $\unicode[STIX]{x1D705}=(2\tilde{\unicode[STIX]{x1D70E}}_{0})^{1/2}$ . Solution (8.5) corresponds to the case $\unicode[STIX]{x1D705}=0$ in Garagash & Detournay (Reference Garagash and Detournay2000). For the latter case, the dimensionless lag is calculated as $\unicode[STIX]{x1D6EC}=0.36$ and this is the largest value of the normalized lag since $\unicode[STIX]{x1D6EC}$ decreases with $\unicode[STIX]{x1D705}$ . The first equation in (8.5), on the other hand, predicts much smaller fluid lag of $\unicode[STIX]{x1D6EC}=0.028$ , but the scaling is captured precisely. This demonstrates that the estimates (8.5) can substantially underestimate the value of the fluid lag, while, at the same time, they properly capture scaling with respect to problem parameters. As a result, they still can be used, but at least $O(10)$ safety factor needs to be applied.

To estimate the size of the fluid lag, let us consider $E=20$  GPa, $\unicode[STIX]{x1D708}=0.2$ , $K_{Ic}=1$  MPa m $^{1/2}$ , $K=0.1$  Pa s, $n=1$ , $V=1$  m s $^{-1}$ , $C^{\prime }=10^{-5}~\text{m}~{\text{s}^{-1}}^{1/2}$ and $\unicode[STIX]{x1D70E}_{0}=20$  MPa. The estimates in (8.5) predict the fluid lag of almost 2 mm and under 1 mm respectively. Given the ‘safety factor’ of the order of 10, the expected lag size is a few centimetres long, which is still much smaller than the typical element size of one metre. The fluid lag is very sensitive to the value of $n$ and decreases by nearly two orders of magnitude if $n$ is reduced to 0.7 (all other parameters are kept the same). However, for this result to be valid, one should also check validity of the power-law model for approximating the fluid rheology at the shear rates that occur in the vicinity of the fluid lag.

8.3 Time dependence

One of the assumptions of the model is that the fracture propagates steadily with the velocity $V$ . In practical applications, however, the fracture front velocity evolves in time, i.e. $V=V(t)$ . To estimate velocity rates, at which the model is still applicable, note that the time derivative term $\unicode[STIX]{x2202}w/\unicode[STIX]{x2202}t$ is neglected in (2.5). If velocity changes slowly, then this term should be small, i.e.

(8.6) $$\begin{eqnarray}{\displaystyle \frac{\unicode[STIX]{x2202}w}{\unicode[STIX]{x2202}t}}\ll V{\displaystyle \frac{\unicode[STIX]{x2202}w}{\unicode[STIX]{x2202}x}}.\end{eqnarray}$$

The velocity rate can be estimated for each regime of propagation by substituting the vertex solutions into (8.6). The toughness solution does not depend on the velocity and hence automatically satisfies (8.6) for any time dependence of the velocity. Note that quasi-static elasticity is used, in which case dynamic effects are still considered negligible. Substitution of the viscous solution into (8.6) gives

(8.7) $$\begin{eqnarray}{\displaystyle \frac{\unicode[STIX]{x2202}w_{m}}{\unicode[STIX]{x2202}V}}{\displaystyle \frac{\text{d}V}{\text{d}t}}={\displaystyle \frac{n}{2+n}}{\displaystyle \frac{w_{m}}{V}}{\displaystyle \frac{\text{d}V}{\text{d}t}}\ll V{\displaystyle \frac{\unicode[STIX]{x2202}w_{m}}{\unicode[STIX]{x2202}x}}={\displaystyle \frac{2}{2+n}}{\displaystyle \frac{Vw_{m}}{x}}\quad \Rightarrow \quad {\displaystyle \frac{\text{d}V}{\text{d}t}}\ll {\displaystyle \frac{2}{n}}{\displaystyle \frac{V^{2}}{x}}.\end{eqnarray}$$

Substitution of the leak-off solution into (8.6) gives

(8.8) $$\begin{eqnarray}{\displaystyle \frac{\unicode[STIX]{x2202}w_{\tilde{m}}}{\unicode[STIX]{x2202}V}}{\displaystyle \frac{\text{d}V}{\text{d}t}}={\displaystyle \frac{n}{4(1+n)}}{\displaystyle \frac{w_{\tilde{m}}}{V}}{\displaystyle \frac{\text{d}V}{\text{d}t}}\ll V{\displaystyle \frac{\unicode[STIX]{x2202}w_{\tilde{m}}}{\unicode[STIX]{x2202}x}}={\displaystyle \frac{4+n}{2(1+n)}}{\displaystyle \frac{Vw_{\tilde{m}}}{x}}\quad \Rightarrow \quad {\displaystyle \frac{\text{d}V}{\text{d}t}}\ll {\displaystyle \frac{2(4+n)}{n}}{\displaystyle \frac{V^{2}}{x}}.\end{eqnarray}$$

Results in (8.7) and (8.8) demonstrate that the fracture velocity difference that is calculated over time it takes the fracture to propagate the distance $x$ should be much smaller than the velocity itself. There is also an explicit dependence on $n$ , which shows that more rapid changes in the velocity are allowed for smaller $n$ . The time dependence can be completely ignored for perfectly plastic fluids with $n=0$ .

If one considers a radial fracture, then the radius varies with time as $R\propto t^{\unicode[STIX]{x1D6FC}}$ , where $\unicode[STIX]{x1D6FC}=O(1)$ changes based on the regime of propagation, e.g. Madyarova (Reference Madyarova2003), Dontsov (Reference Dontsov2016a ). If we then take $n=O(1)$ and neglect all coefficients that are $O(1)$ , then (8.7) and (8.8) translate into

(8.9) $$\begin{eqnarray}{\displaystyle \frac{R}{t^{2}}}\ll {\displaystyle \frac{R^{2}}{xt^{2}}}\Rightarrow {\displaystyle \frac{x}{R}}\ll 1.\end{eqnarray}$$

Here $x$ is the size of the tip element, for which the asymptotic solution would be used. Equation (8.9) demonstrates that the condition (8.6) is automatically satisfied for radial fractures since many elements are typically used per fracture radius. Similar considerations can be applied to fractures with different geometry, such as plane strain fractures.

If a fracture encounters a layer boundary, then it can accelerate or decelerate rapidly. However, the current model does not capture fracture shape change due to the presence of spatial variation of properties and therefore is not directly applicable within a vicinity of a layer boundary.

8.4 Fluid leak-off

Carter’s model is used in this study and it assumes one-dimensional leak-off normal to the fracture surface. This is an assumption that is generally valid for long fractures, for which the diffusion length scale is much smaller than the fracture size. In terms of semi-infinite hydraulic fractures, the work (Kovalyshen & Detournay Reference Kovalyshen and Detournay2013) studies the semi-infinite hydraulic fracture in a poroelastic medium and considers fully three-dimensional leak-off. As indicated there, one should consider a competition between the crack front and the diffusion front. The characteristic parameter for such a competition is the ratio between the diffusion length scale and the fracture size, i.e.

(8.10) $$\begin{eqnarray}{\displaystyle \frac{\sqrt{ct}}{x}}=\sqrt{{\displaystyle \frac{c}{Vx}}}\ll 1,\quad \Rightarrow \quad x_{d}={\displaystyle \frac{c}{V}},\quad x\gg x_{d},\end{eqnarray}$$

where $c$ is the diffusion coefficient for the poroelastic problem. Similar to the plastic zone size, there is a diffusion zone size $x_{d}$ , which should be much smaller than the distance to the fracture tip for the model to apply. The value of the diffusion coefficient may vary significantly from one formation to another. By considering the case $c=O(10^{-3})~\text{m}^{2}~\text{s}^{-1}$ and $V=0.1~\text{m}~\text{s}^{-1}$ , the diffusion length scale becomes of the order of a centimetre and therefore Carter’s model can be used for such parameters provided that the typical element size is one metre.

9 Summary

This paper investigates the problem of a semi-infinite hydraulic fracture that is driven by a power-law fluid in a permeable medium. In particular, the goal of this paper is to develop a numerical solution for the problem and to construct an accurate approximation that allows us to evaluate the solution rapidly. To develop an efficient numerical solution, the governing equations are first reduced to a single integral equation with a non-singular kernel, which greatly simplifies the original problem that involves a singularity. The solution features three limiting regimes of propagation that are associated with dominance of either toughness, viscosity or leak-off. Numerical solutions for different values of the fluid rheology exponent $n$ are obtained and applicability regions of the limiting solutions are analysed. In addition to calculating the bounds of applicability numerically, they are also calculated by estimating the contribution of the higher-order term relative to the leading behaviour. One peculiar feature of the problem is that the parametric space degenerates for perfectly plastic fluids $n=0$ since the viscosity and leak-off dominated limiting solutions coincide for this case and there is no influence of leak-off on the solution. An additional outcome of this paper is the development of an accurate approximation that allows us to rapidly evaluate the solution inside the whole parametric space. This solution is thoroughly tested against the numerical solution, which showed that the maximum error of approximately 1 % occurs for small values of the power-law exponent $n$ , i.e. $0.1\lesssim n\lesssim 0.25$ . This error drops rapidly for larger values of $n$ . Since the problem of a semi-infinite hydraulic fracture describes the tip region of a hydraulic fracture, the developed approximate solution can be used for rapid evaluation of the near-tip solution, which in turn can be used as a propagation condition in hydraulic fracturing simulators.