2.1. Governing equations

We first take into consideration the compressibility effect of the foam flow. Denoting the profile shape as $h(x,t)$ for the film of foam, which is also the profile shape of the foam-driven fracture, the requirement of local continuity is written as

(2.1) \begin{equation} \frac {\partial }{\partial t} \left (\int _0^{h(x,t)} \rho (x,y,t) \,\textrm{d}y\right ) + \frac {\partial }{\partial x} \left (\int _0^{h(x,t)} \rho (x,y,t)\, u(x,y,t)\,\textrm{d}y\right ) = 0, \end{equation}

where $\rho (x,y,t)$ and $u(x,y,t)$ denote the density and the velocity fields, respectively. The profile shape $h(x,t)$ of a fracture is also termed the half ‘fracture opening profile’ or half ‘fracture aperture profile’ in the field or in the solid mechanics literature. It is also commonly denoted by $w(x,t)$ . If we further neglect the influence of buoyancy, which is plausible in laboratory experiments since the film of foam is typically mm in thickness (e.g. Lai et al. Reference Lai, Rallabandi, Perazzo, Zheng, Smiddy and Stone2018), then the density field reduces to $\rho (x,t)$ without the $y$ dependence. Since foam flow typically evolves into a plug flow pattern with wall slip in one-dimensional tubes (e.g. Denkov et al. Reference Denkov, Subramanian, Gurovich and Lips2005, Reference Denkov, Tcholakova, Golemanov, Subramanian, Gurovich and Lips2006; Lai et al. Reference Lai, Rallabandi, Perazzo, Zheng, Smiddy and Stone2018), the velocity field also reduces to $u(x,t)$ without the $y$ dependence. In this case, (2.1) further reduces to

(2.2) \begin{equation} \frac {\partial }{\partial t} \left [\rho (x,t)\,h(x,t) \right ] + \frac {\partial }{\partial x} \left [ \rho (x,t)\, u(x,t)\, h(x,t) \right ] = 0. \end{equation}

Compared with the incompressible version of the local continuity equation, we now include in (2.2) an additional weighting function of the density of foam $\rho (x,t)$ in both the time and space derivatives. This form has been employed in a previous study of gas-driven fracturing (e.g. Wang et al. Reference Wang, Chen, Jin and Bunger2018).

The (net) pressure distribution due to elastic deformation is given by the standard integral expression

(2.3) \begin{equation} p(x,t) = - \frac {E}{2(1-\sigma ^2)}\frac {1}{\pi }\int _0^{x_f(t)} \frac {2s}{s^2-x^2}\frac {\partial h}{\partial s}\, \textrm{d}s, \end{equation}

taking into account the contribution from both the negative and positive halves of the domain. Here, $E$ is the Young’s modulus of the elastic solid, and $\sigma$ is the Poisson ratio (often denoted by $\nu$ in the solid mechanics literature). This remains the same as the standard description of hydraulic fracturing using Newtonian fluids (e.g. Spence & Sharp Reference Spence and Sharp1985). The integral expression (2.3) can be obtained by analysing the normal stress distribution at any specific location $x$ , as induced by displacement continuities in the regime of linear elastostatics without the influence of body forces, followed by a superposition of the contribution of all of them throughout the entire domain (see e.g. the Appendix of Peirce Reference Peirce2025). We note again that $p(x,t)$ is the net pressure that is measured on the basis of any externally imposed far-field normal stress $\sigma _\infty$ (e.g. Lister Reference Lister1990b ; Detournay Reference Detournay2004; Garagash Reference Garagash2006), and by definition the atmospheric pressure is already absorbed in $\sigma _\infty$ . Detailed discussions are also available in Garagash (Reference Garagash2006) on the influence of $\sigma _\infty$ , which becomes negligibly small at sufficiently large times, and the transient time depends on its value.

For an ideal gas (such as air or CO $_2$ ), a linear law $\rho = p/(RT)$ typically applies between density and pressure, where $T$ is temperature, and $R$ is a gas constant. We consider only the isothermal system in this work. For foam, we employ the leading-order Taylor expansion next to a reference state of pressure $p_0$ and density $\rho _0$ as

(2.4) \begin{equation} \rho = \rho _0 + \frac {\partial \rho }{\partial p}\bigg |_{p_0} (p-p_0) + \mbox{h.o.t.} \sim \rho _0 + c^{-2} (p-p_0), \end{equation}

where the speed of sound $c$ is introduced. To be self-consistent, the pressures $p(x,t)$ and $p_0$ in (2.4) are also the net pressure on the basis of the far-field normal stress $\sigma _\infty$ . We expect that (2.4) provides a good estimate when the operation is around the reference pressure $p_0$ , i.e. in the weakly compressible regime. Too far away from $p_0$ , more appropriate nonlinear forms should be considered. Based on the form of (2.4), we see that the density $\rho (x,t)$ is dependent on the fracture’s profile shape $h(x,t)$ . A similar linear law has been employed recently to study the influence of flow compressibility on gas-driven fracturing (e.g. Wang et al. Reference Wang, Chen, Jin and Bunger2018).

For foam flow in a straight tube, plug flow patterns would typically develop due to wall slip. The motion of the plug of foam is lubricated by a thin liquid film attaching to the wall. In the bulk region away from the wall, the foam would advance uniformly with a constant slip velocity $u$ . The following result has been employed previously (e.g. Denkov et al. Reference Denkov, Subramanian, Gurovich and Lips2005, Reference Denkov, Tcholakova, Golemanov, Subramanian, Gurovich and Lips2006; Lai et al. Reference Lai, Rallabandi, Perazzo, Zheng, Smiddy and Stone2018) for the magnitude of the slip velocity:

(2.5) \begin{equation} |u| = \frac {r^2 a^2}{\gamma \mu \kappa ^2} \left (\frac {\partial p}{\partial x}\right )^2, \end{equation}

where $r$ is the (average) bubble size, $a$ is the radius of the tube, $\mu$ is the viscosity of the liquid, $\gamma$ is the interfacial tension, and $\kappa (\epsilon )$ is a (dimensionless) resistance constant that depends on the liquid volume fraction $\epsilon$ . For a sense of scales, typical values for the parameters are summarised in table 1 when fracturing a solid gelatine with Gillette foam as the fracturing fluid (Lai et al. Reference Lai, Rallabandi, Perazzo, Zheng, Smiddy and Stone2018). There is also experimental evidence of (2.5); see e.g. figure 8 of Denkov et al. (Reference Denkov, Subramanian, Gurovich and Lips2005), or figure 3 in the supplementary information of Lai et al. (Reference Lai, Rallabandi, Perazzo, Zheng, Smiddy and Stone2018). Compared with the laminar flow of incompressible Newtonian fluids in a tube where $|u| \propto \partial p/\partial x$ , the plug flow of foam now follows $|u| \propto (\partial p/\partial x)^2$ .

Table 1.

Typical parametric values for the fracturing process of a block of solid gelatine using a Gillette foam as the fracturing fluid in a laboratory experiment (e.g. Lai et al. Reference Lai, Rallabandi, Perazzo, Zheng, Smiddy and Stone2018).

An important assumption of solution (2.5) is that the capillary number is sufficiently small (i.e. $Ca \equiv \mu\, |u|/\gamma \ll 1$ ) such that the wall shear stress follows ${\tau _w}/({\gamma /r}) \propto Ca^{1/2}$ . Based on the stress balance $\tau _w = -({a}/{2}) ({\partial p}/{\partial x})$ , we then arrive at solution (2.5), introducing a resistance constant $\kappa (\epsilon )$ . The transient dynamics towards this steady state (2.5) has also been neglected. We must note that the result $\tau _w \propto Ca^{1/2}$ for the wall shear stress is different from the classic result $\tau _w \propto Ca^{2/3}$ for standard viscous capillary films. The physical picture of the latter corresponds to a propagating Taylor bubble in straight tubes (e.g. Bretherton Reference Bretherton1961). Brief comments are also provided in Appendix B regarding the leading-order physics behind the $\tau _w \propto Ca^{1/2}$ and $\tau _w \propto Ca^{2/3}$ behaviours for the wall shear stress.

To describe the film flow of foam inside a fracture of half-thickness $h(x,t)$ , we modify (2.5) by replacing the radius $a$ of the tube by $h(x,t)$ , which leads to

(2.6) \begin{equation} |u(x,t)| = \frac {r^2 h(x,t)^2}{\gamma \mu \hat \kappa } \left (\frac {\partial p}{\partial x}\right )^2 \end{equation}

for the magnitude of the slip velocity. The two major physical components remain the same: The wall shear stress follows ${\tau _w}/({\gamma /r}) \propto Ca^{1/2}$ for the plug flow of foam, and the stress balance follows $\tau _w = -({h}/{2})({\partial p}/{\partial x})$ in the confined space of half-opening $h(x,t)$ . Since $h(x,t)$ is now allowed to vary slowly with $t$ and $x$ , for the latter we have also made the assumption of lubrication approximation. We also recognise the difference in configuration between (2.5) (for tube flow) and (2.6) (for planar flow), so we expect appropriate corrections being provided by a different resistance constant $\hat \kappa (\epsilon )$ . We also neglect the influence of clogging/jamming for viscoplastic flow in converging conduits (e.g. Rallabandi et al. Reference Rallabandi, Zheng, Winton and Stone2017). To understand more about the theory of viscoplastic thin film flows and the appearance of pseudo-yield surfaces, an interested reader can also refer to e.g. Balmforth & Craster (Reference Balmforth and Craster1999) and Taylor-West & Hogg (Reference Taylor-West and Hogg2021) for more recent developments.

We then arrive at a coupled evolution system for the profile shape $h(x,t)$ of a foam-driven fracture and the pressure field $p(x,t)$ within the foam film, and at leading order this is given by

(2.7a) \begin{align}\frac {\partial }{\partial t} \big(\big[\rho _0 + c^{-2} (p-p_0)\big]h \big ) - \frac {r^2}{\gamma \mu \hat \kappa ^2} \frac {\partial }{\partial x}\left [\big[\rho _0 + c^{-2} (p-p_0)\big]\, h^3 \left (\frac {\partial p}{\partial x}\right )^2\right ] = 0, \end{align}

(2.7b) \begin{align}p(x,t) = - \frac {E}{2(1-\sigma ^2)}\frac {1}{\pi }\int _0^{x_f(t)} \frac {2s}{s^2-x^2}\frac {\partial h}{\partial s}\,\textrm{d}s. \end{align}

Compared with hydraulic fracturing with incompressible Newtonian fluids (e.g. Spence & Sharp Reference Spence and Sharp1985), the Reynolds equation (2.7a ) has been modified accordingly to absorb the influence of foam compressibility based on (2.4) and foam rheology based on (2.6). Equation (2.7b ) for the elastic pressure remains unchanged.

The governing system (2.7a , b ) is ready to be solved, providing the following appropriate initial and boundary conditions (IBC):

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

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

(2.8c) \begin{align}\frac {r^2}{\gamma \mu \kappa ^2} \big[\rho _0 + c^{-2} (p-p_0)\big] h^3 \left (\frac {\partial p}{\partial x}\right )^2\bigg |_{x=0} &= q(t), \end{align}

