2.1. Governing equations

During hydraulic fracturing of an infinite elastic solid, as shown in figure 1(a), the coupled evolution of the fracture’s profile shape $h(x,t)$ and the pressure distribution $p(x,t)$ within the viscous fluid film can be described by (e.g. Spence & Sharp Reference Spence and Sharp1985)

(2.1a) \begin{equation} \frac {\partial h}{\partial t} - \frac {1}{3\mu }\frac {\partial }{\partial x}\left (h^3\frac {\partial p}{\partial x}\right ) = 0, \end{equation}

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

under the assumption of lubrication approximation and neglecting the influence of buoyancy. In (2.1a , b ), $\pm x_f(t)$ denotes the location of the propagating fronts, $\mu$ denotes the viscosity of the fracturing fluid, and $E$ and $\sigma$ denote the Young’s modulus and Poisson ratio of the elastic solid, respectively. Physically, (2.1a ) describes the lubricating thin film flow within a hydroelastic cavity, while (2.1b ) describes the distribution of fluid pressure due to elastic deformation of the solid. Equation (2.1a ) can also be modified appropriately to describe the dynamics of hydraulic fracturing driven by turbulent flow (e.g. Lister Reference Lister1990a ). Equation (2.1b ) can be derived based on the theory of linear elasticity and applies also for the region ahead of the fluid fronts $\pm x_f(t)$ along $h(x,t)=0$ (e.g. Spence & Sharp Reference Spence and Sharp1985).

The nonlinear partial differential–integral system (2.1a , b ) is symmetric about $x=0$ and can be solved subject to appropriate initial and boundary conditions (IBCs) such as

(2.2a) \begin{align} h(x,0) &= 0, \end{align}

(2.2b) \begin{align} h(\pm x_f(t),t) &= 0, \end{align}

(2.2c) \begin{align} \int _{-x_f(t)}^{x_f(t)} h(x,t)\, \textit { d}x &= 2qt, \end{align}

with $q$ denoting a constant area injection rate into the first quadrant (the upper right-hand corner of the plane). Physically, (2.2a ) indicates that initially there exists no fluid within the elastic solid, (2 .2b ) indicates the existence of a propagating front with zero thickness and (2.2c ) is a global constraint that represents the conservation of total fluid volume within a fracture, which can be converted into an equivalent flux boundary condition at $x=0$ as we show later. Conditions (2.2b , c ) apply for all time $t$ .

We have neglected the influence of material toughness $K_c$ in the aforementioned problem formulation. It is also seen that a scaling analysis of (2.1a , b ) subject to boundary conditions (2.2b ) and (2.2c ) leads to

(2.3) \begin{equation} x \propto t^{2/3},\,\,\,h \propto t^{1/3},\,\,\, \mbox {and}\,\,\,p \propto t^{-1/3}, \end{equation}

at intermediate times for the growth of a plane fracture, which is consistent with those obtained in Spence & Sharp (Reference Spence and Sharp1985). Self-similar solutions can also be obtained for $h(x,t)$ and $p(x,t)$ by defining a similarity variable $\xi \equiv x/t^{2/3}$ accordingly. The physics is then very much analogous to that of a viscous gravity current propagating over a horizontal substrate (e.g. Huppert Reference Huppert1982). The main difference is that the pressure distribution is now given by an integral expression rather than an algebraic expression, but scaling analysis can still be employed when the size of the elastic solid is infinite.

The physics could be altered, however, if the size of the elastic solid is finite, compared with that of a fracture, since the stress boundary condition imposed in the far-field could impact the overall stress distribution in the elastic solid and hence the flow. We have neglected the influence of any externally imposed far-field stresses in the present work. In the case that the influence of an imposed far-field normal stress is not negligible, it has been shown that a lag region would possibly appear ahead of the fluid front, filled with either the vapour phase of the fracturing fluid or the ambient fluid leaked in from the porous elastic solid. Detailed asymptotic analysis of this lag region under the influence of a far-field normal stress has been performed (e.g. Garagash Reference Garagash2006), with experimental evidence available when fracturing a block of polymethyl methacrylate by injecting glycerine (e.g. Bunger et al. Reference Bunger, Detournay and Jeffrey2005).

Within the framework of Garagash (Reference Garagash2006), where a propagation condition is employed considering the contribution of both the fluid and lag regions, imposing a positive far-field normal stress, the length of the lag region would decrease with time compared with the fracture in the viscous regime (when material toughness $K_c = 0$ ), and asymptotically a zero-lag solution with a pressure singularity would be obtained. That said, all pressure here can be considered as net pressure, i.e. total pressure minus far-field compressive stress, including the extraction pressure. One can also refer to the analogous problem of elastic Hele-Shaw flows for similar physics with regard to the regularisation of the pressure singularity (e.g. Flitton & King Reference Flitton and King2004; Lister et al. Reference Lister, Peng and Neufeld2013; Hewitt et al. Reference Hewitt, Balmforth and De Bruyn2015). In the present work, we allow the existence of such a pressure singularity at the fluid front without a lag region, similar to the scope of Spence & Sharp (Reference Spence and Sharp1985), Lister (Reference Lister1990a ,Reference Listerb) and Savitski & Detournay (Reference Savitski and Detournay2002), with experimental verifications from Lai et al. (Reference Lai, Zheng, Dressaire, Wexler and Stone2015) and OKeeffe et al. (Reference OKeeffe, Huppert and Linden2018) for leading-order physics of the fracturing process.

The primary focus of the current problem is to study the coupled evolution of the profile shape and the pressure distribution during fluid removal through a narrow conduit that locates at the centre ( $x=0$ ) of a fracture, as shown in figure 1. We assume that fluid removal occurs through a narrow conduit of width $d$ and length $L$ . We also assume that fluid flow is viscous within the conduit with a parabolic velocity distribution, driven by the pressure difference between the elastic pressure $p_0(t) \equiv p(0,t)$ at the inlet of the conduit and the ambient pressure $p_e$ at the outlet. Then, based on Poiseuille’s law, we obtain an area rate of fluid removal,

(2.4) \begin{equation} q_0(t) = \frac {d^3}{12 \mu } \frac {p_0(t)-p_e}{L}, \end{equation}

