2.1. Governing equations

The model problem is sketched in figure 1, with a Newtonian liquid being injected into an elastic solid confined between two parallel plates of gap thickness $b$. Neglecting the influence of body forces and interfacial tension, standard lubrication theory of viscous flow in a Hele-Shaw cell then indicates that the profile shape $h(x,t)$ of a fracture is governed by a partial differential equation

(2.1)\begin{equation} \frac{\partial h}{\partial t} - \frac{b^2}{12\mu}\,\frac{\partial}{\partial x} \left(h\,\frac{\partial p}{\partial x}\right) = 0, \end{equation}

with $p(x,t)$ denoting the distribution of pressure within the liquid film that is to balance the normal stress of the elastically deformed solid at the fluid–solid interface. Equation (2.1) is different from (2.25) of Spence & Sharp (Reference Spence and Sharp1985), since the viscous drag within the liquid film is due predominantly to the no-slip boundary condition at the inner surface of the parallel plates rather than at the fluid–solid interface.

We constrain the model problem in the positive half-domain $x \in [0,\infty ]$. The pressure distribution $p(x,t)$ within the liquid film due to elastic deformation is then given by

(2.2) \begin{equation} p(x,t) ={-}\frac{E}{2(1-\sigma^2)}\,\frac{1}{\rm \pi} \int_0^{x_f(t)} \frac{1}{s-x}\,\frac{\partial h}{\partial s} \,{\rm d}s, \end{equation}

where $E$ is the Young's modulus of the elastic material, $\sigma$ is the Poisson ratio, and $x_f(t)$ is the location of the propagating front. The pressure distribution (2.2) due to elastic deformation is consistent with those in previous studies of plane and radial fractures in infinite elastic solids (e.g. Spence & Sharp Reference Spence and Sharp1985; Lister Reference Lister1990a,Reference Listerb). Since we focus only on the positive half-domain $x \in [0,\infty ]$, to be consistent with the experiment of e.g. Weitz (Reference Weitz2020), the lower limit of integration is $x = 0$ in (2.2) rather than $x = -x_f(t)$, which is the only difference from the well-cited literature (e.g. Spence & Sharp Reference Spence and Sharp1985; Lister Reference Lister1990a,Reference Listerb). The prefactor $1/{\rm \pi}$ is introduced from the definition of the Hilbert transform.

The differential–integral system (2.1) and (2.2) is to be solved providing appropriate initial and boundary conditions (IBCs), e.g. given by

(2.3a)$$\begin{gather} h(x,0) = 0, \end{gather}$$

(2.3b)$$\begin{gather}h(x_f(t),t) = 0, \end{gather}$$

(2.3c)$$\begin{gather}\int_0^{x_f(t)}h(x,t)\,{\rm d}\kern0.7pt x = qt. \end{gather}$$

The initial condition (2.3a) indicates that there exists no fluid initially within the elastic solid, and fluid injection (i.e. fracturing) starts instantaneously at $t=0$. The first boundary condition (2.3b) is a standard frontal condition (with compact support). The second boundary condition (2.3c) is a global constraint for the increase of total area within the plane of the paper, as covered by the fracturing fluid due to injection that proceeds at a constant rate $q$.

The global volume conservation condition (2.3c) can also be shown to be equivalent to a flux condition at the injection point ($x=0$):

(2.4)\begin{equation} \left.- \frac{b^2}{12\mu}\,h\,\frac{\partial p}{\partial x} \right|_{x=0} = q. \end{equation}

The treatment is similar to that for the propagation of a gravity current (e.g. Zheng et al. Reference Zheng, Guo, Christov, Celia and Stone2015). One can integrate the evolution equation (2.1) from $x=0$ towards $x=x_f(t)$, imposing a zero-flux condition at the location of the propagating front $x=x_f(t)$ based on the assumption that there is no fluid loss or entrainment locally at $x=x_f(t)$. The integral constraint (2.3c) then reduces to a nonlinear flux condition (2.4). It is more convenient to solve the coupled evolution equations (2.1) and (2.2) numerically, subject to initial condition (2.3a), frontal condition (2.3b) and flux condition (2.4), e.g. with finite-volume schemes.

2.2. Scaling analysis

