2. Model

We analyse the flow that arises when a buoyant fluid is injected into a porous layer that is bounded above by an impermeable seal; see figure 1. The flow is assumed to be axisymmetric, i.e. it is independent of the azimuthal angle. The radial coordinate is denoted by $R$, the vertical coordinate by $Z$, which is measured in the downwards direction in figure 1, and time by $T$. We assume that there is a sharp interface at $Z=H(R,T)$ between the input and ambient fluids, which we refer to as the ‘free surface’ (Lyle et al. Reference Lyle, Huppert, Hallworth, Bickle and Chadwick2005; Pegler et al. Reference Pegler, Huppert and Neufeld2014); see figure 1(b).

The flow is assumed to have a small aspect ratio (it is much thinner in the vertical direction than its extent in the radial direction) and so the radial velocity is much larger than the vertical velocity (Bear Reference Bear1988; Huppert & Woods Reference Huppert and Woods1995). The validity of this assumption for the present model is discussed a posteriori in §§ 3 and 4.

At later times, the pressure becomes approximately hydrostatic (Huppert & Woods Reference Huppert and Woods1995). The radial pressure gradient within the gravity current is then (Lyle et al. Reference Lyle, Huppert, Hallworth, Bickle and Chadwick2005)

(2.1)\begin{equation} {\frac{\partial P}{\partial R} = \rho g \frac{\partial H}{\partial R}}, \end{equation}

where $\rho =\rho _a-\rho _i>0$ is the density-difference between the ambient and input fluids.

For a uniform porous medium, Huppert & Pegler (Reference Huppert and Pegler2022) showed that the transition time to the gravity current regime, in which the flow is predominantly radial and the pressure hydrostatic, is given by

(2.2)\begin{equation} {T_{transition} \sim \varPhi \left( \frac{\mu^3 Q_{in}}{2 {\rm \pi}(\rho g K)^3} \right)^{1/2},} \end{equation}

where $\mu$ is the viscosity of the input fluid, $Q_{in}$ is the input flux, and $\varPhi$ and $K$ are the porosity and permeability, respectively. Our model requires that $T \gg T_{transition}$.

For a porous layer with non-uniform permeability, the Darcy velocity in the radial direction is obtained from (2.1) (Lyle et al. Reference Lyle, Huppert, Hallworth, Bickle and Chadwick2005; Hinton & Woods Reference Hinton and Woods2018),

(2.3)\begin{equation} U_R ={-}\frac{\rho g K(R,Z,T)}{\mu} \frac{\partial H}{\partial R}, \end{equation}

where $K(R,Z,T)$ is the horizontal permeability, which can vary in space and time. The permeability is related to the porosity, $\varPhi$, via (1.5). Throughout this section, upper case letters represent unscaled quantities except for the constant parameters, $\mu$, $\rho$ and $g$.

Initially, the porous medium has uniform porosity, $\varPhi (R,Z,0)=\varPhi _0$ and uniform permeability $K(R,Z,0)=K_0$, related via (1.5), and the medium is filled with ambient fluid so $H(R,0) = 0$. We assume that the introduction of the input fluid provides the required nutrients to stimulate biofilm growth. The porosity decreases in time according to the law