as $x \to 0^+$ , when the width of the conduit is negligibly small compared with the length of the fracture, i.e. when $d \ll x_f(t)$ . This is the regime we consider in the current work. For the coupled evolution of $h(x,t)$ and $p(x,t)$ during fluid removal, (2.1a , b ) can then be solved subject to an appropriate initial condition, provided at the end of the fracturing process, together with a flux boundary condition (2.4) at $x = 0$ and a frontal condition (2.2b ) at $x = x_f(t)$ .

We consider only the situation of $p_e \leqslant 0$ in this work through the flux condition (2.4), and the value of $|p_e|$ indicates how strong is the influence of the ambient-pressure-driven backflow through extraction (i.e. suction). In particular, $p_e = 0$ represents the special situation of elasticity-driven fluid removal from a hydraulic fracture without extraction. Later we show that the thickness of a fracture eventually approaches zero at the centre ( $x=0$ ), and the fracture evolves into a self-similar shape of the dipole type at intermediate time. Finally, we note also that $p_e \gt p_0(t)$ corresponds to hydraulic fracturing under pressure control, which could be of interest for a future work.

2.2. Non-dimensionalisation

It is standard first to non-dimensionalise the model problem (2.1a , b ) for the coupled evolution of thin film flow and elastic boundary deformation before providing time-dependent numerical solutions and similarity solutions. Since there exist no natural time, length or pressure scales for this moving boundary problem, we choose $t_c = t_i$ based on the time period $t_i$ of initial fracturing, which corresponds also to a target fracture length or thickness. In this case, the initial fracturing process proceeds on $\bar t \in [0,1]$ . The characteristic length, thickness and pressure scales are then chosen based on the model problem of hydraulic fracturing under fluid injection as

(2.5) \begin{align} x_c = t_i^{2/3} \left [\frac {Eq^3}{6\pi \mu (1-\sigma ^2)}\right ]^{1/6},\,\,\, h_c = t_i^{1/3}\left [\frac {6\pi \mu q^3(1-\sigma ^2)}{E}\right ]^{1/6},\nonumber \\ p_c = t_i^{-1/3}\left [\frac {3\mu E^2}{4\pi ^2(1-\sigma ^2)^2}\right ]^{1/3}. \end{align}

That is, fluid removal starts after a time period of $t = t_i$ of hydraulic fracturing under a constant injection rate of $2q$ , which generates an initial fracture of length $x \approx x_c$ , thickness $h \approx h_c$ and pressure $p \approx p_c$ . It is convenient to see that the choice of $t_i$ does not influence the dynamics of fluid removal in the dimensionless system. Scaling results (2.5) are also consistent with those obtained by Spence & Sharp (Reference Spence and Sharp1985). Also, it is interesting on note that lubrication approximation requires that requires that $h_c/x_c \ll 1$ . Based on the scaling results (2.5) for the characteristic length and thickness scales, we obtain an equivalent requirement on the fracturing time that $t_i \gg 6\pi \mu (1-\sigma ^2)/E$ . Physically, this requirement indicates that the initial fracturing time must be sufficiently large to generate a long and thin elastic cavity.

We then define dimensionless variables based on this initial state as

(2.6) \begin{equation} \bar t \equiv t/t_i,\,\,\, \bar x \equiv x/x_c,\,\,\, \bar h \equiv h/h_c\,\,\, \mbox {and}\,\,\, \bar p \equiv p/p_c, \end{equation}

and the model problem (2.1a , b ) with IBCs (2 . 2a –c ) for hydraulic fracturing for the period of $\bar t \in [0,1]$ then becomes

(2.7a) \begin{align} \frac {\partial \bar h}{\partial \bar t} - \frac {\partial }{\partial \bar x}\left (\bar h^3\frac {\partial \bar p}{\partial \bar x}\right ) = 0, \end{align}

(2.7b) \begin{align} \bar p(\bar x,\bar t) = - \int _{-\bar x_f(\bar t)}^{\bar x_f(\bar t)} \frac {1}{s - \bar x}\frac {\partial \bar h}{\partial s}\, {\rm d} s, \end{align}

subject to

(2.8a) \begin{align} \bar h(\bar x,0) &= 0, \end{align}

(2.8b) \begin{align} \bar h(\pm \bar x_f(\bar t), \bar t) &= 0, \end{align}

(2.8c) \begin{align} \int _{-\bar x_f(\bar t)}^{\bar x_f(\bar t)} \bar h(\bar x, \bar t)\, \textit { d}\bar x &= 2\bar t. \end{align}

It could be more convenient to focus only on the positive half-domain since the model problem is symmetric about the axis of $\bar x=0$ . In this case, the elasticity equation (2.7b ) is equivalent to

(2.9) \begin{equation} \bar p(\bar x,\bar t) = - \int _{0}^{\bar x_f(\bar t)} \frac {2s}{s^2 - \bar x^2}\frac {\partial \bar h}{\partial s}\,{\rm d} s, \end{equation}

including appropriately the contribution of elastic deformation in the negative half-domain. Meanwhile, the global volume constraint (2.8c ) is equivalent to a flux boundary condition imposed at $\bar x=0$ as $-\bar h^3 \frac {\partial \bar p}{\partial \bar x}\big |_{\bar x = 0} = 1$ , which is more convenient to impose in finite volume schemes to numerically solve the nonlinear non-local system (2.7a ) and (2.9) for the coupled evolution of $\bar h(\bar x, \bar t)$ and $\bar p(\bar x, \bar t)$ . Indeed, we can integrate the Reynolds equation (2.7a) from $x=0$ to $x=x_f(t)$ , and by applying the volume constraint (2.8c ) and that flux is zero at the front $-\bar h^3 \frac {\partial \bar p}{\partial \bar x}\big |_{\bar x = \bar x_f(\bar t)} = 0$ , we obtain the flux condition at the origin. Similar treatments have been employed in a series of previous studies of gravity current injection into confined porous layers (e.g. Zheng et al. Reference Zheng, Guo, Christov, Celia and Stone2015; Guo et al. Reference Guo, Zheng, Celia and Stone2016) and hydraulic fracturing of elastic solids in Hele-Shaw cells (e.g. Zheng Reference Zheng2024). Equivalently, we then solve (2.7a ) and (2.9) subject to IBCs

(2.10a) \begin{align} \bar h(\bar x,0) &= 0, \end{align}

(2.10b) \begin{align} \bar h(\bar x_f(\bar t), \bar t) &= 0, \end{align}