where $q(t)$ is the mass rate of the foam flow that enters the first quadrant of the domain (per unit depth perpendicular to the plane of paper). Physically, the initial condition (2.8 a) indicates that there is no fracturing fluid within the elastic solid before fracturing occurs, and the net pressure is zero. Equation (2.8 b) is a frontal condition, which corresponds to the location of the fracture’s tip. Equation (2.8 c) is a flux condition, which describes how fast fluid accumulates inside the fracture. It is also assumed that the half-width $a$ of the injection tube is much smaller than the half-length $x_f(t)$ of the fracture, such that the tube can be treated as located at approximately $x=0$ . To determine the influx $q(t)$ , we next consider the plug flow of the foam inside the straight tube, through which foam is delivered into the elastic solid to create a fracture.

2.2. Determining the influx $q(t)$

The compressible plug flow of foam in the injection tube is driven by the pressure difference $\Delta p \equiv p_i(t) - p(0,t)$ between the inlet ( $y=0$ ) and outlet ( $y=L$ ) of the tube. Here, $p_i(t)$ denotes the net pressure imposed at the inlet of the injection tube, while $p(0,t)$ denotes the net pressure at the centre of the fracture (which is also the outlet of the injection tube). Based on (2.5), we obtain $q(t) = \rho u a/2$ into the first quadrant without considering any transient dynamics, and this is a quarter of the total flux $2 \rho u a$ into the fracture. Furthermore, we obtain

(2.9) \begin{equation} q(t) = \frac {r^2 a^3}{2\gamma \mu \hat \kappa ^2}\, \big[\rho _0 + c^{-2} (p-p_0)\big] \left (\frac {\partial p}{\partial y}\right )^2, \end{equation}

the right-hand-side of which must be $y$ -independent for any quasi-steady operations.

Solution (2.9) can be further rearranged to provide the (nonlinear) pressure distribution along the tube:

(2.10) \begin{equation} p(y,t) = \left [p_i(t)^{3/2} - \left (\frac {9 \gamma \mu \hat \kappa ^2 c^2\, q(t)}{2 r^2 a^3 } \right )^{1/2} y \right ]^{2/3} - \big(\rho _0 c^2 - p_0\big). \end{equation}

The influx $q(t)$ can then be determined:

(2.11) \begin{equation} q(t) = \frac {9 r^2a^3 c^4}{8 \gamma \mu \hat \kappa ^2 L^2} \left ( \big[\rho _0 + c^{-2} (p_i(t) - p_0)\big]^{3/2} - \big[\rho _0 + c^{-2} (p(0,t) - p_0)\big]^{3/2} \right )^2, \end{equation}

which applies when $\Delta p \equiv p_i(t) - p(0,t) \geqslant 0$ .

It is easily seen that (2.10) and (2.11) are both different from the corresponding results when the fracturing fluid is an incompressible Newtonian fluid, in which case $p(y,t)$ depends linearly on $y$ (e.g. Zheng Reference Zheng2025a ,Reference Zheng b ). Another useful insight from solution (2.11) is that to enforce a constant injection rate $q$ into the fracture, the inlet pressure $p_i(t)$ must adjust accordingly with the outlet pressure $p(0,t)$ , the value of which must be determined by solving the full nonlinear non-local evolution system. In this work, we consider only the constant-pressure operations when the inlet pressure $p_i$ is an imposed constant.

2.3. Non-dimensionalisation

We have now formulated a complete problem of foam-driven fracturing with a set of governing equations (2.7a , b ) for the coupled evolution of the profile shape $h(x,t)$ and pressure field $p(x,t)$ . Equations (2.7a , b ) are to be solved subject to appropriate IBC (2.8a – c ). The influx $q(t)$ in (2.8 c) is given further by (2.11). The evolution is driven by an externally imposed inlet pressure $p_i(t)$ at the entrance of the injection tube, and in this work, $p_i$ is an imposed constant. The density of foam is denoted as $\rho _i$ at this (net) pressure. In the equation of state (2.4), we now impose $(p_0,\rho _0) = (p_i,\rho _i)$ as the reference (net) pressure and density, which allows us to obtain the density values subject to foam compression/expansion.

It is standard first to non-dimensionalise the model problem before solving it analytically or numerically. The dimensionless variables are defined as $\bar t \equiv t/t_c$ , $\bar \rho \equiv \rho /\rho _c$ , $\bar p \equiv p/p_c$ , $\bar x \equiv x/x_c$ and $\bar h \equiv h/h_c$ , and the characteristic scales are chosen as

(2.12a) \begin{align}t_c &=\frac {9E^5\gamma \mu \hat \kappa ^2 \rho _i^2 c^4 a^3}{256 \pi ^5(1-\sigma ^2)^5 r^2 L^2 p_i^9},\quad \rho _c = \rho _i, \ p_c = p_i, \end{align}

(2.12b) \begin{align}x_c &= \frac {9E^3 \rho _i^2 c^4 a^3}{64 \pi ^3(1-\sigma ^2)^3 L^2 p_i^5},\quad h_c = \frac {9E^2 \rho _i^2 c^4 a^3}{32 \pi ^2(1-\sigma ^2)^2 L^2 p_i^4}. \end{align}

The characteristic pressure and density scales are chosen as $(p_i, \rho _i)$ based on the injection condition. The characteristic time, length and thickness scales are chosen based on the Reynolds equation for the foam flow, the elastic equation for pressure, and the influx boundary condition for foam accumulation, such that the coefficients become 1 in their dimensionless form. Physically, this suggests that at dimensionless time $\bar t = 1$ , all physical effects come into play at leading order, including both the complex flow and the elastic deformation.

We then arrive at the dimensionless system

(2.13a) \begin{align}\frac {\partial }{\partial \bar t}\left [(1+\varPi _c(\bar p - 1))\bar h\right ] - \frac {\partial }{\partial \bar x}\left [(1+\varPi _c(\bar p - 1))\bar h^3 \left (\frac {\partial \bar p}{\partial \bar x}\right )^2\right ] = 0, \end{align}