(2.4)\begin{equation} \frac{\partial \varPhi}{\partial T}=\left\{ \begin{array}{@{}ll@{}} - B F(\varPhi) & \mbox{if } Z\leqslant H(R,T),\\ 0 & \mbox{otherwise}, \end{array} \right. \end{equation}

where $B$ is the initial rate of porosity reduction and the dimensionless function $F({\cdot })$ is an input to the model and the choice of $B$ means that $F(\varPhi _0)=1$. There is no microbial growth in the ambient fluid. As an example, (1.3) furnishes

(2.5a,b)\begin{equation} {B=\beta(\varPhi_0-\varPhi_{\infty}), \quad F(\varPhi)=\frac{\varPhi-\varPhi_{\infty}}{\varPhi_0-\varPhi_{\infty}},} \end{equation}

whilst (1.4) gives

(2.6a,b)\begin{equation} {B=\hat{\beta} \varPhi_0^{1+{1}/{\gamma}}, \quad F(\varPhi)=\left(\frac{\varPhi}{\varPhi_0} \right)^{1+{1}/{\gamma}}.} \end{equation}

The radial flux of fluid (the depth-integrated velocity) is

(2.7)\begin{equation} Q_R = \int_0^{H} U_R \, \mathrm{d} Z ={-}\frac{\rho g }{\mu} \frac{\partial H}{\partial R} \int_0^{H} K(R,Z,T) \, \mathrm{d} Z. \end{equation}

To obtain a governing equation for the evolution of the free surface, we consider continuity over a thin vertical slice of the flow (cf. Lyle et al. Reference Lyle, Huppert, Hallworth, Bickle and Chadwick2005; Hinton & Woods Reference Hinton and Woods2018),

(2.8)\begin{equation} \frac{\partial }{\partial T}\left( \int_0^H \varPhi(R,Z,T) \, \mathrm{d}Z \right) = \frac{\rho g}{\mu R} \frac{\partial}{\partial R} \left[ R \frac{\partial H}{\partial R} \int_0^{H} K(R,Z,T) \, \mathrm{d} Z \right]. \end{equation}

Fluid is injected at the origin with volume flux, $Q_{in}$, evenly distributed over a finite distance in the $Z$ direction. Once the flow becomes predominantly radial (at times given by (2.2)), this distance over which the flux $Q_{in}$ is distributed becomes unimportant (Huppert & Pegler Reference Huppert and Pegler2022). Integrating the radial velocity (2.3) over the depth, $Z \in (0,H)$, at the origin furnishes the following boundary condition:

(2.9)\begin{equation} 2 {\rm \pi}\lim_{R\to 0} \left[-\frac{\rho g}{\mu} R \frac{\partial H}{\partial R} \int_0^{H} K(R,Z,T) \, \mathrm{d} Z \right]= Q_{in}. \end{equation}

Since injection begins at $T=0$, global volume conservation of the input fluid is given by

(2.10)\begin{equation} 2 {\rm \pi}\int_0^{R_f} \int_0^H \varPhi(R,Z,T) \, \mathrm{d} Z \, R \, \mathrm{d} R = Q_{in} T, \end{equation}

where $R_f(T)$ is the frontal contact point with $H(R_f(T),T) =0$.

In the case where $F(\varPhi )=0$ (corresponding to axisymmetric flow in a uniform layer with no microbial growth), the governing equations reduce to those of Lyle et al. (Reference Lyle, Huppert, Hallworth, Bickle and Chadwick2005) and Huppert & Pegler (Reference Huppert and Pegler2022). These researchers verified the model with analogue laboratory experiments.

2.1. Non-dimensionalisation

To non-dimensionalise the system, we introduce characteristic scales for the dimensional variables. The time scale for microbial growth is given by

(2.11)\begin{equation} \mathcal{T} = \frac{\varPhi_0}{B}, \end{equation}

and the velocity scale associated with the buoyant slumping of the fluid is

(2.12)\begin{equation} \mathcal{U} = \frac{\rho g K_0}{\mu}. \end{equation}

If the flow has radial length scale $\mathcal {L}$ and thickness scale $\mathcal {H}$, then a balance in (2.8) and (2.10) gives (Lyle et al. Reference Lyle, Huppert, Hallworth, Bickle and Chadwick2005; Dudfield & Woods Reference Dudfield and Woods2012)

(2.13a,b)\begin{equation} \frac{\varPhi_0 \mathcal{H}}{\mathcal{T}} \sim \frac{\mathcal{U} \mathcal{H}^2}{\mathcal{L}^2}, \quad \varPhi_0 \mathcal{H} \mathcal{L}^2 \sim \frac{Q_{in}}{2 {\rm \pi}} \mathcal{T}. \end{equation}

We obtain the following length and thickness scales:

(2.14a,b)\begin{equation} \mathcal{H} = \left( \frac{Q_{in}}{2 {\rm \pi}\mathcal{U}} \right)^{{1}/{2}}, \quad \mathcal{L}=\left( \frac{Q_{in} \mathcal{U}}{2 {\rm \pi}B^2} \right)^{{1}/{4}}. \end{equation}

As the input fluid invades the porous layer, the porosity and permeability are reduced, and so these quantities are scaled relative to their uniform initial values,

(2.15a,b)\begin{equation} \phi = \frac{\varPhi}{\varPhi_0}, \quad k=\frac{K}{K_0} = \phi^3, \end{equation}

where we have used (1.5). The governing equations and boundary conditions can now be non-dimensionalised using the following relations:

(2.16a–c)\begin{equation} t = \frac{T}{\mathcal{T}}, \quad (z,h) = \frac{(Z,H)}{\mathcal{H}}, \quad r=\frac{R}{\mathcal{L}}, \end{equation}

where the lower-case letters represent dimensionless quantities. The rate of change of the porosity (2.4) is re-expressed as

(2.17)\begin{equation} \frac{\partial \phi}{\partial t}=\left\{ \begin{array}{@{}ll@{}} - f(\phi) & \mbox{if } z\leqslant h(r,t),\\ 0 & \mbox{otherwise}, \end{array} \right. \end{equation}

where $f(\phi )= F(\varPhi )$ and $f(1)=1$. In the region uninvaded by the input fluid ($z>h(r,t)$), $\phi =1$ and $k= 1$. The governing equation (2.8) becomes

(2.18)\begin{equation} \frac{\partial }{\partial t}\left( \int_0^{h} \phi(r,z,t) \, \mathrm{d}z \right) = \frac{1}{r} \frac{\partial}{\partial r} \left( r \frac{\partial h}{\partial r} \int_0^h k(r,z,t) \, \mathrm{d} z \right). \end{equation}

The initial conditions are $\phi (r,z,0) \equiv 1$, $k(r,z,0) \equiv 1$ and $h(r,0) = 0$. The dimensionless form of the boundary condition at $r=0$, (2.9), is

(2.19)\begin{equation} \lim_{r\to 0} \left[{-}r \frac{\partial h}{\partial r} \int_0^{h} k(r,z,t) \, \mathrm{d} z \right]= 1, \end{equation}

and global volume conservation (2.10) becomes

(2.20)\begin{equation} \int_0^{r_f} \int_0^{h} \phi(r,z,t) \, \mathrm{d}z \, r \, \mathrm{d} r = t, \end{equation}

where $r_f(t)=R_f(T)/\mathcal {L}$. The system consisting of (2.17)–(2.20) is parameter-free with the exception of the functional form of the porosity variation, $f(\phi )$.

The evolution of the interface shape, $h(r,t)$, is obtained by solving (2.18) with initial condition $h(r,0) =0$, boundary condition (2.19), and the quantities $\phi (r,z,t)$ and $k(r,z,t)$ are obtained as follows. Initially, $\phi (r,z,0)=1$ everywhere and it then evolves according to (2.17), and the permeability is given by $k=\phi ^3$ (2.15b). The numerical method for approximating the solution, $h(r,t)$, is given in Appendix A (see also § 2.3).

2.2. Functional form of the porosity variation, $f(\phi )$

As discussed in § 1.2, we restrict our attention to two simple functions, $f(\phi )$, for how the rate of porosity reduction within the gravity current depends on the porosity (see (2.5a,b), (2.6a,b), (2.17)),

(2.21a,b)\begin{equation} f(\phi) = \frac{\phi - \phi_{\infty}}{1-\phi_{\infty}}, \quad f(\phi) = \phi^{1+{1}/{\gamma}}, \end{equation}

where $0 \leqslant \phi _{\infty }<1$ and $\gamma >0$. These two functions are plotted in figures 2(a) and 2(b), respectively, for different values of $\phi _{\infty }$ and $\gamma$. For $\phi _{\infty } = 0$ and $\gamma = \infty$, both expressions in (2.21a,b) become identical.

Figure 2. Forms of the porosity variation; see (2.17). (a) Function $f(\phi )$ from (2.21a) for three values of the late time porosity, $\phi _{\infty }$. (b) Function $f(\phi )$ from (2.21b) for three values of the decay rate $\gamma$. (c,d) Corresponding dependence of the porosity on the residence time of the input fluid, $t_R$; see (2.22), (2.23).

In both cases, (2.17) may be integrated to obtain the porosity as a function of the residence time of the input fluid, $t_R=t_R(r,z,t)$,

(2.22)$$\begin{gather} \phi(t_R) =\phi_{\infty} +(1-\phi_{\infty}) \exp\left(\frac{-t_R}{1-\phi_{\infty}}\right), \end{gather}$$

(2.23)$$\begin{gather}\phi(t_R) = \left(1+\frac{t_R}{\gamma}\right)^{-\gamma}, \end{gather}$$

with $t_R(r,z,t) = t-t_P(r,z)$, where $t_P(r,z)$ is the time at which the free surface passed the location $(r,z)$, and $t_R=0$ on the free surface. Provided that $\phi _{\infty }>0$, (2.22) corresponds to self-limiting microbial growth with $\phi \to \phi _{\infty }$ as $t_R \to \infty$; see figure 2(c). Equation (2.23) corresponds to progressively increasing bioclogging at long residence times; see figure 2(d).

2.3. Numerical results

The governing equations (2.17)–(2.20) are integrated numerically using finite differences. To handle the evolving permeability and porosity, the numerical method keeps track of the past evolution of the free surface, $h(r,\tau )$ for $\tau \in [0,t]$. The history of the free surface determines the residence time of the input fluid, $t_R(r,z,t)$, which in turn furnishes the porosity, $\phi (r,z,t)=\phi (t_R)$, and permeability, $k(r,z,t)=k(t_R)$, within the flow. Full details of the numerical method are given in Appendix A. Throughout this paper, the solutions are shown with the $z$ axis in the upwards direction noting that a simple reflection in $z=0$ relates the figures to the geological motivation shown in figure 1.

Results for the evolution of the free surface, $h(r,t)$, are shown in figure 3 at $t=0.1,1,10$ with the porosity decreasing according to (2.22). The three panels in figure 3 correspond to different limiting values of the porosity at long residence times ($\phi _{\infty }=0.5,0.25,0$). At early times (e.g. $t=0.1$), the porosity variation has a relatively small influence on the flow and the shape of the free surface is similar across the three panels. At later times, the greater reduction in porosity and permeability near the origin for $\phi _{\infty }=0.25$ and $\phi _{\infty }=0$ (panels b,c) causes the flow to invade a much greater area of rock and the gravity current has an increased aspect ratio; these features are discussed in more detail in §§ 3 and 4.

Figure 3. Shape of the free surface $z=h(r,t)$ at $t=0.1, 1, 10$ for an exponentially decaying porosity within the input fluid (2.22) with late-time porosity given by (a) $\phi _\infty =0.5$, (b) $\phi _\infty =0.25$ and (c) $\phi _\infty =0$.

Figure 4 shows the residence time of the input fluid, $t_R$, the porosity, $\phi$, and the permeability, $k$, within the gravity current for the case shown in figure 3(a). The free-surface history determines the porosity and permeability, which in turn influence the flow structure and hence the future free-surface evolution. Figure 4 demonstrates that the gradual increase in residence time away from the free surface has a small influence on the porosity at early times but a large influence at late times, with the porosity eventually becoming approximately $\phi _{\infty }$ almost everywhere within the porous gravity current.

Figure 4. Properties of the porous rock at $t=0.1, 1, 10$ for an exponentially decaying porosity within the input fluid (2.22) with late-time porosity $\phi _\infty =0.5$; see also figure 3(a). Dashed lines show the free surface. (a–c) Residence time of the input fluid at $t=0.1, 1, 10$. (d–f) Corresponding porosity, $\phi$, which tends to $\phi _{\infty }=0.5$ at long residence times. (g–i) Corresponding permeability, $k= \phi ^3$, which tends to $k_{\infty }=0.125$ at long residence times.

3. Effect of limited microbial growth on the gravity current ($\phi \to \phi _{\infty } >0$)

In this section, the case in which the microbial growth becomes self-limiting with $\phi \to \phi _{\infty } >0$ as $t \to \infty$ is analysed; see figure 4. The case in which $\phi$ becomes progressively smaller at late times is qualitatively different and discussed in § 4.

3.1. Early time

At early times, $t \ll 1$, the porosity and permeability have changed little from their initial values (e.g. see figures 4a,4d,4g). Hence, we write

(3.1a,b)\begin{equation} \phi(r,z,t) = 1 + O (t), \quad k(r,z,t)=1+ O (t). \end{equation}

The leading-order terms in the governing equation (2.18) and global volume conservation (2.20) are

(3.2)\begin{equation} \frac{\partial h}{\partial t} = \frac{1}{r} \frac{\partial}{\partial r} \left(r h \frac{\partial h}{\partial r} \right) \end{equation}

and

(3.3)\begin{equation} \int_0^{r_f(t)} h r \, \mathrm{d} r = t, \end{equation}

respectively. The porous layer is effectively uniform and the evolution of the free surface is self-similar with solution (Lyle et al. Reference Lyle, Huppert, Hallworth, Bickle and Chadwick2005)

(3.4a,b)\begin{equation} h(r,t) = G\left(\eta\right), \quad \eta = \frac{r}{t^{1/2}}, \end{equation}

where $G(\eta )$ satisfies

(3.5a–c)\begin{equation} -\frac{1}{2} \eta \frac{\mathrm{d} G}{\mathrm{d} \eta} = \frac{1}{\eta} \frac{\mathrm{d}}{\mathrm{d} \eta} \left( \eta G \frac{\mathrm{d} G}{\mathrm{d} \eta} \right), \quad \int_0^{\eta_f} \eta G \, \mathrm{d} \eta =1, \quad G(\eta_f)=0, \end{equation}

with $\eta _f=r_f(t)/t^{1/2}$. The shape function, $G(\eta )$, is obtained via numerical integration (Lyle et al. Reference Lyle, Huppert, Hallworth, Bickle and Chadwick2005). This self-similar solution is shown with blue dots in figure 5 in the $r/t^{1/2}$ coordinates. Figure 5 also shows the numerical results for the free surface at five times (black lines) in the case where the porosity decays exponentially with $\phi \to \phi _{\infty }=0.5$. There is excellent agreement between the similarity solution and the full numerical results at $t=0.01$ and $t=0.1$.

Figure 5. Position of the free surface in rescaled coordinates, $r/t^{1/2}$, at $t=0.01, 0.1, 1, 10, 100$ for exponential porosity decay (2.22) with $\phi _\infty =0.5$. The blue dots show the early-time self-similar solution, (3.4a,b), and the red dots show the late-time self-similar solution, (3.15).

Contours of the residence time of the input fluid, $t_R(r,z,t)$, are given by the shape of the free surface at early times. Hence, the solution $G(\eta )$ furnishes the following implicit equation for $t_R$:

(3.6)\begin{equation} z = G \left(\frac{r}{\left(t-t_R\right)^{1/2}} \right), \end{equation}

which is valid for $(t-t_R) \ll 1$, i.e. the time of the free surface passing $(r,z)$ is small. Similar expressions can be obtained for contours of the porosity and permeability by rewriting (2.22) as $\phi =1-t_R$ for $t_R\ll 1$ and using (2.15b),

(3.7a,b)\begin{equation} z = G \left(\frac{r}{\left(t-1+\phi\right)^{1/2}} \right), \quad z = G \left(\frac{r}{\left(t-1+k^{1/3} \right)^{1/2}} \right). \end{equation}

The porosity contours (3.7a) are shown as red dashed lines in figure 6 for $t=0.05$ and they show good agreement with the numerical results (continuous lines); parameter values are as in figure 5.

Figure 6. Contours of the porosity, $\phi =0.96,0.97,0.98,0.99$, at $t=0.05$ for exponential porosity decay and $\phi _\infty =0.5$ (continuous lines). The red dashed lines show the predictions of (3.7a).

Next, we discuss the validity of the hydrostatic pressure assumption, which is associated with the radial extent of the flow being much smaller than the characteristic thickness. We denote by $\lambda$ the ratio of the characteristic thickness of the flow to its radial length scale. For the early-time behaviour, $r\sim t^{1/2}$ and $h\sim 1$, and so

(3.8)\begin{equation} \lambda \sim \frac{\mathcal{H}}{\mathcal{L} t^{1/2}}, \end{equation}

where $\mathcal {H}$ and $\mathcal {L}$ are given by (2.14a,b). We require $\lambda \ll 1$ and for the early-time approximations, we also require $t\ll 1$. Together, this gives the following condition for the present early-time analysis to apply:

(3.9)\begin{equation} B \left(\frac{Q_{in}}{2{\rm \pi}}\right)^{1/2} \mathcal{U}^{{-}3/2} \ll t \ll 1. \end{equation}

At these times, (2.2) is also satisfied.

Finally, we note that the shape function $G(\eta )$ has a logarithmic singularity at the origin, which violates the assumption of hydrostatic pressure, but it can be removed by incorporating the fully three-dimensional pressure-driven flow in this region. This is outside the scope of the present study and the interested reader is referred to Benham et al. (Reference Benham, Neufeld and Woods2022).

3.2. Long-time behaviour for $\phi \to \phi _{\infty }>0$

At late times, $t \gg 1$, the microbial growth has become limited ($\partial \phi /\partial t \approx 0$) within most of the pore space occupied by the gravity current. The porosity and permeability within the current are approximately uniform with

(3.10a,b)\begin{equation} \phi=\phi_{\infty}, \quad k= k_{\infty}, \end{equation}

as can be seen in figures 4(f) and 4(i). The leading-order terms in the governing equation (2.18) and global volume conservation (2.20) are

(3.11)$$\begin{gather} \phi_{\infty} \frac{\partial h}{\partial t} = k_{\infty} \frac{1}{r} \frac{\partial}{\partial r} \left(r h \frac{\partial h}{\partial r} \right), \end{gather}$$

(3.12)$$\begin{gather}\int_0^{r_f} \phi_{\infty} h r \, \mathrm{d} r = t, \end{gather}$$

where $\phi _{\infty }>0$. Under the rescaling

(3.13a,b)\begin{equation} h=k_{\infty}^{{-}1/2} \tilde{h}, \quad r= k_{\infty}^{1/4}\phi_{\infty}^{{-}1/2} \tilde{r}, \end{equation}

we recover the early-time equation (3.2) in the tilde variables. Hence, the solution is given by $\tilde {h}=G(\tilde {r}/t^{1/2})$, where $G(\eta )$ is as in § 3.1.

Equation (3.13a,b) demonstrates that the reduction in permeability by a factor of $k_\infty <1$ leads to a thicker current with lesser radial extent, whilst the reduction in porosity leads to an unchanged thickness and a greater radial extent. It should be noted that the change in thickness and radial extent given by (3.13a,b) is independent of the choice of relation between the porosity and permeability. Using the particular relation $k_\infty =\phi _\infty ^3$, the influence of the porosity and permeability change can be amalgamated so that the rescaling (3.13a,b) becomes

(3.14a,b)\begin{equation} h=\phi_{\infty}^{{-}3/2} \tilde{h}, \quad r= \phi_{\infty}^{1/4} \tilde{r}. \end{equation}

This demonstrates that the aggregate effect of the microbial growth is that the flow becomes substantially thicker and slightly shorter. We can write the self-similar solution as

(3.15)\begin{equation} h = \phi_{\infty}^{{-}3/2} G \left( \frac{r}{\phi_{\infty}^{1/4} t^{1/2}} \right), \end{equation}

which is shown as a red dotted line in figure 5 for $\phi _\infty =0.5$. There is excellent agreement with the numerical result at $t=100$. Approximate contours for the residence time, $t_R(r,z,t)$, can be obtained in an analogous fashion to (3.6), and contours for the porosity and permeability then follow using (2.22) and (2.15b).

The condition for the validity of the assumption of hydrostatic pressure (3.9) is re-expressed for the present late-time behaviour as

(3.16)\begin{equation} t \gg B \left(\frac{Q_{in}}{2{\rm \pi}}\right)^{1/2} \mathcal{U}^{{-}3/2} \phi_{\infty}^{{-}7/2}. \end{equation}

At these times, (2.2) is also satisfied.

4. Gradual clogging of the porous medium

We analyse the case in which the porous layer becomes gradually clogged by the biomass at long residence times, i.e. $\phi \ll 1$ and $k \ll 1$ for $t_R \gg 1$. We show that the late-time behaviour depends qualitatively on the rate at which $\phi$ decreases with two distinct regimes: relatively slow bioclogging and relatively fast bioclogging. Throughout this section, we assume that the porosity and permeability are relatively small but non-zero.

To expound the distinction between ‘slow’ and ‘fast’ bioclogging, we focus on algebraic decay in time of the porosity and permeability (this case is also convenient as it is associated with self-similar evolution of the flow). The evolution of the porosity as a function of the residence time (2.23) is $\phi (t_R)=(1+t_R/\gamma )^{-\gamma }$, which is shown graphically in figure 2(d) for various values of $\gamma$. The permeability is given by $k(t_R)=\phi (t_R)^3$. Exponential decay of the porosity (and permeability), $\phi =\exp (-t_R)$, is recovered in the limit $\gamma \to \infty$.

Numerical results for $\gamma =0.2$ and $t=1,10,100$ are shown in figure 7 illustrating the variation in the residence time of the input fluid and the rock properties. The minimum value of the porosity at $t=100$ is $0.288$ and the minimum value of the permeability is $0.024$, both of which are attained at the origin.

Figure 7. Properties of the porous rock at $t=1, 10, 100$ for an algebraically decaying porosity within the input fluid (2.23) with exponent $\gamma =0.2$. Dashed lines show the free surface. (a–c) Residence time of the input fluid, $t_R(r,z,t)$, at $t=1, 10, 100$. (d–f) Corresponding porosity, $\phi$. (g–i) Corresponding permeability, $k= \phi ^3$. The minimum value of $k$ in (i) is $0.024$.

The dimensionless flow speed in the radial direction within the input fluid, ${u=-k \partial h/\partial r}$, is shown in figure 8 at $t=10$ (corresponding to the permeability in figure 7h). In a uniform porous layer, the flow speed would be independent of $z$. Here, there is flow rerouting owing to the change in the permeability. However, not all the flow goes through the high-permeability zone near the free surface. Buoyancy also plays a role in driving the flow, especially near the source, where there is significant flux through the lower permeability regions (cf. Hinton & Woods Reference Hinton and Woods2018).

Figure 8. Dimensionless flow speed in the radial direction within the input fluid, $u=-k \partial h/\partial r$ at $t=10$ for $\gamma =0.2$ (corresponding to figures 7b,7e,7h).

We seek a solution to the governing equation (2.18) at late times, $t \gg 1$, that accounts for the decrease in porosity and permeability away from the free surface. First, we introduce a new variable, $\tau _R=t_R/t$, the relative residence time of the input fluid. The integral of the permeability, $k=(1+t_R/\gamma )^{-3\gamma }$, over the flow thickness is re-expressed as an integral over the past evolution of the free surface, $h(r,t-t_R)=h(r,t(1-\tau _R))$ with $\tau _R \in [0,\tau _{max}]$, where $\tau _{max} <1$ is the relative residence time at $z=0$, which depends on $r$ and $t$. The depth-integrated permeability in the governing equation (2.18) becomes

(4.1)\begin{align} \int_0^h k \, \mathrm{d} z &={-}\int_0^{\tau_{max}} \left(1+\frac{t \tau_R}{\gamma} \right)^{{-}3 \gamma} \frac{ \mathrm{d} z}{\mathrm{d} \tau_R} \, \mathrm{d} \tau_R \end{align}

(4.2)\begin{align} &=\left.\int_0^{\tau_{max}} \left(1+\frac{t \tau_R}{\gamma} \right)^{{-}3 \gamma} t \frac{ \partial h}{\partial t}\right\vert_{(r,t(1-\tau_R))} \, \mathrm{d} \tau_R. \end{align}

The depth-integrated porosity is given by an equivalent expression with the exponent $-3 \gamma$ replaced by $- \gamma$.

At late times, the depth-integrated permeability (4.2) can be approximated by considering the contributions from two different regions: (a) the global contribution; and (b) the contribution from a neighbourhood of the free surface where $\tau _R \ll 1$. We assume that $\tau _{max}$ is $ O (1)$, which is verified later. The contribution to the integral (4.2) from each region is given by the magnitude of the integrand multiplied by the width of the region (see § 3.4 of Hinch Reference Hinch1991),

(4.3)$$\begin{gather} \text{Global: if}\ \tau_R = O \left(1 \right), \quad \left(1+\frac{t \tau_R}{\gamma} \right)^{{-}3 \gamma}= O \left(t^{{-}3\gamma}\right), \quad \int_0^h k \, \mathrm{d} z = O \left(t^{1-3\gamma} \frac{ \partial h}{\partial t} \right). \end{gather}$$

(4.4)$$\begin{gather}\text{Local: if}\ \tau_R = O \left(t^{{-}1} \right), \quad \left(1+\frac{t \tau_R}{\gamma} \right)^{{-}3 \gamma}= O \left(1\right), \quad \int_0^h k \, \mathrm{d} z = O \left(\frac{ \partial h}{\partial t} \right). \end{gather}$$

For $0<\gamma <1/3$, the global contribution ($\tau _R = O (1)$) is dominant, whilst for $\gamma >1/3$, the local contribution near the free surface ($\tau _R = O (t^{-1})$) is dominant. A similar result applies to the depth-integrated porosity with the global contribution dominating when $0< \gamma <1$ and the local contribution dominating when $\gamma >1$.

This distinction divides the flow behaviour into two regimes, with ‘slow’ bioclogging in the case where $0< \gamma < 1/3$, for which both the depth-integrated porosity and permeability consist of contributions from the entire flow thickness. The self-similar behaviour for this regime is described in the present section. Figure 8 ($\gamma =0.2$) demonstrates that the depth integrated flux $q=\int _0^h u \, \mathrm {d}z$ will include contributions from across the flow thickness.

The case of ‘fast’ bioclogging is split into two subregimes: (i) for $1/3 <\gamma <1$, the depth-integrated permeability is dominated by the region near the free surface but the depth-integrated porosity is not; (ii) for $\gamma >1$, both quantities are dominated by the free-surface region. For these regimes, the porosity and permeability can become small very quickly and so the injection pressure may need to be very high to continue supplying fluid. Their analysis is included for completeness in Appendices B and C.

The shape of the free surface in the case of algebraic decay of the porosity is shown at $t=10$ in figure 9(a) and at $t=100$ in figure 9(b) for $\gamma =0.1, 0.2, 0.5, 1, 2, \infty$ (exponential decay). As expected, there is a qualitative difference in the evolution of the free surface for $0< \gamma < 1/3$ (the red and blue curves).

Figure 9. Shape of the free surface for complete bioclogging ($\phi \to 0$ as $t \to \infty$) at (a) $t=10$ and (b) $t=100$. The six curves correspond to the following exponents for the algebraic decay of the porosity: $\gamma =0.1, 0.2, 0.5, 1, 2$ and the case of exponential decay ($\gamma \to \infty$); see (2.23). There is a qualitative change in behaviour across $\gamma =1/3$; see discussion in the text.

In the present section, we derive the self-similar solution for the flow in the case of relatively slow bioclogging, $0<\gamma < 1/3$. The depth-integrated permeability is given by (4.2) with the approximation (4.3) furnishing the scaling for the dominant contribution, $\int _0^h k \, \mathrm {d} z \sim t^{-3\gamma } h$. The depth-integrated porosity scales with $t^{-\gamma } h$. The time exponents for the self-similar scaling of the radial coordinate $r \sim t^a$ and the flow thickness $h \sim t^b$ are obtained by balancing the terms in the governing equation (2.18) and global volume conservation (2.20). This furnishes the self-similar form,

(4.5a,b)\begin{equation} h(r,t)=t^{{3\gamma}/{2}} \gamma^{{-}2\gamma} C_0^2 \psi(\xi), \quad \xi = \frac{r}{C_0 t^{({2-\gamma})/{4}}}, \end{equation}

where the constant $C_0$ is defined so that $\xi = 1$ at the contact point, $r=r_f(t)$ (i.e. ${\psi (1) = 0}$). The shape function $\psi (\xi )$ and the constant $C_0$ are to be determined.

To obtain an integro-differential equation governing $\psi (\xi )$, the depth-integrated permeability (and porosity) must be expressed in terms of the self-similar coordinates, (4.5a,b),

(4.6)\begin{equation} \int_0^h k \, \mathrm{d} z = \int_1^{\xi} k \frac{\mathrm{d} z}{\mathrm{d} s} \, \mathrm{d} s ={-} \int_{\xi}^1 \left( \frac{ t \tau_R}{\gamma} \right)^{{-}3\gamma} \, \frac{\mathrm{d}}{\mathrm{d} s}\left[ (t(1-\tau_R))^{3\gamma/2} \gamma^{{-}2\gamma} C_0^2 \psi\left(s\right)\right] \mathrm{d} s, \end{equation}

where

(4.7)\begin{equation} s= \frac{r}{C_0 \left(t-t_R \right)^{({2-\gamma})/{4}} }= \frac{\xi}{\left(1-\tau_R\right)^{({2-\gamma})/{4}}} \end{equation}

encapsulates the ‘history’ of the free surface in the self-similar coordinates, and so $\tau _R$ and $s$ are interdependent. Note that $s=\xi$ at $z=h$ and $s = 1$ at $z=0$ because $\psi (1) =0$.

At a point $(\xi, \psi (\xi ))$ in similarity space, the prior evolution of the free surface, $h(r,t\tau _R)$ with $\tau _R \in [0,\tau _{max}]$, is encoded in $\psi (s)$ with $s \in [\xi,1]$; see figure 10.

Figure 10. Relationship between the past evolution of the free surface in real space $\tau _R \in [0,\tau _{max}]$ and the past evolution in similarity space $s \in [\xi,1]$; see (4.7).

Equation (4.7) can also be used to obtain the relative residence time, $\tau _R$, at $z=0$: $\tau _{max} = 1- \xi ^{{4}/({2-\gamma })}$, which is $ O (1)$ as assumed earlier. The integrand in (4.6) is re-written in terms of $s$ and $\xi$ only,

(4.8)\begin{equation} \int_0^h k\, \mathrm{d} z ={-} \gamma^{\gamma} t^{{-}3\gamma/2} C_0^2 \int_{\xi}^1 \left[1-\left(\frac{\xi}{s} \right)^{{4}/({2-\gamma})} \right]^{{-}3\gamma} \frac{\mathrm{d}}{\mathrm{d} s} \left[ \left(\frac{\xi}{s} \right)^{{6 \gamma}/({2-\gamma})} \psi(s) \right] \mathrm{d} s. \end{equation}

The depth-integrated porosity is obtained via an almost identical calculation,

(4.9)\begin{equation} \int_0^h \phi\, \mathrm{d} z ={-}\gamma^{-\gamma} t^{\gamma/2} C_0^2 \int_{\xi}^1 \left[1-\left(\frac{\xi}{s} \right)^{{4}/({2-\gamma})} \right]^{-\gamma} \frac{\mathrm{d}}{\mathrm{d} s} \left[ \left(\frac{\xi}{s} \right)^{{6 \gamma}/({2-\gamma})} \psi(s) \right] \mathrm{d} s. \end{equation}

Using these formulae and (4.5a,b), the governing equation (2.18) is recast as an integro-differential equation for the shape function $\psi (\xi )$,

(4.10)\begin{equation} \frac{\gamma}{2} J_1 - \left(\frac{1}{2} -\frac{\gamma}{4} \right) \xi \frac{\mathrm{d} J_1}{\mathrm{d} \xi} = \frac{1}{\xi} \frac{\mathrm{d}}{\mathrm{d}\xi} \left( \xi J_3 \frac{\mathrm{d} \psi}{\mathrm{d} \xi} \right), \end{equation}

where the functions $J_n(\xi )$ are associated with the integrals of the porosity ($n=1$) and permeability ($n=3$) over the thickness,

(4.11)\begin{equation} J_n(\xi) ={-} \int_{\xi}^1 \left( 1- \left( \frac{\xi}{s} \right)^{{4}/({2-\gamma})} \right)^{{-}n \gamma} \frac{\mathrm{d} }{\mathrm{d} s} \left[ \left( \frac{\xi}{s} \right)^{{6 \gamma}/({2-\gamma})} \psi(s) \right] \mathrm{d} s. \end{equation}

Although the integrand in $J_n(\xi )$ is singular as $s \to \xi$, this is integrable provided that $\gamma < 1/n$, which is satisfied for $J_1$ and $J_3$ in the ‘slow’ regime. Global mass conservation (2.20) becomes

(4.12)\begin{equation} C_0 = \left(\gamma^{-\gamma} \int_0^1 J_1(\xi) \xi \, \mathrm{d} \xi \right)^{{-}1/4}. \end{equation}

If we set $\gamma =0$, then $J_n(\xi ) = \psi (\xi )$ and the similarity solution of Lyle et al. (Reference Lyle, Huppert, Hallworth, Bickle and Chadwick2005) is recovered for flow in a uniform porous layer. For $0<\gamma <1/3$, (4.10) is solved by numerically integrating backwards from the contact point $\xi =1$; for details, see Appendix A.2. The similarity solution is compared with the numerical integration of the full governing system in figure 11 for $\gamma =0.1,0.2,0.3$, and there is good agreement at late times.

Figure 11. Self-similar evolution of the free surface in the case of ‘slow’ bioclogging. The numerical results are shown at $t= 1, 10, 100$ (black, green and blue lines) and compared with the similarity solution (4.5a,b) for (a) $\gamma =0.1$, (b) $\gamma =0.2$ and (c) $\gamma =0.3$ (red dashed lines).

The similarity solution (4.5a,b) also determines the spatial evolution of the residence time of the input fluid, $t_R(r,z,t)$. The porosity and permeability can then be inferred using $\phi = (t_R/\gamma )^{-\gamma }$ and $k= (t_R/\gamma )^{-3\gamma }$, which are valid at late times. For example, contours of the porosity are given by

(4.13)\begin{equation} z = \left(t-\gamma\phi^{-{1}/{\gamma}} \right)^{{3 \gamma}/{2}} \gamma^{{-}2 \gamma} C_0^2 \psi \left( \frac{r}{C_0 \left(t-\gamma\phi^{-{1}/{\gamma}} \right)^{{(2- \gamma)}/{4}}} \right). \end{equation}

The characteristic thickness of the gravity current increases in proportion to $t^{3\gamma /2}$ whereas its radial extent grows as $t^{(2-\gamma )/4}$ meaning that the aspect ratio $h/r$ is proportional to $t^{(7\gamma -2)/4}$. The present late-time analysis is thus valid provided that $\gamma <2/7$. For $2/7< \gamma <1/3$, the thickness of the gravity current grows faster than the radial extent and so the assumption that the pressure is hydrostatic is not valid.

5. Incorporating consumption of the input fluid

To account for the loss of input fluid owing to microbial activity, we assume that each unit volume of biomass gained required $\alpha \geqslant 0$ units (volume) of the input fluid to be consumed. The parameter $\alpha$ represents the aggregate effect of a range of reactions (see § 1.2).

The total increase in the volume of biomass since the start of injection is equal to the total reduction in the porosity, given by

(5.1)\begin{equation} {2 {\rm \pi}\int_0^{R_f} \int_0^H (\varPhi_0 -\varPhi) \,\mathrm{d}Z \, R \, \mathrm{d}R,} \end{equation}

where $\varPhi _0$ is the uniform initial porosity (see § 2). Thus, upon incorporating the consumption of the input fluid, dimensionless global volume conservation (2.20) becomes

(5.2)\begin{equation} {\int_0^{r_f} \int_0^h \phi + \alpha (1-\phi) \, \mathrm{d} z \, r \, \mathrm{d} r = t,} \end{equation}

where the first term in the integrand is associated with the input volume that is still mobile fluid and the second term accounts for the input volume that has become biomass. Similarly, the dimensionless governing equation (2.18) becomes

(5.3)\begin{equation} {\frac{\partial}{\partial t} \left( \int_0^h \phi + \alpha (1-\phi) \, \mathrm{d} z \right) = \frac{1}{r} \frac{\partial}{\partial r} \left(r \frac{\partial h}{\partial r} \int_0^h k \, \mathrm{d} z \right).} \end{equation}

The boundary conditions, initial conditions, and the laws for the evolution of the porosity and permeability are unchanged from § 2. The case $\alpha =0$ corresponds to no consumption (studied earlier) whilst $\alpha =1$ is associated with a balance whereby the volume loss of input fluid is matched by the increase in volume of the biomass.

5.1. Limited microbial growth ($\phi \to \phi _\infty >0$)

We first consider the effect of consumption on the regime in which the microbial growth becomes self-limiting ($\phi \to \phi _\infty >0$). The decay rate of the porosity is given by (2.21a) and the permeability is $k=\phi ^3$. The shape of the free surface at $t=1$ and $t=100$ is shown in figures 12(a) and 12(b), respectively, with $\phi _{\infty }=0.5$ and five different values of the stoichiometric parameter $\alpha$.

Figure 12. Effect of consumption of the input fluid with decreasing porosity and permeability, $\phi _{\infty }=0.5$. (a) Shape of the free surface at $t=1$ for various values of the stoichiometric parameter $\alpha =0, 0.5, 1, 2, 4$ (calculated numerically; see Appendix A.1). (b) Corresponding shape of the free surface at $t=100$. The dotted lines show the similarity solution (5.6a,b).

At early times, $t \ll 1$, there has been little microbial activity and so the results of § 3.1 for a uniform porous medium apply. At late times, $\phi \to \phi _{\infty }$, $k \to k_{\infty }$ and the governing equation (5.3) and global mass conservation (5.2) become

(5.4)$$\begin{gather} {\left[ \phi_{\infty} + \alpha (1-\phi_{\infty}) \right]\frac{\partial h}{\partial t} = k_{\infty}\frac{1}{r} \frac{\partial}{\partial r} \left(r \frac{\partial h}{\partial r} h \right),} \end{gather}$$

(5.5)$$\begin{gather}{\left[ \phi_{\infty} + \alpha (1-\phi_{\infty}) \right] \int_0^{r_f} h r \, \mathrm{d} r = t.} \end{gather}$$

The evolution is accurately described by a rescaling of the similarity solution for flow in a uniform medium, $\tilde {h}(\tilde {r},t)$ (see § 3.2),

(5.6a,b)\begin{equation} {h=k_{\infty}^{{-}1/2} \tilde{h}, \quad r= k_{\infty}^{1/4} \left[ \phi_{\infty} + \alpha (1-\phi_{\infty}) \right]^{{-}1/2} \tilde{r},} \end{equation}

where $k_{\infty }= \phi _{\infty }^3$. Figure 12(b) shows a comparison between the numerical results for the free surface (continuous lines) and the rescaled similarity solution (dotted lines) at $t=100$ for five different values of $\alpha$. In general, consumption of the input fluid reduces the radial extent of the current but not the vertical extent. The rescaling (5.6a,b) is similar to that in § 3.2 for the case of no consumption but with $\phi _{\infty }$ replaced by the effective late-time porosity $\phi _{{eff},\infty } = [ \phi _{\infty } + \alpha (1-\phi _{\infty }) ]>\phi _{\infty }$ owing to adjusted volume conservation associated with the consumption.

The sensitivity of the the effective late-time porosity, $\phi _{{eff},\infty }$, to the late-time porosity, $\phi _{\infty }$, depends on the value of $\alpha$. For $0 \leqslant \alpha < 1$, the effect of the porosity reduction associated with the growth in biomass dominates the loss of input fluid and so lower values of $\phi _{\infty }$ are associated with lower values of $\phi _{{eff},\infty }$. For $\alpha >1$, the reduction in the volume of input fluid owing to consumption is dominant and so $\phi _{{eff},\infty }$ is larger at lower values of $\phi _{\infty }$ (the volume of mobile fluid reduces faster than the pore throats constrict).

5.2. Gradual clogging of the porous medium

We revisit the analysis of § 4 but account for loss of volume of the input fluid by consumption. The porosity is given by $\phi = (1+t_R/\gamma )^{-\gamma }$ and the permeability by $k=\phi ^3$. For $\alpha >0$, at late times, the governing equation (5.3) and global mass conservation (5.2) reduce to

(5.7)$$\begin{gather} {\alpha \frac{\partial h}{\partial t} = \frac{1}{r} \frac{\partial}{\partial r} \left(r \frac{\partial h}{\partial r} \int_0^h k \, \mathrm{d} z \right),} \end{gather}$$

(5.8)$$\begin{gather}{\int_0^{r_f} \alpha h r \, \mathrm{d} r = t.} \end{gather}$$

At long residence times, the dominant contribution to the volume conservation term on the left-hand side of (5.3) is the biomass rather than the mobile fluid because the latter has mostly been consumed and converted into biomass.

The evolution is self-similar (cf. § 4) and for relatively slow clogging, $0< \gamma < 1/3$, the solution is given by

(5.9a,b)\begin{equation} {h(r,t)=t^{{3\gamma}/{2}} \alpha \gamma^{{-}3\gamma} C_0^2 \psi(\xi), \quad \xi = \frac{r}{C_0 t^{({2-3\gamma})/{4}}}.} \end{equation}

It should be noted that the time exponent of the radial extent is smaller than in the case of $\alpha =0$; see (4.5a,b).

The governing equation is recast as an integro-differential equation for $\psi (\xi )$,

(5.10)\begin{equation} {\frac{3 \gamma}{2} \psi - \left(\frac{1}{2} -\frac{3\gamma}{4} \right) \xi \frac{\mathrm{d} \psi}{\mathrm{d} \xi} = \frac{1}{\xi} \frac{\mathrm{d}}{\mathrm{d}\xi} \left( \xi J_3 \frac{\mathrm{d} \psi}{\mathrm{d} \xi} \right),} \end{equation}

which is integrated numerically; for details, see Appendix A.2. Global mass conservation (5.8) becomes

(5.11)\begin{equation} {C_0 = \left(\gamma^{{-}3\gamma} \alpha^2 \int_0^1 \psi(\xi) \xi \, \mathrm{d} \xi \right)^{{-}1/4},} \end{equation}

and so the flow thickness, $h$, is independent of $\alpha$; see (5.9a,b). The radial extent is singular as $\alpha \to 0$ and in this limit, the results of § 4 apply instead. As in § 5.1, the radial extent of the gravity current is generally reduced by consumption but the thickness is somewhat insensitive to consumption.

The aspect ratio of the gravity current becomes small at late times provided that $\gamma < 2/9$. For $\gamma >2/9$ and $\alpha >0$, the assumption of hydrostatic pressure does not apply.

6. Application

In this section, the modelling is applied to an example subsurface porous layer that could be used for underground hydrogen storage. In general, a reduction in permeability (with porosity fixed) is associated with a thicker gravity current with shorter radial extent. In contrast, a reduction in porosity (with permeability fixed) has negligible effect on the thickness but leads to a greater radial extent. Microbial growth causes a reduction in both permeability and porosity, and the total effect consists of a competition between the influence of each. For the simplified Kozeny–Carman relation used in the present work, the permeability reduction is proportional to the porosity reduction cubed ($K \sim \varPhi ^3$) and so the two effects become coupled. Consumption of the input fluid by the microbes reduces the volume of mobile fluid, but this arises primarily through a lesser radial extent rather than a change in the thickness. These results are demonstrated through the example discussed below; see also figure 13.

Figure 13. (a) Radial extent and (b) mean thickness of stored hydrogen after one month of injection. Results shown as a function of the porosity reduction factor $\phi _{\infty }$ associated with microbial growth, and the stoichiometric parameter $\alpha$ quantifying the consumption of the input fluid. The parameter values for the subsurface layer are given in table 1.

To explore the influence of hydrogen-induced microbial activity, we consider a subsurface layer with porosity that reduces with hydrogen residence time according to (2.22). Results for the radial extent and the mean thickness of the hydrogen flow after one month of injection are shown in figure 13. The properties of the flow are shown as a function of the long-time porosity reduction factor, $\phi _{\infty }$, and the stoichiometric consumption parameter $\alpha$; recall that $k_{\infty } = \phi _{\infty }^3$. The other subsurface parameters are given in table 1.

These results have some similarities to the numerical simulations of Eddaoui et al. (Reference Eddaoui, Panfilov, Ganzer and Hagemann2021) who studied the buoyant rise of hydrogen in the presence of microbial growth. They considered an unconfined porous medium without an overlying impermeable layer so the hydrogen rises vertically rather spreading horizontally as a gravity current. Nonetheless, they found that the permeability decreased most near the injection well and that this leads to the hydrogen plume having a greater cross-flow extent, which is analogous to the present result where the gravity current becomes thicker owing to the permeability reduction (see also the experiments of Ham et al. Reference Ham, Kim and Park2007).