(2.10c) \begin{align} -\bar h^3 \frac {\partial \bar p}{\partial \bar x}\bigg |_{\bar x = 0} &= 1, \end{align}

for the period of constant fluid injection ( $\bar t \in [0,1]$ ).

By solving the fracturing problem from $\bar t = 0$ to $\bar t = 1$ , we obtain a set of initial conditions $\bar h_i(\bar x) = \bar h(\bar x, 1)$ and $\bar p_i(\bar x) = \bar p(\bar x, 1)$ for the subsequent process of fluid removal, during which time period we focus also on the positive half-domain and deal with the nonlinear non-local system (2.7a ) and (2.9). We choose the same characteristic scales as (2.5 a-c), and the dimensionless IBCs for this period are

(2.11a) \begin{align} \bar h_i(\bar x) &= \bar h(\bar x, 1), \end{align}

(2.11b) \begin{align} \bar h(\bar x_f(\bar t), \bar t) &= 0, \end{align}

(2.11c) \begin{align} -\bar h^3 \frac {\partial \bar p}{\partial \bar x}\bigg |_{\bar x = 0} &= -\Pi _c\,\bar p_0(\bar t) - \Pi _e, \end{align}

where two dimensionless parameters are recognised and defined as

(2.12) \begin{equation} \Pi _c \equiv \frac {E^{2/3} d^3}{[48\pi \mu (1-\sigma ^2)]^{2/3}\,q L\, t_i^{1/3}}, \,\,\,\mbox {and}\,\,\, \Pi _e \equiv -\frac {p_e d^3}{12\mu q L}. \end{equation}

The flux condition (2.11c ) is the dimensionless form of Poiseuille’s law (2.4) for the viscous backflow through the narrow conduit, and note that only half of the flow rate (2.4) through the conduit is considered for fluid removal from the positive half-domain.

Physically, $\Pi _c \gt 0$ measures the relative strength of elasticity-driven fluid removal within a fracture, while $\Pi _e \geqslant 0$ measures the influence of ambient-pressure-driven fluid removal (i.e. the influence of extraction, when $p_e \leqslant 0$ ). Physically, $\Pi _c$ and $\Pi _e$ can also be understood as two different ratios between the three pressure scales in the present problem: the elastic pressure due to deformation of the cavity; the inflation over-pressure due to fluid injection through the conduit; the imposed under-pressure $p_e$ for extraction. In particular, $\Pi _c$ is a pressure ratio, which measures the relative strength of the elastic pressure $ (3\mu E^2/[4\pi ^2(1-\sigma ^2)^2 t_i] )^{1/3}$ due to cavity deformation compared with the inflation over-pressure $\mu \rm qL\,d^-{^3}$ due to fluid injection through the conduit at constant rate $q$ . In contrast, $\Pi _e$ is another pressure ratio, which measures the competition between an imposed under-pressure $p_e$ and the fluid injection induced inflation over-pressure $\mu \rm qL\,d^-{^3}$ through the narrow conduit. The definition and physical meaning of $\Pi _c$ and $\Pi _e$ are summarised in table 1. In particular, $\Pi _e = 0$ represents the situation of fluid removal driven completely by elasticity, while $\Pi _e \gt 0$ represents fluid removal that is also under the influence of extraction.

Table 1.

Definition and physical meaning of the two dimensionless parameters $\Pi _c$ and $\Pi _e$ identified for the process of fracture deflation during fluid withdrawal through a narrow conduit at the centre.

2.3. Numerical solutions

We are now ready to numerically solve the full partial differential–integral system (2.7a ) and (2.9), subject to appropriate IBCs for the particular process of interest (i.e. either fracturing or backflow). The influence of $\Pi _c$ and $\Pi _e$ can also be investigated. In particular, a finite volume scheme is described in appendix A, which allows us to numerically solve (2.7a ) and (2.9), subject to IBCs (2.10) for the period of hydraulic fracturing, or IBCs (2.11a – c ) for the period of fluid withdrawal following fracturing. We note again that in the present work we assume that backflow starts at the moment when injection stops. If there exists any time lag between injection and backflow, other interesting dynamics can possibly occur (e.g. Tanikella et al. Reference Tanikella, Sigallon and Dressaire2023), which could also be of interest for future explorations. Also, with the numerical scheme presented in appendix A, there is no need to introduce a thin prewetting film ahead of the fracture’s front to obtain numerical solutions for the coupled evolution of the profile shape $\bar h(\bar x, \bar t)$ and pressure distribution $\bar p(\bar x, \bar t)$ .

To illustrate the time evolution of the profile shape and pressure distribution during fluid removal, representative numerical solutions are shown in figure 2(a), where we have set $(\Pi _c,\Pi _e) = (1,40)$ , and also in figure 2(b), where we have set $(\Pi _c,\Pi _e) = (1,0)$ for purely elasticity-driven fluid withdrawal without the influence of extraction. In particular, we have observed in figure 2(a) that with extraction ( $\Pi _e \gt 0$ ), the thickness at the location of the conduit ( $\bar x = 0$ ) approaches zero fairly quickly, i.e. $\bar h(0,\bar t) \to 0^+$ , after a fairly short initial transition period. The local pressure at the conduit becomes negative, i.e. $\bar p(0,\bar t) \lt 0$ , due to extraction. But in contrast, without extraction ( $\Pi _e = 0$ ), for the purely elasticity-driven backflow in figure 2(b), the pressure at the location of the conduit approaches zero after an initial transition, i.e. $\bar p(0,\bar t) \to 0^+$ , while the profile thickness $\bar h(0, \bar t)$ approaches zero at a much slower rate (but indeed $\bar h(0, \bar t) \to 0^+$ after a sufficiently long period of time of backflow). Also, in both cases, we have observed the existence of a maximum pressure between the location of the conduit $x=0$ and the fluid front $x=x_f(t)$ , which indicates that the fracturing fluid continues to move forward as $x_f(t)$ increases with time, while it also flows backwards towards the conduit and leads to a decrease of total fluid volume with time. In fact, when we impose different values of $(\Pi _c,\Pi _e)$ with extraction such as $(\Pi _c,\Pi _e) = \{(10,40);(1,20)\}$ , as shown in figure 3, we found that in all cases we tested the thickness approaches zero at the centre after an initial transition period, the length of which depends on the values of $\Pi _c$ and $\Pi _e$ . More details are provided in appendix C on the pinching time $\bar t_d$ for the fracture’s thickness at the centre to approach zero and its dependence on $(\Pi _c,\Pi _e)$ . Eventually, the dynamics of the coupled flow and deformation can be described by a self-similar solution of the dipole type, as we discuss in the following section.