(2.13b) \begin{align}\bar p = - \int _0^{\bar x_f(\bar t)} \frac {2s}{s^2 - \bar x^2} \frac {\partial \bar h}{\partial s}\, \textrm{d}s, \end{align}

subject to dimensionless IBC

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

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

(2.14c) \begin{align}(1+\varPi _c(\bar p - 1))\bar h^3 \left (\frac {\partial \bar p}{\partial \bar x}\right )^2 \bigg |_{\bar x=0} &= \left (1 - [1+\varPi _c(\bar p (0,\bar t) - 1)]^{3/2}\right )^2.\end{align}

We are then ready to solve (2.13a , b ) subject to IBC (2.14a - c ) for the dynamics of foam-driven fracturing under constant pressure operations.

We have also identified a dimensionless parameter

(2.15) \begin{equation} \varPi _c \equiv \frac {p_i}{\rho _i c^2}. \end{equation}

Physically, $\varPi _c^{1/2}$ is the Mach number at speed $(p_i/\rho _i)^{1/2}$ , hence $\varPi _c$ measures the compressibility of the plug flow of foam. We only consider the subsonic regime $\varPi _c \in (0,1)$ of foam compression/expansion in the present work. A greater $\varPi _c$ corresponds to more significant compressible effects in both the injection tube and the fracture, while $\varPi _c \to 0^+$ corresponds to the incompressible limit. The value of $\varPi _c$ distinguishes the ‘compressible’ and ‘incompressible’ regimes of foam-driven fracturing.

Figure 2.

Representative numerical solutions of foam-driven fracturing at dimensionless times $\bar t = 400,450,500$ , when we impose $\varPi _c = 0.001$ . The $\bar p(\bar x, \bar t)$ values ahead of the front represent the normal elastic stress along the axis, while the $\bar \rho (\bar x, \bar t)$ values ahead of the front contain no physical meaning because there exists no fluid. We have imposed $N=200$ space grids on a domain with $\bar L_d = 0.2$ in this simulation; see § A.1. In the bulk region there is reasonable convergence, even though in a narrow region close to the fracture’s tip, the pressure singularity is still under the influence of grid refinement. The density solutions ahead of the fracture’s tip have no physical meaning and are marked with a different colour.

Figure 3.

Representative numerical solutions of foam-driven fracturing at dimensionless times $\bar t = 400,450,500$ , when we impose $\varPi _c = 0.01$ . We have imposed $N=200$ space grids on a domain with $\bar L_d = 1.0$ in this simulation; see § A.1. In the bulk region, there is reasonable convergence, even though in a narrow region close to the fracture’s tip, the pressure singularity is still under the influence of grid refinement. The density solutions ahead of the fracture’s tip have no physical meaning and are marked with a different colour.

We can also focus on the fracture itself and evaluate its global mass/volume evolution:

(2.16) \begin{equation} \bar M(\bar t) \equiv \int _0^{\bar x_f(\bar t)} [1+\varPi _c(\bar p - 1)] \bar h(\bar x, \bar t)\, \textrm{d}\bar x,\,\,\, \bar V(\bar t) \equiv \int _0^{\bar x_f(\bar t)} \bar h(\bar x, \bar t)\, \textrm{d}\bar x, \end{equation}

where $\bar M(\bar t)$ is the dimensionless mass of the foam inside a fracture, while $\bar V(\bar t)$ is the dimensionless volume. Their evolution depends on many physical aspects of the model problem, including the plug flow of the foam, compressibility of the foam flow, elastic deformation of the infinite solid, and the dynamic coupling of all of them.

2.4. Numerical observations

We first make some representative calculations with different values of $\varPi _c$ to illustrate the different dynamics of foam-driven fracturing, as shown in figures 2 and 3. Key insights obtained are as follows.

1. (i) Influence of foam compressibility. It is observed in the numerical simulations that the fluid pressure decreases from the origin towards the fracture’s tip location, and so as the density of the foam (due to the expansion of the gas phase). An increase in $\varPi _c$ leads to more serious decrease of the fluid pressure and the foam density. Correspondingly, this leads to a more significant increase in fluid volume, and hence a more significant enhancement of film spreading and fracture propagation, as shown in figures 2 and 3. Furthermore, the volume accumulation departs from the linear law in time due to foam expansion, as shown in figure 4, and the departure becomes more significant as $\varPi _c$ increases.

2. (ii) Influence of foam rheology. The plug flow nature of the foam flow leads to modified profile shapes and pressure distributions, compared with standard results of hydraulic fracturing using incompressible Newtonian fluids (e.g. Spence & Sharp Reference Spence and Sharp1985). In particular, we have observed a linear asymptotic shape towards the fracture’s tip in the weakly compressible regime in the numerical simulations, as shown in figures 2 and 3. This is to be contrasted with the well-known $h \propto (x_f-x)^{2/3}$ viscous asymptote when the fracturing fluid is Newtonian. Later, we prove that the tip asymptotes follow $h \propto (x_f-x)$ and $p \propto \mbox{log}\,(x_f-x)$ as $x \to x_f^-(t)$ as $\varPi _c \to 0^+$ in the weakly compressible regime.

Figure 4.