The form of the governing equations (2.1), (2.2) and IBCs (2.3) indicates the following scaling behaviours for the evolution of the length scale $x$, thickness scale $h$, and pressure scale $p$:

(2.5a) $$\begin{gather} x \propto t^{1/2} \left[\frac{b^2 E q}{24 {\rm \pi}(1-\sigma^2) \mu}\right]^{1/4}, \end{gather}$$

(2.5b)$$\begin{gather}h \propto t^{1/2} \left[\frac{24 {\rm \pi}(1-\sigma^2) \mu q^3}{b^2 E}\right]^{1/4}, \end{gather}$$

(2.5c)$$\begin{gather}p \propto \left[\frac{6 E \mu q}{{\rm \pi} (1-\sigma^2) b^2}\right]^{1/2}. \end{gather}$$

These scaling behaviours apply at intermediate times when $x \gg b$ and $h \gg b$, such that the resistance for the elongation of a fracture comes predominantly from the viscous drag due to the no-slip condition at the surface of the parallel plates rather than at the fluid–solid interface. An estimate for the appropriate time scale is immediately available as $t \gg \max ([24{\rm \pi} (1-\sigma ^2)b^2 \mu /Eq]^{1/2}, [Eb^6/(24{\rm \pi} (1-\sigma ^2)\mu q^3)]^{1/2})$.

It is also seen that these scaling laws (2.5a–c) during hydraulic fracturing in a Hele-Shaw cell (e.g. $x \propto t^{1/2}$ and $h \propto t^{1/2}$) are different from those (e.g. $x \propto t^{2/3}$ and $h \propto t^{1/3}$) for the growth of two-dimensional plane fractures of infinite depth (e.g. Spence & Sharp Reference Spence and Sharp1985). This is due to the difference of the velocity field of the lubricating flow and the viscous drag, which leads to a different scaling exponent of $h$ in the thin film equation (2.1).

2.3. Non-dimensionalisation

The scaling results also suggest the existence of self-similar solutions for the evolution of the fracture shape $h(x,t)$ and pressure distribution $p(x,t)$. It is standard first to rescale the governing system (2.1) and (2.2) and the IBCs (2.3a–c). Since there exist no natural time or length scales within the plane of the fracture, we chose the gap thickness $b$ as a reference length scale, and define $h_c = b$. This would correspond to a time scale $t_c$ when the thickness of a fracture reaches $b$. Based on this choice, our model would work in the late-time regime $t \gg t_c$, when the length and thickness scales are both much greater than $b$. Accordingly, we define dimensionless variables as

(2.6a–d)\begin{equation} \bar{t} \equiv \frac{t}{t_c}, \quad \bar{x} \equiv \frac{x}{x_c}, \quad \bar{h} \equiv \frac{h}{b}, \quad \bar{p}\equiv \frac{p}{p_c}, \end{equation}

where the characteristic time, length and pressure scales are chosen as

(2.7a–c) \begin{align} t_c = \left[\frac{b^6 E}{24 {\rm \pi} (1-\sigma^2) \mu q^3}\right]^{1/2}, \quad x_c = \left[\frac{b^4 E}{24 {\rm \pi}(1-\sigma^2) \mu q}\right]^{1/2},\quad p_c = \left[\frac{6 \mu q E }{{\rm \pi} (1-\sigma^2) b^2}\right]^{1/2}. \end{align}

We then arrive at the dimensionless form of (2.1), (2.2) and IBCs (2.3), and for simplicity we continue to use $t$, $x$, $h$ and $p$ to represent dimensionless variables in this section and § 3 from now on:

(2.8a)$$\begin{gather} \frac{\partial h}{\partial t} = \frac{\partial}{\partial x} \left(h\,\frac{\partial p}{\partial x}\right), \end{gather}$$

(2.8b)$$\begin{gather}p(x, t) ={-}\int_0^{x_f(t)} \frac{1}{s-x}\,\frac{\partial h}{\partial s} \,{\rm d}s, \end{gather}$$

and

(2.9a)$$\begin{gather} h(x,0) = 0, \end{gather}$$

(2.9b)$$\begin{gather}h(x_f(t), t) = 0, \end{gather}$$