Figure 2.

Numerical solutions for the time evolution of the profile shape $\bar h(\bar x, \bar t)$ and the pressure distribution $\bar p(\bar x, \bar t)$ for different $(\Pi _c,\Pi _e)$ . (a) Numerical solutions at $\bar t = 0.5\times \{0-6\}$ for $(\Pi _c,\Pi _e) = (1,40)$ : the thickness approaches zero at $\bar x = 0$ . (b) Numerical solutions at $\bar t = 0.3\times \{0-5\}$ for $(\Pi _c,\Pi _e) = (1,0)$ : the thickness at $\bar x = 0$ remains finite, while the local pressure approaches zero at a faster rate. The insets show the evolution of the pressure singularity at the fracture’s front, but note that the accuracy in this local region (but not the bulk region) depends significantly on the numerical resolution. We have imposed $\{N,L_d\}=\{300,3\}$ . The red dot–dashed curves represent the initial conditions. Thereafter, the solid curves evolve from cyan to red as time progresses.

Figure 3.

Numerical solutions for the time evolution of the profile shape $\bar h(\bar x, \bar t)$ for (a) $(\Pi _c,\Pi _e) = (1,20)$ and (b) $(\Pi _c,\Pi _e) = (10,40)$ at $\bar t = 0.5\times \{0-6\}$ . The thickness at $\bar x = 0$ approaches zero after an initial transition period, the length of which depends on $\Pi _c$ and $\Pi _e$ . We have imposed $\{N,L_d\}=\{300,3\}$ here.

3.1. Conservation of the dipole moment

As described before, our numerical solutions indicate that, under the influence of extraction (when $\Pi _e \gt 0$ ), the thickness of a fracture eventually approaches zero at the centre, i.e.

(3.1) \begin{equation} h_0(t) \equiv h(0,t) \to 0^+, \end{equation}

as $t \to \infty$ , as shown in figures 2 and 3. This is an important condition, as we show next, that leads to the conservation of the dipole moment at intermediate times. To some extent, such a zero-thickness condition (3.1) is analogous to the edge condition during buoyancy-driven drainage of a viscous gravity current at intermediate times, with experimental observations available from draining flows in vertically placed Hele-Shaw cells from finite-volume releases (e.g. King & Woods Reference King and Woods2003; Zheng et al. Reference Zheng, Soh, Huppert and Stone2013). There exists also a similar situation in the context of seasonal thermal energy storage, when the working fluid is extracted from a propagating gravity current through a narrow conduit at the origin (e.g. Dudfield & Woods Reference Dudfield and Woods2012, Reference Dudfield and Woods2014). In their mimicking laboratory experiments, it is also observed that eventually the thickness of the fluid film approaches zero at the location of the conduit.

We can also show that the dipole moment approaches a constant at this stage,

(3.2) \begin{equation} M(t) \equiv \int _0^{x_f(t)} xh(x,t)\,{\rm d}x \to M_{\infty }, \end{equation}

where $M_{\infty }$ is a constant (with dimension [ $L^3$ ]) that depends possibly on $(\Pi _c,\Pi _e)$ . In particular, we can show that the time derivative of the dipole moment goes to zero at intermediate times. That is, assuming that $h_0(t) \to 0^+$ and $h(x_f(t),t) = 0$ as $t \to \infty$ , as supported by the representative numerical solutions with extraction ( $\Pi _e\gt 0$ ), it is then seen that

(3.3) \begin{align} &\frac {{\rm d}}{{\rm d}t} \int _0^{x_f(t)} x h(x,t)\,{\rm d}x = x_f(t)\,h(x_f,t)\frac { {\rm d}x_f(t)}{ {\rm d}t} - \int _0^{x_f(t)} x \frac {\partial h}{\partial t} {\rm d}x\nonumber \\ &= - \frac {1}{3\mu }\int _0^{x_f(t)} x \frac {\partial }{\partial x} \left (h^3\frac {\partial p}{\partial x}\right ) \, {\rm d}x = - \frac {1}{3\mu } \left [ xh^3\frac {\partial p}{\partial x}\bigg |_{0}^{x_f(t)} - \int _0^{x_f(t)} h^3 \frac {\partial p}{\partial x}{\rm d}x \right ] \nonumber \\ &= \frac {1}{3\mu } \int _0^{x_f(t)} h^3 \frac {\partial p}{\partial x}\,{\rm d}x = \frac {1}{3\mu }\left [ph^3 \bigg |_0^{x_f(t)} - \int _{h_0^3(t)}^0 p(x(h,t),t)\, {\rm d}h^3 \right ]\quad \quad \quad \quad \,\, \nonumber \\ &\sim -\frac {1}{3\mu } \int _0^{0} p(x(h,t),t)\, {\rm d}h^3 \sim - \frac {p_1(t)}{3\mu } \int _0^{0} p_2(h)\, {\rm d}(h^3) = 0,\quad \quad \quad \quad \quad \quad \quad \,\,\,\,\, \end{align}

where we have assumed that $p$ and $\partial p/\partial x$ are not too singular at $x = 0$ and $x = x_f(t)$ , such that $xh^3\partial p/\partial x|_{x=0} = 0$ and $xh^3\partial p/\partial x|_{x=x_f(t)} = 0$ . Also, the last step is based on the assumption of the existence of self-similar behaviours at intermediate times (as $t \to \infty$ ) such that $p(x(h,t),t) \sim p_1(t)\,p_2(h)$ . We later-on self-consistently confirm this statement also by comparing the predictions of a dipole self-similar solution with time-dependent numerical solutions. In particular, we have tested a range of $(\Pi _c, \Pi _e)$ , and in all cases we found that the dimensionless dipole moment $\bar M(\bar t)$ indeed approaches a constant after an initial transition period, as shown in figure 4. In fact, this seems also a good approximate even without extraction (when $\Pi _e = 0$ ), but in which case the time transition is longer.

Figure 4.