The influence of $\varPi _c$ (which measures foam compressibility) on the time evolution of the total volume $\bar V(\bar t)$ of a foam-driven fracture. Here, we include the appropriately rescaled numerical data of the fracture’s volume for the period $\bar t \in [400,500]$ for three values of $\varPi _c$ (i.e. $\varPi _c = 0.001,0.01,0.02$ ). The solutions $\hat V = \hat t$ and $\hat V = \hat t^{6/5}$ are also included for comparison: an increase in $\varPi _c$ leads to more significant influence of foam compressibility, as it departs more significantly from the linear accumulation law. Specifically, power-law fitting of these numerical data provides $\hat V \sim 1.00\, \hat t^{1.08}$ when $\varPi _c=0.001$ , $\hat V \sim 1.00\, \hat t^{1.15}$ when $\varPi _c=0.01$ , and $\hat V \sim 1.00\, \hat t^{1.18}$ when $\varPi _c=0.02$ . The $R^2$ values are 1.00 for all three fittings.

3.1. Asymptotic structure near the tip of the fracture

It is of use to have some insights on the asymptotic structure near the fracture’s tip for both the profile shape $\bar h(\bar x, \bar t)$ and the pressure distribution $\bar p(\bar x, \bar t)$ . Unfortunately, there seems to be no self-similar behaviour for all time, since the injection rate (3.3) depends on $(\bar p(0,\bar t) - 1)^2$ . However, the tip shape can still be studied based on the assumption of quasi-steady evolution (e.g. Lister, Peng & Neufeld Reference Lister, Peng and Neufeld2013; Hewitt, Balmforth & Bruyn Reference Hewitt, Balmforth and Bruyn2015), which must be verified, however, based on appropriate comparisons with numerical solutions of the full time-dependent problem.

We define $\chi \equiv \bar x_f(\bar t) - \bar x$ as the distance from the tip location $\bar x_f(\bar t)$ , and we now evaluate $\bar h(\chi )$ and $\bar p(\chi )$ as $\chi \to 0^+$ at any specific time $\bar t$ of interest. Equations (3.2a , b ) further reduce to

(3.4a) \begin{align}c \bar h - \bar h^3 \left (\frac {\textrm{d} \bar p}{\textrm{d} \chi }\right )^2= 0, \end{align}

(3.4b) \begin{align}\bar p = - \int _0^{\infty } \frac {2s}{s^2 - \chi ^2} \frac {\textrm{d} \bar h}{\textrm{d} s}\, \textrm{d}s, \end{align}

where $c \equiv \textrm{d}\bar x_f/\textrm{d}\bar t$ is the speed of advancement of the fracture’s tip.

We now impose the power-law form of the profile shape near the tip,

(3.5) \begin{equation} \bar h(\chi ) \sim A \chi ^{\alpha }, \end{equation}

as $\chi \to 0^+$ , and the integral (3.4 b) for elastic pressure leads to (e.g. Spence & Sharp Reference Spence and Sharp1985; Peirce Reference Peirce2025)