(2.9c)$$\begin{gather}\int_0^{x_f(t)} h(x, t)\,{\rm d}\kern0.7pt x = t. \end{gather}$$

Meanwhile, the dimensionless flux condition (2.4) is given by

(2.10)\begin{equation} \left.- h\,\frac{\partial p}{\partial x} \right|_{x=0} = 1, \end{equation}

which is equivalent to the integral constraint (2.9c). A finite-volume scheme is described in Appendix A that is employed to solve numerically the dimensionless system (2.8a,b) subject to (2.9a–c) for the coupled evolution of the dimensionless profile shape $h(x, t)$ and pressure distribution $p(x, t)$.

Representative numerical solutions are shown in figure 2 for the time evolution of the profile shape $h(x, t)$ and pressure distribution $p(x, t)$ for a hydraulic fracture within a Hele-Shaw cell. In particular, it is observed that the fracture develops into the cusp shape in a Hele-Shaw cell, which is consistent with the experimental observation in figure 1(a). This is also completely different from that of plane fractures of infinite depth with profiles in the form $h \propto (x_f - x)^{2/3}$ as $x \to x_f^-$. (In the viscous regime, see e.g. Spence & Sharp Reference Spence and Sharp1985.) The inset of figure 2(a) also shows that the normalised profile solutions collapse onto a universal shape, which suggests the existence of self-similarity, as we discuss later. Meanwhile, numerical solutions for the pressure distribution indicate that the pressure approaches a finite value as $x \to x_f^-$, which is also fundamentally different from the well-known pressure singularity of plane and radial fractures in the Cartesian and radial configurations (e.g. Spence & Sharp Reference Spence and Sharp1985; Spence et al. Reference Spence, Sharp and Turcotte1987; Lister Reference Lister1990a,Reference Listerb). Later, in § 3.1, these numerical solutions are also rescaled appropriately and compared with the self-similar solutions that we explore next.

Figure 2.

(a) Time evolution of cusp-like profile shapes $h(x,t)$ near the front, and (b) time evolution of the pressure field $p(x, t)$, based on solving numerically the full partial differential–integral system (2.8a,b). The inset of (a) shows also that the normalised profiles collapse onto a universal shape. The profiles were taken at time $t = 250 \times \{1\unicode{x2013}9\}$. The domain length is $L=600$, with 1200 space grids employed here.

2.4. Self-similar solutions

To look for self-similar solutions of (2.8a,b) and (2.9a–c) at intermediate times ($t \gg 1$), we start by defining a similarity variable as

(2.11)\begin{equation} \xi \equiv x/t^{1/2}. \end{equation}

We then look for $h(x, t) = t^{1/2}\,k(\xi )$ and $p(x, t) = g(\xi )$ by solving

(2.12a)$$\begin{gather} \frac{1}{2}\,k-\frac{1}{2}\,\xi\,\frac{{\rm d} k}{{\rm d} \xi} = \frac{{\rm d}}{{\rm d} \xi} \left(k\, \frac{{\rm d} g}{{\rm d} \xi}\right), \end{gather}$$

(2.12b)$$\begin{gather}g(\xi) ={-} \int_0^{\xi_f} \frac{1}{s-\xi}\,\frac{{\rm d}k}{{\rm d} s} \,{\rm d}s, \end{gather}$$

subject to

(2.13{a,b})\begin{equation} k(\xi_f) = 0 \quad \mbox{and}\quad \int_0^{\xi_f} k(\xi)\,{\rm d}\xi = 1. \end{equation}

By definition, $\xi _f \equiv x_f/t^{1/2}$ represents the frontal location. Meanwhile, the influence of initial condition (2.9a) no longer exists, as we are looking for intermediate time behaviours when $t \gg 1$.

To solve numerically the ordinary differential–integral system (2.12a,b) for the self-similar profile shape $k(\xi )$ and pressure distribution $g(\xi )$, we keep the time dependence by introducing a mathematical transform

(2.14a,b)\begin{equation} \tau = \ln t\quad \mbox{and}\quad \xi = x/t^{1/2}. \end{equation}

The original partial differential-integral system (2.8a,b) then becomes