Numerical solutions for the time evolution of the dipole moment $\bar M(\bar t)$ for different values of $(\Pi _c,\Pi _e)$ . In all cases $\bar M(\bar t)$ appears to approach a constant after an initial transition time period. We have imposed $\{N,L_d\}=\{300,3\}$ in all of the calculations for a sense of numerical error.

Such a feature is again similar to the dipole flow of a viscous gravity current, with experimental evidence provided from dipole currents in vertically placed Hele-Shaw cells that the dipole moment $M(t)$ indeed approaches a constant at intermediate times (e.g. King & Woods Reference King and Woods2003; Mathunjwa & Hogg Reference Mathunjwa and Hogg2007). Again, all these features further demonstrate that the physics of such a free interface problem of coupled flow and deformation in an infinite elastic solid is closely related to that of a viscous gravity current, when the influence of material toughness is negligible.

Further, scaling analysis of (2.1a , b ) subject to boundary conditions (2.2b ) and (3.2) (rather than (2.2c ) for hydraulic fracturing) leads to a distinct dipole scaling as

(3.4) \begin{equation} x \propto t^{1/9},\,\,\,h \propto t^{-2/9}\,\,\, \mbox {and}\,\,\,p \propto t^{-1/3}, \end{equation}

at intermediate times. Accordingly, the total fluid volume $V(t)$ within a fracture (in the positive half-domain) decreases with time and so is the rate of backflow $q(t)$ , which now scale with time as

(3.5) \begin{equation} V \equiv \int _{0}^{x_f(t)}h(x,t)\, {\rm d}x \propto t^{-1/9}\,\,\, \mbox {and}\,\,\, q \propto t^{-10/9}. \end{equation}

Self-similar solutions can also be explored for $h(x,t)$ and $p(x,t)$ by defining a similarity variable $\xi \equiv x/t^{1/9}$ , which we discuss later. These features might be of interest in practical contexts such as subsurface environment remediation: The thickness of a fracture is related to the effective permeability of a fractured porous medium, and the rate of backflow is related to the rate of release of pollutants and chemicals potentially brought back from the porous rock systems.

3.2. Self-similar solution of the dipole type

We now construct self-similar solutions of the dipole type during fracture deflation at intermediate times. Our numerical solutions indicate that this situation is arrived at sufficiently late times, in particular, when $\Pi _e \gt 0$ , when the influence of extraction becomes strong enough. Here in § 3.2 we focus only on the situation with sufficient extraction ( $\Pi _e \gt 0$ ), while in § 3.4 we discuss the situation of purely elasticity-driven backflow without extraction ( $\Pi _e=0$ ). Again we focus only on the positive half-domain since the problem is symmetric about the $z$ axis. Knowing that the thickness of a fracture approaches zero at the centre, we then solve the dimensionless system (2.7a ) and (2.9) subject to boundary conditions

(3.6) \begin{equation} \bar h(0, \bar t) = 0, \,\,\,\mbox {and}\,\,\, \bar h(\bar x_f(\bar t), \bar t) = 0. \end{equation}

By defining a similarity variable of $\xi \equiv \bar x/\bar t^{1/9}$ , we look for $\bar h(\bar x, \bar t) = \bar t^{-2/9} k(\xi )$ and $\bar p(\bar x, \bar t) = \bar t^{-1/3} g(\xi )$ . The governing partial differential–integral system (2.7a ) and (2.9) is then rearranged to provide an ordinary differential–integral system of

(3.7a) \begin{align} \frac {2}{9} k +\frac {1}{9} \xi \frac {{\rm d}k}{\textit { d}\xi } + \frac {{\rm d}}{{\rm d} \xi }\left (k^3\frac {{ \rm d} g}{ {\rm d} \xi }\right ) = 0, \end{align}

(3.7b) \begin{align} g(\xi ) = - \int _{0}^{\xi _f} \frac {2s}{s^2-\xi ^2}\frac {{\rm d} k}{ {\rm d} s} {\rm d} s, \end{align}

where $\xi _f \equiv \bar x_f(\bar t)/\bar t^{1/9}$ is a frontal constant. System (3.7) is then solved with boundary conditions,

(3.8) \begin{equation} k(0) = 0, \,\,\,\mbox {and}\,\,\, k(\xi _f) = 0. \end{equation}

Eventually, $k(\xi )$ and $g(\xi )$ are self-similar solutions of the dipole type for the coupled evolution of the profile shape and pressure distribution during fluid removal at intermediate times. A numerical scheme is provided in appendix B to solve (3.7a , b ), which can also be used to study the stability of these self-similar solutions.

It is found that $k(\xi )$ and $g(\xi )$ are nonlinear attractors indeed for the initial conditions we imposed. The normalised self-similar solutions $k(\xi )/k_m$ and $g(\xi )/g_m$ we obtained are shown in figure 5 as the dashed curves, where $k_m$ and $g_m$ are the corresponding maximum values. In addition, in figure 5, the raw numerical solutions of $\bar h(\bar x, \bar t)$ and $\bar p(\bar x, \bar t)$ in figure 2(a) when $(\Pi _c,\Pi _e) = (1,40)$ are rescaled based on the maximum values $\bar h_m(\bar t)$ and $\bar p_m(\bar t)$ , which demonstrates very good data collapse for both the profile shape and pressure distribution after an initial period of time transition and also very good agreement with the normalised self-similar solutions $k(\xi )/k_m$ and $g(\xi )/g_m$ . We also performed the same analysis for the time-dependent numerical solutions when $(\Pi _c,\Pi _e) = (1,20)$ and $(\Pi _c,\Pi _e) = (10,40)$ , as shown in figure 3, and arrive at the same conclusion. The comparison is similar to figure 5 for $(\Pi _c,\Pi _e) = (1,40)$ and is not shown here.

Figure 5.

The rescaled numerical solutions for the time evolution of both the profile shape $\bar h(\bar x, \bar t)/\bar h_m(\bar t)$ in (a) and pressure distribution $\bar p(\bar x, \bar t)/\bar p_m(\bar t)$ in (b) as a function of dimensionless horizontal location $\eta = \bar x/\bar x_f(\bar t)$ for $(\Pi _c,\Pi _e) = (1,40)$ at $\bar t = 0.5\times \{1-6\}$ . The time dependent numerical solutions collapse onto a universal curve eventually, which also agrees reasonably well with the prediction of the self-similar solutions from solving (3.7) subject to (3.8).