(3.6) \begin{equation} \bar p(\chi ) \sim \left \{\begin{array}{ll} \displaystyle A \alpha \cot(\pi \alpha )\, \chi ^{\alpha -1}, & \alpha \in (1/2,1), \\ A \pi ^{-1} \log \chi , & \alpha =1. \end{array} \right .\\[3pt] \end{equation}

Substituting (3.5) and (3.6) into the Reynolds equation  (3.4 a), now considering only the influence of plug flow of foam in the weakly compressible limit, the correct balance leads to $\alpha = 1$ . Hence we must admit

(3.7) \begin{equation} \bar h(\chi ) \sim A \chi , \quad\bar p(\chi ) \sim A \pi ^{-1} \log \chi , \end{equation}

where $A = c^{1/4} \pi ^{1/2}$ .

This $\bar h(\chi ) \propto \chi$ tip asymptote in the viscosity-dominated regime is fundamentally different from the well-known $\bar h(\chi ) \propto \chi ^{2/3}$ tip asymptote of hydraulic fracturing using incompressible Newtonian fluids (e.g. Spence & Sharp Reference Spence and Sharp1985; Detournay Reference Detournay2004). This asymptote is also different from the $\bar h(\chi ) \propto \chi ^{1/2}$ tip asymptote in the material toughness-dominated regime. Its appearance is not due to foam compressibility, but due to the plug flow nature of foam rheology (2.6). The result (3.7) is verified based on a comparison with the numerical solutions of the full time-dependent problem, as shown in figure 5, normalised by the maximum fracture opening in each case. We show the range $\chi \in [0.1,1]$ here for $\varPi _c = 0.001,0.01$ , where the numerical solutions are more accurate (sufficiently far away from the tip region with pressure singularity, which requires more space grids to clearly capture numerically).

Figure 5.

The profile shape, normalised by its maximum value, near the fracture’s tip in the weakly compressible regime (as $\varPi _c \to 0^+$ ) from the numerical solutions of the full time-dependent problem. Numerical data are taken from figure 2(a) for $\varPi _c=0.001$ and figure 3(a) for $\varPi _c=0.01$ . The theoretical prediction (3.7) (slope 1) is also included for comparison. The slopes $1/2$ (toughness-dominated regime) and $2/3$ (viscosity-dominated regime) are also included, as the standard results of hydraulic fracturing using incompressible Newtonian fluids. The inset indicates that $\bar h/\chi$ indeed approaches a constant moving close to the tip region.

It may be of interest to note that this linear asymptotic profile $\bar h(\chi ) \propto \chi$ has also been obtained in a recent study on hydraulic fracturing using Herschel–Bulkley fluids (Bessmertnykh & Dontsov Reference Bessmertnykh and Dontsov2019), by setting $n=0$ in the imposed viscoplastic constitutive law $\tau = \tau _0 + k \dot \gamma ^n$ for perfectly plastic fluids. Note also that the model is formulated in a moving coordinate system along with the fracture’s tip, which is different from the present work (except in § 3.1 for the tip asymptotes), which maintains the time evolution in the Reynolds equation. Furthermore, the linear asymptotic profile $\bar h(\chi ) \propto \chi$ is also achievable using an incompressible power-law fluid ( $\tau = k \dot \gamma ^m$ ) for hydraulic fracturing, imposing $m=0$ (e.g. Dontsov & Kresse Reference Dontsov and Kresse2018).

Furthermore, $\bar p(\chi ) \propto \log \chi$ is weaker than $\bar p(\chi ) \propto \chi ^{-1/3}$ for tip asymptotes in the viscosity-dominated regime of hydraulic fracturing with Newtonian fluids. It is also weaker than the $\bar p(\chi ) \propto \chi ^{-1/2}$ tip asymptote in the toughness-dominated regime. The tip asymptote for pressure is more difficult to capture numerically due to the existence of the pressure singularity, which is well known. With some $N=200$ space grids imposed in the present work, we found that the $\bar p(\chi ) \propto \log \chi$ singularity is not clearly captured in the numerical solutions.

3.2. Self-similar evolution when $\bar p(0,\bar t)$ remains approximately constant

Because the injection rate (3.3) depends on $(\bar p(0,\bar t) - 1)^2$ , there seems to be no self-similar structure throughout the entire period of time evolution. However, if pressure $\bar p(0,\bar t)$ remains approximately unchanged during some finite time period of interest, then we expect the evolution to follow approximately a self-similar fashion in the viscosity-dominated regime.

Standard scaling analysis of (3.2a , b ) and (3.3), under the assumption that $\bar p(0,\bar t)$ is constant, indicates that

(3.8) \begin{equation} \frac {\bar x}{\bar t} \propto \frac {\bar h^3 \bar p^2}{\bar x^3},\quad \bar p \propto \frac {\bar h}{\bar x},\quad \bar h \bar x \propto \bar t, \end{equation}

and after some rearrangement, we arrive at

(3.9) \begin{equation} \bar x \propto \bar t^{5/9},\quad \bar h \propto \bar t^{4/9},\quad \bar p \propto \bar t^{-1/9}, \end{equation}

which shows that pressure $\bar p(\bar x, \bar t)$ evolves really with time. However, the evolution is quite slow and follows a $\bar t^{-1/9}$ time dependence, which suggests that the assumption that $\bar p(0,\bar t)$ remains approximately unchanged could be a reasonable assumption for quite a period of time. For example, when the time evolution covers one order of magnitude, the pressure scale becomes $10^{-1/9} \approx 0.77$ times the initial value. When the time evolution covers two orders of magnitude, the pressure scale becomes $100^{-1/9} \approx 0.60$ times the initial value. When the time evolution covers three orders of magnitude, the pressure scale becomes $1000^{-1/9} \approx 0.46$ times the initial value. All these calculations suggest that the assumption of neglecting the time evolution of pressure $\bar p(0,\bar t)$ could be a reasonable one for time periods that cover up to several orders of magnitude. This is verified based on the numerical solutions of the full time-dependent system, as we show later.

We next define a similarity variable as $\xi \equiv \bar x/\bar t^{5/9}$ , and explore self-similar solutions for the evolution of the profile shape $\bar h(\bar x, \bar t)$ and pressure field $\bar p(\bar x, \bar t)$ in the forms $\bar h(\bar x, \bar t) = \bar t^{4/9}\,\mathcal{K}(\xi )$ and $\bar p(\bar x, \bar t) = \bar t^{-1/9}\,\mathcal{G}(\xi )$ , respectively. We obtain

(3.10a) \begin{align}\frac {4}{9}\mathcal{K} -\frac {5}{9}\xi \frac {\textrm{d}\mathcal{K}}{\textrm{d}\xi } - \frac {\textrm{d}}{\textrm{d} \xi }\left [\mathcal{K}^3 \left (\frac {\textrm{d} \mathcal{G}}{\textrm{d} \xi }\right )^2\right ] = 0, \end{align}

(3.10b) \begin{align}\mathcal{G}(\xi ) = - \int _0^{\xi _f} \frac {2s}{s^2 - \xi ^2} \frac {\textrm{d} \mathcal{K}}{\textrm{d} s}\, \textrm{d}s, \end{align}

where $\xi _f \equiv \bar x_f(\bar t)/ \bar t^{5/9}$ is the frontal location. The ordinary differential system (3.10a , b ) is solved subject to boundary conditions

(3.11a) \begin{align}-\mathcal{K}^3 \left (\frac {\textrm{d} \mathcal{G}}{\textrm{d} \xi }\right )^2 \bigg |_{\xi =0} = \frac {9}{4} \varPi _c^{2} \left [\bar p(0,\bar t) - 1 \right ]^2, \end{align}

(3.11b) \begin{align}\mathcal{K}(\xi _f) = 0. \end{align}

Now the inlet pressure $\bar p(0,\bar t)$ is an imposed constant without any time evolution. This is an approximation, but it could be quite a reasonable one when the time evolution spans up to several orders of magnitude.

The self-similar system (3.10a , b ) can be solved numerically with (3.11) using the same numerical scheme that we have used in several previous studies (e.g. Zheng Reference Zheng2024, Reference Zheng2025a , Reference Zhengb ). The central idea is to maintain the time dependence in the similarity transform, defining $\tau \equiv -\ln \bar t$ and $\xi \equiv \bar x/\bar t^{5/9}$ , and transform the original system (3.2a , b ) in the $(\bar x, \bar t)$ space into the $(\tau , \xi )$ space. We then numerically solve the system in the $(\tau , \xi )$ space until there is no time evolution, and the steady-state solution provides the self-similar solution of interest. Brief remarks are also provided in § A.3.

Figure 6.

Comparison between the normalised numerical solutions $\bar h(\bar x, \bar t)$ and $\bar p(\bar x, \bar t)$ of the original time-dependent system (2.13a , b ), and the self-similar solutions $\mathcal{K}(\xi )$ and $\mathcal{G}(\xi )$ from solving (3.10a , b ). The numerical solutions are the appropriated rescaled ones in figure 2 at $\bar t = 400,450,500$ when $\varPi _c = 0.001$ . The self-similar solutions from solving (3.10a , b ) are included as the dashed curves for comparison. The self-similar solutions of hydraulic fracturing using incompressible Newtonian fluids are included also as square symbols (e.g. Spence & Sharp Reference Spence and Sharp1985; Detournay Reference Detournay2004). In the flux condition (3.11 a), the pressure at the origin is imposed as a constant $\bar p(0,\bar t) = 0.270$ for this time period based on the numerical data in figure 7(d).

We then provide a comparison in figure 6 between the normalised numerical solutions $\bar h(\bar x, \bar t)$ and $\bar p(\bar x, \bar t)$ of the partial differential integral system (2.13a , b ) and the self-similar solutions $\mathcal{K}(\xi )$ and $\mathcal{G}(\xi )$ of the ordinary differential integral system (3.10a , b ), when the pressure $\bar p(0,\bar t)$ remains approximately constant. Here, we have imposed $\varPi _c = 0.001$ in the weakly compressible and viscosity-dominated regime. The numerical solutions of (3.2a , b ) are provided at $\bar t = 400,450,500$ , and $\bar p(0,\bar t) = 0.27$ is taken as its initial value at $\bar t = 400$ in (3.11 a). We have observed reasonably good agreement between the normalised time-dependent solutions and the self-similar solutions in this regime. For comparison, the well-known self-similar solutions of hydraulic fracturing using incompressible Newtonian fluids are included also as the square symbols in the same figure (e.g. Spence & Sharp Reference Spence and Sharp1985; Detournay Reference Detournay2004).

Figure 7.

The propagation laws in the weakly compressible and viscosity-dominated regime of foam-driven fracturing. Here, we have provided numerical solutions of the frontal location $\bar x_f(\bar t)$ , inlet thickness $\bar h_m(\bar t) = \bar h(0,\bar t)$ , maximum pressure $\bar p_m(\bar t) = \bar p(0,\bar t)$ , and maximum density $\bar \rho _m(\bar t) = \bar \rho (0,\bar t)$ . The numerical solutions are obtained by solving the original system (2.13a , b ) for the period $\bar t \in [400,500]$ when $\varPi _c = 0.001$ . The scaling results (3.9a – c ) are also included for comparison. The insets of (a–c) show that the appropriately rescaled values by their corresponding power-law time dependences become time-independent.

Furthermore, the corresponding frontal laws are shown in figure 7, imposing $\varPi _c=0.001$ , along with a comparison with the theoretical prediction (3.9) in the viscosity-dominated and weakly compressible regime of foam-driven fracturing. We have also observed fairly good agreement between the theoretical and numerical results.

3.3. Remarks on the experimental observations of Lai et al. (2018)

One useful note is that the tip asymptote (3.7) as obtained in the Cartesian configuration should also apply in the radial configuration. Hence it could be used to be compared with experimental observations of penny-shaped foam-driven fractures as well. As an example, we deal with the final-time profile shape at $t = 173$ s in the experiment in figure 1, as performed by Lai et al. (Reference Lai, Rallabandi, Perazzo, Zheng, Smiddy and Stone2018). We then obtain a comparison with solution (3.7) in the $(\chi ,h/h_m)$ space for the fracture’s profile shape in the tip region in figure 8. Arguably, there is a reasonable degree of agreement between the theoretical prediction and the experimental observation. But the agreement is not as good as that in figure 5.

Figure 8.

The profile shape near the fracture’s tip at the final time ( $t = 173$ s) of the experiment shown in figure 1, as performed by Lai et al. (Reference Lai, Rallabandi, Perazzo, Zheng, Smiddy and Stone2018). The theoretical prediction (3.7) (slope 1) is also included for a comparison, now with $\chi \equiv (R_f-R)/R_f$ , where $R_f(t)$ denotes the frontal location of the penny-shaped fracture. The slopes $1/2$ (toughness-dominated regime) and $2/3$ (viscosity-dominated regime) are also included, as the standard results of hydraulic fracturing using incompressible Newtonian fluids.

Meanwhile, it is observed in the experiments of Lai et al. (Reference Lai, Rallabandi, Perazzo, Zheng, Smiddy and Stone2018) that the time evolution of the fracture’s volume could significantly depart from the $V \propto t$ linear law, using foam as the fracturing fluid. This was attributed to the influence of foam compressibility. This is consistent with the numerical observations in figure 4, even though foam is now injected under constant-pressure operations, which is different from how it is supplied in the experiments of Lai et al. (Reference Lai, Rallabandi, Perazzo, Zheng, Smiddy and Stone2018). Furthermore, we observe that an increase in the value of $\varPi _c$ leads to more significant departure from the linear accumulation law. All these observations illustrate the influence of foam compressibility.

It may also be of use to comment on one more interesting observation: when we stretch the self-similar solutions in figure 6 for the profile shape in both the $x$ and $h$ directions, we are able to obtain very promising agreement with the experimentally measured profiles, as shown in figure 9. This seems to suggest that for this experiment, the profile shape of a foam-driven fracture is predominantly affected by the specific foam rheology. Indeed, we have imposed the plug flow rheology in the model based on the experimental measurement of Lai et al. (Reference Lai, Rallabandi, Perazzo, Zheng, Smiddy and Stone2018). The profile shape depends only weakly on the radial geometry and other experimental features such as how the foam is supplied.

Figure 9.

The profile shapes in the experiment of Lai et al. (Reference Lai, Rallabandi, Perazzo, Zheng, Smiddy and Stone2018), along with a comparison with the numerical solutions in the viscosity-dominated and weakly compressible regime. Here, the theoretical profile shape is taken from the universal shape in figure 6 and stretched in both the $x$ and $h$ directions for the best match with the experimental profiles.

4.1. Brief summary

We have investigated the dynamics of foam-driven fracturing of an infinite elastic solid. Compared with the standard models of hydraulic fracturing using incompressible Newtonian fluids, we now take into consideration the influence of foam compressibility and foam rheology. We arrived at a modified Reynolds equation that describes the dynamics of foam-driven fracturing, which is coupled with the standard integral equation for elastic pressure. The influence of the plug flow of foam in the injection tube is integrated appropriately as a flux boundary condition.

Upon appropriate non-dimensionalisation, a dimensionless control parameter $\varPi _c$ has been identified, which affects the dynamics of foam-driven fracturing. Physically, $\varPi _c$ measures the relative influence of compressibility, and its value distinguishes the compressible and incompressible regimes. The dynamics of foam-driven fracturing are then demonstrated through numerical solutions of the nonlinear non-local system. Asymptotic analysis has also been attempted to obtain more insights into the weakly compressible regime.

Regarding the influence of the foam compressibility, because fluid pressure always decreases from the origin towards the tip location of a fracture, introducing compressibility always leads to an expansion of the gas phase and hence a decrease in the density of the foam. This feature enhances the spreading of the foam film and the propagation of the foam-driven fracture. Time evolution of the fracture’s total volume departs also from the linear accumulation law due to gas expansion, and this departure becomes more significant as $\varPi _c$ increases.

Regarding the influence of the foam rheology, we obtain a linear asymptotic profile shape $h \propto \chi$ as $\chi \to 0^+$ in the weakly compressible regime, which is to be compared with the well-known $h \propto \chi ^{2/3}$ asymptote for hydraulic fracturing with incompressible Newtonian fluids. This asymptote has been verified based on a comparison with the time-dependent numerical solutions. This asymptote is also, arguably, consistent with the experimental observations made by Lai et al. (Reference Lai, Rallabandi, Perazzo, Zheng, Smiddy and Stone2018), even though the fracture is penny-shaped in the experiments. Meanwhile, the pressure singularity becomes weaker as $p \propto \log\chi$ , to be compared with the standard result $p \propto \chi ^{-1/3}$ .

We have also observed that the inlet pressure $p(0,t)$ remains approximately constant in the weakly compressible regime during the time evolution for up to several decades. This indicates that the inlet flux of foam into the fracture remains constant. Correspondingly, there exist $x \propto t^{5/9}$ , $h \propto t^{4/9}$ and $p \propto t^{-1/9}$ (which is indeed very slow) scaling laws for the evolution of the length and pressure scales. The fracture’s profile shape and pressure distribution can also be rescaled to collapse onto universal solutions during the target period of time evolution.

We note also that we have neglected the influence of material toughness $K_c$ in this work, which typically requires integrating an additional crack propagation condition into the model problem (e.g. Lister Reference Lister1990b ; Savitski & Detournay Reference Savitski and Detournay2002; Zheng Reference Zheng2025a ). For hydraulic fracturing using Newtonian fluids, a fracture could evolve into an elliptic shape and propagate in a quasi-steady fashion, when the influence of material toughness is dominant. For foam-driven fracturing, its influence remains to be explored. Brief remarks are included in Appendix C.

4.2. Potential implications

Hydraulic fracturing has been known to consume a large amount of water in the context of energy resource recovery. Using foam as an alternative fracturing fluid has the potential to significantly reduce the consumption of water. In the future, it would be of interest to study foam-driven fracturing, e.g. in the context of shale gas recovery. This requires appropriate inputs of the foam rheology and compressibility at reservoir conditions. It would also be of interest to study the backflow dynamics following foam-driven fracturing, which is related to the risk of drinking water pollution.

Before we close this paper, we provide brief comments regarding the influence of foam compressibility at reservoir conditions. To fracture a typical shale rock of Young’s modulus $E = 15$ GPa, when a typical fracture of length $x_f = 1.0$ km and thickness $h_m = 10$ mm is generated, the elastic stress due to crack opening is estimated to be $p_f = E h_m/x_f \approx 0.15$ MPa. Presuming that the shale formation locates at some 1–3 km underneath the ground surface, a reasonable estimate of the external normal stress $\sigma _\infty$ (due to gravity) would be $\mathcal{O}$ (10) MPa. That said, the density change of foam along a fracture, which is related to the pressure change, is likely very small (say, 1 %, neglecting the pressure singularity in tip region). Hence the present model, imposing a linear law between $\rho$ and $p$ under the assumption of small pressure variations, is likely practically relevant for quantitative insights, and so are the example calculations in the weakly compressible regime. Also, the real impact of using foam as the fracturing fluid in the field is likely regarding its rheology (at high pressure and temperature, probably).

But the major influence of foam compressibility likely occurs during its flow in the injection well before a fracture is initiated, to our understanding. Its influence is also likely non-negligible during a short initial period of fast expansion right after a fracture is generated. These arguments are also based on typical experimental observations of foam/gas-driven fracturing in solid gelatines. The underlying physics, within the framework of linear elastic fracture mechanics, is that for a fracture to initiate, the stress intensity factor (which is $\propto p$ ) must be large enough to overcome the influence of material toughness $K_c$ (i.e. to break the atomic/molecular bonds), a typical value of which is $\mathcal{O}(10^2)$ Pa m $^{1/2}$ for a solid gelatine, and $\mathcal{O}(1)$ MPa m $^{1/2}$ for a shale rock. More discussions are included in Appendix C regarding the initial period of foam compression in the injection well, along with the influence of material toughness.