(2.15a)$$\begin{gather} -\frac{\partial k}{\partial \tau} + \frac{1}{2}\,k - \frac{1}{2}\,\xi\,\frac{\partial k}{\partial \xi} = \frac{\partial}{\partial \xi} \left(k\,\frac{\partial g}{\partial \xi}\right), \end{gather}$$

(2.15b)$$\begin{gather}g ={-} \int_0^{\xi_f} \frac{1}{s-\xi}\,\frac{\partial k}{\partial s} \,{\rm d}s. \end{gather}$$

Compared with (2.12a,b), (2.15a,b) now include an additional term of time evolution $\partial k/\partial \tau$. We then solve this evolution system for $k(\xi,\tau )$ and $g(\xi,\tau )$ using the same scheme as described in Appendix A, now subject to appropriate IBCs:

(2.16a)$$\begin{gather} k(\xi,0) = 0, \end{gather}$$

(2.16b)$$\begin{gather}k(\xi_f, \tau) = 0\quad \mbox{and}\quad \left. k\,\frac{\partial g}{\partial \xi} \right|_{\xi = 0} ={-}1. \end{gather}$$

This system is solved until there is negligible time evolution, as shown in figure 3, due to the existence of a sink term. The steady-state solutions $k(\xi,\infty ) \to k_s(\xi )$ and $g(\xi,\infty ) \to g_s(\xi )$ become, effectively, the self-similar solutions of the ordinary differential–integral system (2.12a,b) that we are looking for. From figure 3, we also obtain the end-point values as $\xi _f \approx 2.42$, $k_s(0) \approx 0.82$, $g_s(0) \approx 1.37$ and $g_s(\xi _f) \approx -1.35$. For $g(\xi )$, relatively more significant numerical error appears around $\xi = \xi _f$, but this does not seem to influence the bulk structure. The stability of the self-similar solutions obtained here can also be studied based on the time-dependent system (2.15a,b) by adding small perturbations. Similar ideas have been employed before to study the development of finite-time singularities for thin fluid threads (e.g. Eggers & Fontelos Reference Eggers and Fontelos2009). To some extent, the transient dynamics here is similar to that of a propagating viscous gravity current that suffers also slow drainage from thin permeable substrates, driven by buoyancy or background flow (e.g. Pritchard, Woods & Hogg Reference Pritchard, Woods and Hogg2001). With fluid supply at a constant rate, the model problem evolves into steady flow solutions after an initial transition period, when fluid injection balances drainage.

Figure 3.

The partial differential–integral system (2.15a,b) is solved until there is negligible time evolution, i.e. $k(\xi,\infty ) \to k_s(\xi )$, and $g(\xi,\infty ) \to g_s(\xi )$. The steady-state solutions are, effectively, the self-similar solutions of (2.12a,b). Here, the domain length is $L=3$ with 300 space grids employed, and we have included solutions at $\tau = 0.8\times \{1\unicode{x2013}8\}$.

It is of particular interest to note that the normalised solutions in figure 3 suggest that pressure decreases linearly towards the fracture's front while the profile shape develops into a cusp shape. It can also be shown that the following analytical solutions provide very good predictions:

(2.17a)$$\begin{gather} h/h_0 ={-}{\rm \pi}^{{-}1}\left[(2\eta - 1)[1-(2\eta-1)^2]^{1/2} + \sin^{{-}1}(2\eta-1) - {\rm \pi}/2\right], \end{gather}$$

(2.17b)$$\begin{gather}p/p_0 = (1 - 2\eta), \end{gather}$$

where $\eta \equiv \xi /\xi _f = x/x_f \in [0,1]$, as we included also in figure 3 as well. In particular, (2.17a) indicates that

(2.18)\begin{equation} h/h_0 \propto (1-x/x_f)^{3/2}, \end{equation}

as $x \to x_f(t)^-$, which describes the cusp shape. More discussions on the frontal structure (2.17a,b) and the universal behaviours of the profile shape and pressure distribution are provided in Appendix B based on series expansions.

3.1. Comparison with time-dependent numerical solutions