The model problem (3.7a , b ) can also be rearranged to provide more insights by further defining $\eta \equiv \xi /\xi _f \in [0,1]$ , $\mathcal {K}(\eta ) = \xi _f\, k(\xi /\xi _f)$ and $\mathcal {G}(\eta ) = g(\xi /\xi _f)$ . Then, (3.7a , b ) and (3.8) can be rewritten as

(3.9a) \begin{align} \frac {2}{9} \mathcal {K} +\frac {1}{9} \eta \frac {{\rm d}\mathcal {K}}{{\rm d}\eta } + \frac { {\rm d}}{ {\rm d} \eta }\left (\mathcal {K}^3\frac {{\rm d} \mathcal {G}}{ {\rm d} \eta }\right ) = 0, \end{align}

(3.9b) \begin{align} \mathcal {G}(\eta ) = - \int _{0}^{1} \frac {2s}{s^2-\eta ^2}\frac {{\rm d} \mathcal {K}}{ {\rm d} s}{\rm d} s \end{align}

and

(3.10) \begin{equation} \mathcal {K}(0) = 0, \,\,\,\mbox {and}\,\,\, \mathcal {K}(1) = 0. \end{equation}

In this case, conservation of the dipole moment (3.2) can be rearranged to provide an explicit solution for the frontal prefactor $\xi _f$ ,

(3.11) \begin{equation} \xi _f = \bar M_{\infty }^{1/3} \left [\int _0^{1} \eta \, \mathcal {K}(\eta )\,{\rm d}\eta \right ]^{-1/3}, \end{equation}

where $\bar M_{\infty } \equiv M_{\infty }/(h_cx_c^2)$ . It is seen that the prefactor $\xi _f$ depends possibly on the initial condition and the time transition. This is because the dimensionless dipole moment $\bar M_{\infty }$ , a constant at late times, could be under the influence of the initial condition and the values of $\Pi _c$ and $\Pi _e$ . When $\Pi _c$ and $\Pi _e$ are the same, the value of $\xi _f$ then measures the difference between the initial profile shape and the self-similar dipole shape $\mathcal {K}(\eta )$ . As noted before, this is because the dipole-type self-similar solutions are incomplete self-similar solutions, for which the scaling exponents are independent on the initial condition and time transition, but the prefactors can be impacted by them (e.g. Barenblatt Reference Barenblatt1996). Furthermore, with the problem reformulated as (3.9a , b ), it can be shown that we must have

(3.12) \begin{equation} \mathcal {K} \propto (1-\eta )^{2/3},\,\,\,\mbox {and}\,\,\, \mathcal {G} \propto -(1-\eta )^{-1/3}, \end{equation}

as $\eta \to 1^-$ , following the same scaling analysis in Spence & Sharp (Reference Spence and Sharp1985). That is, the local structure near the fracture’s front is the same for either fracturing or dipole flows.

3.3. Comparison with time-dependent numerical solutions

Time transition towards the self-similar solutions of the dipole type can also be illustrated based on the time dependent numerical solutions of (2.7a , b ), subject to IBCs (2.11a - c ), as shown in figures 2– 5. For the parameter range of $(\Pi _c,\Pi _e)$ we tested, it is shown that $h(0,t) \to 0^+$ eventually for all cases with extraction $\Pi _e \gt 0$ . Meanwhile, the values of $(\Pi _c,\Pi _e)$ can influence significantly the exact time when the fracture thickness approaches zero at the centre, i.e. the time needed for $h(0,t) \to 0^+$ .

Figure 6.

Numerical solutions for the time evolution of the frontal location $\bar x_f(\bar t)$ , maximum thickness $\bar h_m(\bar t)$ , maximum pressure $\bar p_m(\bar t)$ and total fluid volume $\bar V(\bar t)$ . We have included three cases with $(\Pi _c,\Pi _e) = \{(1,20),(1,40),(10,40)\}$ considered here with extraction, which exhibit the same scaling laws at late times. We have imposed $\{N,L_d\}=\{300,3\}$ in all simulations.

We can also track the dipole moment $\bar M(\bar t)$ , as shown in figure 4, based on the time-dependent numerical solutions, in addition to the frontal location $\bar x_f(\bar t)$ , maximum thickness $\bar h_m(\bar t)$ , maximum pressure $\bar p_m(\bar t)$ and total fluid volume $\bar V(\bar t)$ , as shown in figure 6. Following an initial transition period they all approach the power-law scaling laws of dipole similarity, as identified in (3.2), (3.4) and (3.5),

(3.13) \begin{equation} \bar x_f \propto \bar t^{1/9},\,\,\, \bar h_m \propto \bar t^{-2/9},\,\,\, \bar p_m \propto \bar t^{-1/3},\,\,\,\bar V \propto \bar t^{-1/9}\,\,\,\mbox {and}\,\,\, \bar M \propto 1. \end{equation}

The simulation ends at around $\bar t \approx 4$ for practical reasons, and we expect better agreement at later times, allowing more time to finish the time transition. It is also of interest to note that the prefactors, nevertheless, depend possibly on $(\Pi _c,\Pi _e)$ and the initial condition. This is as expected, because self-similar solutions of the dipole type are incomplete self-similar solutions (e.g. Barenblatt Reference Barenblatt1996).

3.4. Purely elasticity-driven backflow

The special case of $\Pi _e = 0$ corresponds to the backflow of the fracturing fluid driven completely by elasticity and without any external effects of extraction. Our time-dependent numerical solutions indicate that both the fracture’s thickness $\bar h(0,\bar t)$ and the local pressure $\bar p(0,\bar t)$ at the centre approach zero as time progresses, while $\bar p(0,\bar t)$ approaches zero at a faster rate. This aspect is partly consistent with experimental observations during flowing back from penny-shaped hydraulic fractures in solid gelatins (Lai et al. Reference Lai, Zheng, Dressaire, Ramon, Huppert and Stone2016a ), even though the front remains approximately pinned in those backflow experiments due to the influence of material toughness. However, this is different from the dipole situation (3.1) under the influence extraction, as discussed in this section, and the self-similar solution from solving (3.7) subject to (3.8) does not apply.

Figure 7.