We first provide a comparison between the self-similar solution and the time-dependent numerical solutions of the partial differential–integral system (2.8a,b). In particular, the time-dependent numerical solutions for the evolution of the profile shape and pressure distribution in figure 2 are now rescaled appropriately based on the maximum values in figure 4, which exhibits very good collapse onto universal curves. We then include in the same figure the (appropriately stretched) self-similar solutions from solving numerically (2.15a,b), which provides very good agreement with the collapsed numerical solutions.

Figure 4.

The rescaled solutions of the profile shape and pressure field evolution of the numerical data in figure 2. They already collapse onto universal curves, which also exhibits reasonably good agreement with the prediction of the self-similar solutions from solving (2.15a,b) until there is negligible time evolution, shown as the dashed curves.

Meanwhile, we track the time-dependent frontal location $x_f(t)$, inlet thickness $h_0(t) \equiv h(0, t)$, inlet pressure $p_0(t) = p(0, t)$ and tip pressure $p_f(t) \equiv p(x_f(t), t)$ from numerical solutions of (2.8a,b). In all cases, the numerical solutions (the symbols) approach the prediction of the self-similar solutions from solving numerically (2.15a,b) at intermediate times, shown as the solid lines in figure 5, and as described by

(3.1a–d)\begin{equation} x_f(t) \sim 2.42\, t^{1/2},\quad h_0(t) \sim 0.82\, t^{1/2},\quad p_0(t) \sim 1.37\quad \mbox{and}\quad p_f(t) \sim{-}1.35. \end{equation}

There is very good agreement for both the scaling exponents of time dependence and the prefactors between the self-similar and numerical solutions. The very small difference, we believe, is more likely due to numerical error (e.g. on the last digit of solution (3.1a–d)) and is less likely due to unfinished time transition.

Figure 5.

Time evolution of the frontal location $x_f(t)$, fracture thickness $h_0(t) \equiv h(0,t)$, inlet pressure $p_0(t) \equiv p(0, t)$, and tip pressure $p_f(t) \equiv p(x_f(t), t)$, based on solving numerically the partial differential–integral system (2.8a,b). The numerical results (symbols) agree very well with the prediction of the self-similar solutions (solid lines) at intermediate times.

3.2. Comparison with available experimental observations

There are also limited experimental pictures available from the literature (e.g. Weitz Reference Weitz2020; Aime et al. Reference Aime, Sabato, Xiao and Weitz2021). In a typical experiment, a fluid such as air or water is injected through a syringe pump into a thin layer of colloidal gel that is confined between the two parallel plates of a Hele-Shaw cell, as shown in figure 6. The porosity of the gel (1.875 % Ludox Gel, SM 40) is $\approx$70 %, the modulus is ${\approx }100$ kPa, and the yield stress is $\approx$480 Pa from measurements. The Hele-Shaw cell has length 7.5 cm, width 5.0 cm, and thickness $7\ \mathrm {\mu }{\rm m}$. The injection rate is small, and typical propagating speed of a hydraulic fracture is ${\approx }1\ \mathrm {\mu }{\rm m}\ {\rm s}^{-1}$. A typical experiment lasts for 1–10 h.

Figure 6.

Experimental pictures at three different times for the profile shape of a hydraulic fracture in a Hele-Shaw cell (e.g. Weitz Reference Weitz2020), while the appropriately stretched theoretical predictions from solving (2.8a,b) are superimposed as the cyan curves.

A comparison between experimentally recorded (e.g. Weitz Reference Weitz2020) and rescaled (and collapsed) numerical solutions of the full partial differential–integral system (2.8a,b) is shown in figure 6 for the profile shape of a hydraulic fracture in a Hele-Shaw cell. Experimental pictures of the near-tip region of the profile shape at three different times are shown in figures 6(a–c), respectively, with the appropriately stretched numerical solutions superimposed on each of them. In particular, the length scales at two different times $t_1$ and $t_2$ must satisfy $x_1/x_2 = h_1/h_2 = (t_1/t_2)^{1/2}$ at the self-similar stage based on the scaling laws (2.5a,b). We observe fairly good agreement between the theoretical and experimental results of this cusp-shaped hydraulic fracture in a Hele-Shaw cell. We do recognise that the available experimental data are rather limited. So we suggest further experiments being performed in the future with time-dependent profile shapes and possibly pressure fields reported, which would allow more thorough comparison with the theory.

4.1. Summary

A theoretical model is provided to describe the coupled evolution of the profile shape and pressure distribution of a hydraulic fracture in a Hele-Shaw cell. The model is based on the lubrication theory of viscous thin film flow in an elastically deformable cavity, which is coupled with the linear elastic theory of solid deformation for the distribution of normal stresses at the fluid–solid interface (i.e. the distribution of pressure within the liquid film). It is of particular interest to note that the model leads to hydraulic fractures of cusp shapes, which is consistent with recent experimental observations in a Hele-Shaw cell. The model also suggests the existence of self-similar solutions at intermediate times, which is verified through a comparison with rescaled numerical solutions of the full partial differential–integral system (2.8a,b). Compared with the classic results of plane or penny-shaped fractures, there exists no pressure singularity at the tip of a hydraulic fracture in a Hele-Shaw cell, which is a fundamental difference. We suggest that more experiments should be performed in the future to obtain time-dependent data of the fracture shape, frontal location, and pressure distribution (if possible), which can be used to further verify the theory and/or suggest new directions for future explorations.

4.2. Time-dependent injection modes

The theoretical model works also when the injection rate is time-dependent, as already pointed out in earlier studies on the propagation of plane and radial hydraulic fractures in elastic solids (e.g. Spence & Sharp Reference Spence and Sharp1985; Lai et al. Reference Lai, Zheng, Dressaire, Wexler and Stone2015). This is also true for the growth of a hydraulic fracture in a Hele-Shaw cell in the current work. Of particular interest is the situation when the injection rate follows power-law or exponential forms of time dependence, in which case self-similar solutions can also be explored to describe the universal behaviours at intermediate times. Here, we briefly point out that under a power-law form of fluid injection, such that the total volume (area) follows $V(t) = q t^\alpha$, the scaling results (2.5a–c) are extended to

(4.1a) $$\begin{gather} x \propto t^{(\alpha+1)/4} \left[\frac{b^2 Eq}{24 {\rm \pi}(1-\sigma^2) \mu}\right]^{1/4}, \end{gather}$$

(4.1b)$$\begin{gather}h \propto t^{(3\alpha-1)/4} \left[\frac{24 {\rm \pi}(1-\sigma^2) \mu q^3}{b^2 E}\right]^{1/4}, \end{gather}$$

(4.1c)$$\begin{gather}p \propto t^{(\alpha-1)/2} \left[\frac{6 \mu q E }{{\rm \pi} (1-\sigma^2) b^2}\right]^{1/2}, \end{gather}$$

to describe the time evolution of the length and pressure scales. We stop the discussion here by noting simply that self-similar solutions for the profile shape $h(x,t)$ and pressure distribution $p(x,t)$ can also be explored following the procedure in § 2.4 of the current work.

4.3. The influence of material toughness

It is well known that the propagation of a hydraulic fracture can possibly be resisted by material breaking at the fracture tip as well. In the well-cited literature of hydraulic fracturing of an infinite elastic solid that is also brittle (e.g. Spence & Sharp Reference Spence and Sharp1985; Lister Reference Lister1990b; Savitski & Detournay Reference Savitski and Detournay2002; Roper & Lister Reference Roper and Lister2007), a revised model is to introduce a propagation condition

(4.2)\begin{equation} K_{c} \sim \frac{E}{2^{3/2}(1-\sigma^2)}\,\frac{h(x,t)}{(x_f(t)-x)^{1/2}}, \end{equation}

as $x \to x_f(t)^-$. Condition (4.2) is consistent with those derived based on the theory of linear elastic fracture mechanics of a mode I crack without the influence of body forces and that of flow (e.g. Rice Reference Rice1968; Howell, Kozyreff & Ockendon Reference Howell, Kozyreff and Ockendon2008). The critical stress intensity factor $K_{c} > 0$ is typically named ‘material toughness’ and is an important property of a brittle material. Detailed analysis of the influence of material breaking on hydraulic fracturing in a Hele-Shaw cell is left for a future work. In such a context, the description of (2.1a,b) subject to IBCs (2.3) applies, when fracturing a non-brittle elastic material as $K_{c} \to 0^+$.