Numerical solutions for the time evolution of the frontal location $\bar x_f(\bar t)$ , maximum thickness $\bar h_m(\bar t)$ , maximum pressure $\bar p_m(\bar t)$ and the total fluid volume $\bar V(\bar t)$ . We have included three cases with $\Pi _c = \{1,10,100\}$ considered here for purely elasticity-driven fracture deflation without extraction ( $\Pi _e=0$ ), which exhibit the same scaling laws at late times. We have imposed $\{N,L_d\}=\{300,3\}$ in all simulations.

Nevertheless, when we studied numerical solutions of the time-dependent system (2.7a ) and (2.9) for purely elasticity-driven backflow when $\Pi _e = 0$ , we found that the dipole moment $M$ also approaches a constant for all cases we considered, as shown in figure 4. This seems plausible, as one can check that all steps of the derivation (3.3) that lead to the conservation of the dipole moment continue to apply, providing that the simulation time is sufficiently long, such that the condition of zero thickness at the centre $\bar h(0,\bar t) \to 0^+$ applies approximately even without extraction. It is also observed in figure 4 that the lengths of time transition when $\Pi _e = 0$ seem to be longer than those when $\Pi _e \gt 0$ , and this seems also reasonable, because without extraction it takes much more time for the necessary condition of $\bar h(0,\bar t) \to 0^+$ to apply.

In addition, the scaling results of (3.4) also appear to apply for the time evolution of the frontal location $\bar x_f(\bar t)$ , maximum thickness $\bar h_m(\bar t)$ , maximum pressure $\bar p_m(\bar t)$ as well as the total fluid volume $\bar V(\bar t)$ , as shown in figure 7. It is also observed in figure 7 that within the time scale of the numerical simulation, the degree of agreement of the slopes (in the log–log plot) with the prediction of the scaling theory (3.4) seems to be better for greater $\Pi _c$ . This seems reasonable, because a greater $\Pi _c$ corresponds to a shorter time period required for the fracture’s thickness at the centre $\bar h(0,\bar t)$ to approach zero, which is a necessary condition for the dipole moment $M$ to be considered as a constant and for all of the scaling results in (3.4) to be applicable. There is also the observation that the degree of agreement of the slopes with the scaling theory (3.4) in figure 7 seems worse than that in figure 6 for cases with extraction (when $\Pi _e \gt 0$ ). This is also as expected for the same reason that with extraction it takes a shorter amount of time for the condition of $\bar h(0,\bar t) \to 0^+$ to apply.

Figure 8.

The normalised numerical solutions for the time evolution of both the profile shape $\bar h(\bar x, \bar t)/\bar h_m(\bar t)$ and pressure distribution $\bar p(\bar x, \bar t)/\bar p_m(\bar t)$ as a function of dimensionless horizontal location $\eta = \bar x/\bar x_f(\bar t)$ for $(\Pi _c,\Pi _e) = (1,0)$ at $\bar t = 0.5\times \{1-8\}$ . The normalised time dependent numerical solutions collapse onto a universal curve eventually, which also agrees reasonably well with the prediction of the self-similar solutions. We have imposed $\{N,L_d\}=\{200,1.5\}$ for the self-similar solutions provided here.

Furthermore, when we study numerical solutions for the profile shape $\bar h(\bar x, \bar t)$ and pressure distribution $\bar p (\bar x, \bar t)$ , we found that the normalised/rescaled numerical solutions in figure 2(b) for $(\Pi _c,\Pi _e) = (1,0)$ still approach universal solutions at intermediate times, as shown in figure 8. These observations indicate further the existence of another self-similar regime for purely elasticity-driven backflow ( $\Pi _e=0$ ).

To construct this new set of self-similar solutions, based on scaling results (3.4), we define the same similarity transform $\bar h(\bar x, \bar t) = \bar t^{-2/9} k(\xi )$ and $\bar p(\bar x, \bar t) = \bar t^{-1/3} g(\xi )$ , where $\xi \equiv \bar x/\bar t^{1/9}$ . The self-similar solutions for the profile shape $k(\xi )$ and pressure distribution $g(\xi )$ can then be obtained by solving the same ordinary differential–integral system (3.7a , b ), now subject to a different set of boundary conditions

(3.14) \begin{equation} g(0) = 0, \,\,\,\mbox {and}\,\,\, k(\xi _f) = 0, \end{equation}

which indicate zero pressure at the origin rather than zero thickness (as used in (3.8) when $\Pi _e \gt 0$ ) in the similarity space. The corresponding self-similar (monopole) solutions $k_m^{-1}k(\xi )$ and $g_m^{-1}g(\xi )$ are included in figure 8 as the dashed curves, employing the same numerical scheme described in appendix B. We observe very good agreement with the normalised time-dependent numerical solutions. Notably, these self-similar solutions appear to be independent of $\Pi _c$ , the value of which impacts only the transient dynamics. This is as expected because the governing system at intermediate times does not include any information of $\Pi _c$ . We also found numerically that the normalised inlet thickness $k_m^{-1}k(0) \approx 0.89$ (with 200 space grids imposed).

4.1. Summary

A theoretical model was developed to describe the process of fracture deflation, following an initial period of hydraulic fracturing of an infinite elastic solid. The model is based on a nonlinear partial differential–integral system for the coupled process of viscous thin film flow and pressure evolution with elastic deformation. The novel component for the model problem of fracture deflation is to couple the process of viscous backflow through a narrow conduit at the centre of the fracture. This is achieved by imposing a time dependent flux boundary condition at the centre that depends also on the local pressure evolution within the fracture. Two dimensionless parameters ( $\Pi _c$ and $\Pi _e$ ) were identified after appropriate non-dimensionalisation of the model problem, as summarised in table 1. They describe the relative influence of elasticity-driven and ambient-pressure-driven backflows, respectively.

Figure 9.

Self-similar solutions from hydraulic fracturing under constant injection rates towards fracture deflation in the viscous regime: (a) dipole solutions with extraction ( $\Pi _e \gt 0$ ) as discussed in § 3.2, and (b) monopole solutions without extraction ( $\Pi _e = 0$ ) as discussed in § 3.4.

Time dependent numerical solutions indicate that, under the influence of ambient-pressure-driven extraction (when $\Pi _e\gt 0$ ), the fracture’s thickness at the centre approaches zero after an initial period of time transition. In this situation, self-similar solutions of the dipole type can be obtained to describe the key universal feature of the coupled evolution of the profile shape and pressure distribution at intermediate times, when the fracture’s thickness becomes zero at the centre. Notably, the exact value of $\Pi _c$ and $\Pi _e$ do not influence the normalised self-similar solutions, but they influence the length of time transition and the frontal prefactor $\xi _f$ . This is because self-similar solutions of the dipole type are incomplete self-similar solutions. In contrast, for purely elasticity-driven fracture deflation (when $\Pi _e=0$ ), the elastic pressure at the centre approaches zero instead, but the dipole moment is also conserved at intermediate times. We have obtained another set of self-similar solutions to describe the universal profile shape and pressure distribution in this situation. Key features of the fracture deflation process are summarised in figure 9 in both regimes of dipole solutions when $\Pi _e\gt 0$ (with extraction) and monopole solutions when $\Pi _e=0$ (without extraction), and also in table 2 with an emphasis on their time dependence during the time evolution.

Table 2.

Key features of fracture deflation at intermediate times during fluid removal through a narrow conduit at the centre, when the influence of material toughness is negligible. These scaling results apply for backflow with extraction $(\Pi _e \gt 0)$ and without extraction $(\Pi _e = 0)$ for sufficiently large times.

4.2. The influence of material toughness

For hydraulic fracturing of an elastic and brittle solid, it is well understood that besides the coupled description of lubricating flow and elastic deformation, the propagation of a fracture must also satisfy an additional stress condition that is related to a material property named material toughness $K_{c}$ (e.g. Spence & Sharp Reference Spence and Sharp1985; Lister Reference Lister1990a ),

(4.1) \begin{equation} K \equiv \frac {E}{2^{3/2}(1-\sigma ^2)}\frac {h}{(x_f-x)^{1/2}} = K_{c},\,\,\,\mbox {as}\,\,\, x\to x_f^-. \end{equation}

Physically, (4.1) indicates that to break a material at the tip (for a crack to propagate), the stress integral factor $K$ must be large enough to balance material toughness $K_{c}$ . This is an important feature of hydraulic fracturing in brittle materials, which could impact significantly the fracturing dynamics. The influence of material toughness on the dynamics of hydraulic fracturing has also been discussed in previous studies with (e.g. Spence et al. Reference Spence, Sharp and Turcotte1987; Lister Reference Lister1990a ,Reference Lister b ; Roper & Lister Reference Roper and Lister2007) and without (e.g. Spence & Sharp Reference Spence and Sharp1985; Detournay Reference Detournay2004) the influence of buoyancy. In particular, in a neutrally buoyant environment, the influence of flow becomes negligible in the toughness dominated regime, and a fracture propagates in a quasisteady fashion that maintains always an elliptic shape (e.g. Howell et al. Reference Howell, Kozyreff and Ockendon2008). The influence of material toughness on fracture deflation is discussed in another work (Zheng Reference Zheng2025).

4.3. Potential implications

The model and solutions can possibly find applications, for example, in the context of fluid-phase resource recovery. In the field it is sometimes observed that a fracture eventually closes at the well bore after pressure release following fracturing, which is then consistent with the prediction of the current work. In this case, it is of particular interest to understand the time scale of fluid removal and to what extent the fracturing fluid can drain, since they are related to the rate of oil and gas recovery and that of pollutants and chemicals being carried out from subsurface porous rocks. They are all addressed in the current work in the context of fluid flow and elastic deformation. The time-dependent numerical solutions and the self-similar solutions can both find use for appropriate time periods, using appropriate geophysical and operational parameters.

To further illustrate the idea, we provide an example calculation here. We consider fracturing a porous rock of Young’s modulus $E = 15$ GPa, Poisson ratio $\sigma = 0.33$ and the material toughness is assumed to be negligibly small here. The rock is initially fractured with water and we are interested in how a fracture deflates and closes once injection ceases and how the fracturing fluid flows back. The geophysical and operation parameters (e.g. Lister Reference Lister1990b ; Savitski & Detournay Reference Savitski and Detournay2002; Bunger et al. Reference Bunger, Lau, Wright and Schmidt2023) are shown in table 3. The only difference of this example calculation from that in another work (Zheng Reference Zheng2025) is that we have now neglected the influence of material toughness here (i.e. imposing $K_{c} = 0$ ), and an interested reader can have a comparison of these results to obtain insights in different deflation regimes. The initial fracturing process lasts for 150 mins while backflow also proceeds for 150 mins in this model problem. Two cases are considered with and without external extraction, and the dimensionless parameters are $(\Pi _c,\Pi _e) = (0.096,0)$ for case (i) and $(\Pi _c,\Pi _e) = (0.096,14.9)$ for case (ii), respectively. The predictions are summarised also in the table, including the final fracture length and thickness, and the maximum and inlet pressures. In particular, it is seen that with extraction the area covered by the fracture is much smaller, which indicates that more fracturing fluid has been removed. However, with extraction the fracture already closes at the centre/conduit by the end of the backflow process at $\bar t = 1$ , which would significantly reduce the backflow rate to continue to remove the fracturing fluid.

Table 3.

Example calculations of fluid withdrawal from an initial fracture in a rock formation: (i) without extraction, and (ii) with extraction. *The external pressure $p_e$ is estimated based on the atmospheric pressure. We have imposed $N=300$ for a domain of length $L_d=3$ here for a sense of numerical error.

Finally, it is good also to note that fluid drainage (or leakage) into the elastic solid might be relevant in practical contexts and leads possibly to new regimes of fracture deflation (e.g. Garagash et al. Reference Garagash, Detournay and Adachi2011; Peirce Reference Peirce2022), as analogous to drainage from gravity currents (e.g. Acton et al. Reference Acton, Huppert and Worster2001; Pritchard et al. Reference Pritchard, Woods and Hogg2001; Yu et al. Reference Yu, Zheng and Stone2017; Liu et al. Reference Liu, Zheng and Stone2017). This problem might also be related to brain blood clot removal. Biological tissues such as the brain tissue might possibly be modelled as elastic materials for their mechanical responses (e.g. Goriely et al. Reference Goriely, Geers, Holzapfel, Jayamohan, Jerusalem, Sivaloganathan, quier, van Dommelen, Waters and Kuhl2015). In the context of brain blood clot removal, it is then of interest to understand how fast and to what extent blood drains. In particular, the model predicts that, while it is possible for a significant amount of the fluid to be removed through the conduit, it is impossible for the fluid to drain out within a finite time. Nevertheless, we also recognise that fluid absorption could possibly lead to new possibilities.