2.1. Buoyant plume in a semi-infinite porous aquifer

Here, we derive a solution for the steady flow of a two-dimensional buoyant plume in a semi-infinite porous medium with unit thickness in the third dimension. A constant input flux $q_0$ of fluid with density $\rho _0$ is injected into a porous medium with permeability $k$ and porosity $\phi$ saturated with a denser ambient fluid of density $\rho _a$. Both fluids are assumed to be miscible, with equal viscosities such that there are no capillary effects during the flow. The implications of this assumption for application to CO$_2$–water systems is discussed in § 7. The injected fluid mixes with the ambient fluid and forms a plume with a volume flux $Q(z)$ and density $\rho (x,z)$ at a given height $z$ (figure 1a). We follow the analysis of Sahu & Flynn (Reference Sahu and Flynn2015) who considered an unconfined plume in a porous medium for Darcy flow with ${\textit {Pe}} \gg 1$, in contrast to Wooding (Reference Wooding1963). Our analysis differs by considering a half-space model, with an impermeable vertical boundary contacting the injection location. Assuming steady Boussinesq flow (shown in previous studies to be a useful approximation in CO$_2$–water systems (Soltanian et al. Reference Soltanian, Amooie, Dai, Cole and Moortgat2016; Amooie, Soltanian & Moortgat Reference Amooie, Soltanian and Moortgat2018)), the governing equations based on mass continuity, momentum continuity, solute transport and a linear equation of state are

(2.1)\begin{gather} \frac{\partial u}{\partial x} + \frac{\partial w}{\partial z} = 0, \end{gather}

(2.2)\begin{gather}\frac{1}{\rho_a}\frac{\partial P}{\partial x} + \frac{\nu}{k}u = 0, \end{gather}

(2.3)\begin{gather}\frac{1}{\rho_a}\frac{\partial P}{\partial z} + \frac{\nu}{k}w ={-}\frac{g\rho}{\rho_a}, \end{gather}

(2.4)\begin{gather}\frac{1}{\phi}\left( w \frac{\partial C}{\partial z} + u\frac{\partial C}{\partial x} \right) = \frac{\partial }{\partial x}\left(D_T \frac{\partial C}{\partial x}\right) + \frac{\partial }{\partial z}\left(D_L\frac{\partial C}{\partial z}\right), \end{gather}

(2.5)\begin{gather}\rho = \rho_a (1 - \beta C), \end{gather}

where $u$ and $w$ are the fluid velocities in the $x$ and $z$ directions respectively, $P$ is the fluid pressure, $\nu$ is the kinematic viscosity, $g$ is the acceleration due to gravity, $C$ is the solute concentration, $\beta$ is the solutal expansion coefficient and $D_L$ and $D_T$ are the longitudinal and transverse dispersion coefficients, respectively. Since the plume is long and thin we may neglect vertical variations in the horizontal velocity $({\partial w}/{\partial x} \gg {\partial u}/{\partial z})$ and longitudinal dispersion. Hence, the combined momentum equations and the solute transport equation become,

(2.6)\begin{gather} \frac{\nu}{k}\frac{\partial w}{\partial x} ={-}\frac{g}{\rho_a}\frac{\partial \rho}{\partial x}, \end{gather}

(2.7)\begin{gather}u\frac{\partial C}{\partial x} + w \frac{\partial C}{\partial z} = \phi \frac{\partial }{\partial x}\left(D_T\frac{\partial C}{\partial x}\right). \end{gather}

We now introduce a streamfunction $\psi$ such that $w = {\partial \psi }/{\partial x}$ and $u = -({\partial \psi }/{\partial z})$. For ${\textit {Pe}} \gg 1$, the transverse dispersion coefficient may be approximated as

(2.8)\begin{equation} D_T \simeq \alpha w, \end{equation}

where $\alpha$ is the transverse dispersivity (Delgado Reference Delgado2007). On using (2.5) and (2.8), (2.6) and (2.7) become

(2.9)\begin{gather} \frac{\partial^{2} \psi}{\partial x^{2}} = \frac{g\beta k}{\nu}\frac{\partial C}{\partial x}, \end{gather}

(2.10)\begin{gather}\frac{\partial \psi}{\partial x}\frac{\partial C}{\partial z} - \frac{\partial \psi}{\partial z}\frac{\partial C}{\partial x} = \phi\alpha\left( \frac{\partial^{2} \psi}{\partial x^{2}}\frac{\partial C}{\partial x} + \frac{\partial \psi}{\partial x}\frac{\partial^{2} C}{\partial x^{2}}\right). \end{gather}

For a steady plume, the buoyancy flux, or equivalently the solute mass flux, remains constant with height, as the horizontally entraining ambient fluid only increases the volume of the plume but does not alter the solute mass. The buoyancy flux of the plume per unit thickness is conserved with height and is

(2.11)\begin{equation} F_0 = \int_{0}^{\infty}w g'\,{\textrm{d}\kern0.06em x}, \end{equation}

where $g' = g(\rho _a -\rho )/\rho _a \equiv g\beta C$ is the reduced gravity. A scaling analysis of (2.9), (2.10) and (2.11) suggests that

(2.12a)\begin{gather} \frac{\psi}{x^{2}} \sim \frac{g\beta k}{\nu}\frac{C}{x}, \end{gather}

(2.12b)\begin{gather}\frac{\psi C}{xz} \sim \phi \alpha \frac{\psi C}{x^{3}}, \end{gather}

(2.12c)\begin{gather}F_0 \sim wg'x \sim \psi g\beta C. \end{gather}

Equation (2.12b) motivates us to define a self-similar horizontal length scale of the plume

(2.13)\begin{equation} \eta = \frac{x}{\sqrt{\phi\alpha z}} \quad (z > 0), \end{equation}

where we note that $z=0$ is the point of injection and therefore entrainment of the ambient fluid and the corresponding plume equations are applicable only for $z>0$. It is also worth noting that the plume is defined as the region of the flow field where the concentration $C>0$, or equivalently $g'>0$. Now, on considering (2.12a), (2.12c) and (2.13) together, we can define the streamfunction and concentration of the plume in the forms of similarity functions $\mathscr {F}(\eta )$ and $\mathscr {G}(\eta )$ as

(2.14)\begin{equation} \psi = \left[\left(\frac{F_0 k}{ \nu}\right)^{2} \phi\alpha z \right]^{1/4} \mathscr{F}(\eta), \end{equation}

and

(2.15)\begin{equation} C = \frac{1}{g\beta} \left[\left(\frac{F_0\nu}{ k}\right)^{2} \frac{1}{\phi\alpha z} \right]^{1/4} \mathscr{G}(\eta), \end{equation}

respectively. These expressions for $\psi$ and $C$ are related by (2.9) such that $\mathscr {F'}(\eta ) = \mathscr {G}(\eta )$. Similarly solute concentration, (2.10) implies that

(2.16)\begin{equation} \mathscr{F'''}\mathscr{F'}+ \mathscr{F''}\mathscr{F''}+ \tfrac{1}{4}\mathscr{F''}\mathscr{F}+ \tfrac{1}{4}\mathscr{F'}\mathscr{F'}= (\mathscr{F''}\mathscr{F'})' + \tfrac{1}{4}(\mathscr{F'}\mathscr{F})' = 0. \end{equation}

We solve (2.16) subject to the conditions that the solute concentration cannot be negative, and tends to the background concentration in the far field, $C(x,z) = 0$ and $\mathscr {F'}(\eta )=0$ as $x,\eta \to \infty$, whereas inside the plume $C(x,z)>0$ and $\mathscr {F'}(\eta ) > 0$. The transverse velocity against the fault is zero, such that the value of the streamfunction $\psi (0,z) = 0$ so $\mathscr {F}(0) = 0$. With these conditions, it can be shown that

(2.17a,b)\begin{equation} \mathscr{F}(\eta)=\begin{cases} c \sin{\dfrac{\eta}{2}},& \eta < {\rm \pi},\\ c, & \eta > {\rm \pi}, \end{cases}\quad \text{and}\quad \mathscr{G}(\eta)= \mathscr{F'}(\eta)= \begin{cases} \dfrac{c}{2} \cos{\dfrac{\eta}{2}},& \eta < {\rm \pi},\\ 0, & \eta > {\rm \pi}, \end{cases} \end{equation}

where $c=\sqrt {8/{\rm \pi} }$ is a constant of integration which is obtained by applying the solutions for $w$ $(= {\partial \psi }/{\partial x})$ and $g'$ ($=g\beta C$) to (2.11). Therefore, the expression for the volume flux per unit thickness of the plume is

(2.18)\begin{equation} Q = \int_{0}^{\infty} w\,{\textrm{d}\kern0.06em x} = \left[\left( \frac{8 F_0k}{\nu{\rm \pi}} \right)^{2} \phi\alpha z\right]^{1/4}. \end{equation}

This result is a factor of $\sqrt {2}$ smaller than the volume flux of the unconfined plume given by Sahu & Flynn (Reference Sahu and Flynn2015). Using (2.18), the average reduced gravity across the plume can be calculated as a function of height,

(2.19)\begin{equation} \bar{g'} = \frac{F_0}{Q} = \left[\left( \frac{{\rm \pi} F_0\nu}{8 k} \right)^{2} \frac{1}{\phi\alpha z}\right]^{1/4}. \end{equation}

The equations here are for an ideal plume formed purely by a buoyancy flux, where the volume flux $Q \to 0$ as $z \to 0$. For a non-ideal plume with a finite volume flux these assumptions do not hold. We can correct for this by extrapolating the flow to negative values of $z$ and defining a point $z = -z_0$ where the plume flux $Q(-z_0)=0$. The virtual source location $z_0$ is given by

(2.20)\begin{equation} z_0 = \frac{1}{\phi\alpha}\left(\frac{{\rm \pi}\nu}{8F_0k }\right)^{2} q_0^{4}, \end{equation}

where $q_0$ is the volume flux into the system, and the buoyancy flux $F_0 = q_0 g'_0$ where $g'_0$ is the reduced gravity of the injected fluid. Given these corrections for source conditions, expressions for the plume volume flux per unit thickness and mean reduced gravity are

(2.21)\begin{equation} Q = \left[\left(\frac{8 q_0g'_0k }{\nu{\rm \pi}} \right)^{2} \phi\alpha(z+z_0)\right]^{1/4}, \end{equation}

and

(2.22)\begin{equation} \bar{g'} = \left[\left( \frac{{\rm \pi} q_0g'_0\nu}{8 k } \right)^{2} \frac{1}{\phi\alpha (z+z_0)}\right]^{1/4} . \end{equation}

Figure 1.

(a) Fluid with density $\rho _0$ is injected with a constant input flux $q_0$ into a porous medium with permeability $k$ and porosity $\phi$, saturated with an ambient fluid of density $\rho _a$. The fluid forms a plume with a volume flux $Q(z)$ and density $\rho (x,z)$ at a given height $z$. (b) After the buoyant plume reaches the top of the aquifer, it provides a constant input flux $q$ of fluid with density $\rho _1$ into a gravity current which spreads into a porous medium of permeability $k$ and porosity $\phi$ saturated with an ambient fluid of density $\rho _a$. The thickness of the current is given by $h(x,t)$. The current spreads under an impermeable baffle containing a fault of gap width $d_f$, thickness $h_f$, permeability $k_f$ and porosity $\phi _f$ between $x=0$ and $x=d_f$. Injected fluid leaks through the fault with flux $q_F$.

2.2. Gravity current with leaking fault at the origin

Here, we consider the behaviour of a two-dimensional density-driven flow of a fluid in a porous medium of permeability $k$ and porosity $\phi$ saturated with an ambient fluid of higher density $\rho _a$ and bounded on one side by an impermeable vertical boundary (figure 1b). The injected fluid spreads below an impermeable horizontal baffle containing a fault of gap width $d_f$, thickness $h_f$, permeability $k_f$ and porosity $\phi _f$ which abuts the boundary. At the horizontal baffle, a gravity current forms due to fluid input from a rising plume. The fluid density $\rho _1$ and input rate $q$ into the gravity current are calculated using the expressions for the plume volume flux and mean reduced gravity derived in (2.21) and (2.22). The values of $q$ and $\bar {g'}$ are evaluated at a vertical distance $H$ from the initial injection point, which corresponds with the level of the base of the impermeable baffle. We assume that the depth of the ambient fluid is large compared with the thickness of the current $(H\gg h)$. This means we can neglect the effects of flow in the ambient and assume the values of $q$ and $\bar {g'}$ remain constant with time.

Using (2.21) and (2.22), we obtain expressions for the input flux into the current and reduced gravity of the current,

(2.23a,b)\begin{equation} q = \left[\left( \frac{8q_0g'_0k}{\nu{\rm \pi}} \right)^{2} \phi\alpha H(1+\theta)\right]^{1/4} \quad \mathrm{and} \quad g' = \left[\left( \frac{{\rm \pi} q_0g'_0\nu}{8k} \right)^{2} \frac{1}{\phi\alpha H (1 + \theta)}\right]^{1/4}, \end{equation}

where the dimensionless parameter

(2.24)\begin{equation} \theta = \frac{z_0}{H} = \frac{1}{\phi \alpha H}\left(\frac{{\rm \pi}\nu q_0}{8g'_0k}\right)^{2}. \end{equation}

At the impermeable baffle, the injected fluid forms a gravity current, with some fluid leaking through the fault. We assume that the current is long and thin, and the velocity is predominantly parallel to the baffle so that the pressure in the current is hydrostatic. The fault is modelled using the barrier-conduit system described in § 1, with the fault damage zone within the baffle modelled using an average permeability $k_f$ and width $d_f$. Neufeld et al. (Reference Neufeld, Vella and Huppert2009) formulated a model for drainage through a fissure of given permeability and width, which incorporates both flow driven by hydrostatic pressure as well as the buoyancy of the fluid in the fault itself,

(2.25)\begin{equation} w_f(t) = \frac{k_f}{\mu}\frac{{\rm \Delta} \rho g (h_0(t) + h_f)}{h_f} = \frac{k_f g' h_0(t)}{\nu h_f}\left[1+\frac{h_f}{h_0(t)}\right]. \end{equation}

Here, $k_f$ is the permeability of the fault, $h_f$ is the length of the fault, $\mu$ is the dynamic viscosity and $\nu$ is the kinematic viscosity of the leaking fluid, $h_0$ is the thickness of the current at $x=0$ $(h_0 = h(0,t))$ and $g'$ is the reduced gravity of the current. Note that, despite the width of the fault, $d_f>0$, we assume that its effect on the gravity current is localised to the point $x = 0$, and the leakage is driven by the thickness there. A sharp interface between the injected and ambient fluids is described by a thickness $h(x,t)$ below the baffle at $z = H$. Dispersion will occur predominantly at the edges of the gravity current as it propagates into the reservoir and is expected to alter the shape of the current, the implications of which are discussed during comparison with the experimental results in § 5.3. However, where the plume feeds the current directly below the fault, dispersion will be less pronounced. Hence, we neglect dispersion in the gravity current as it does not significantly affect leakage through the fault. The flow in the gravity current is driven by gradients in the hydrostatic pressure, and the evolution of the height is determined by the divergence of the fluid flux,

(2.26)\begin{equation} \phi \frac{\partial h}{\partial t} = \frac{kg'}{\nu}\frac{\partial }{\partial x} \left(h\frac{\partial h}{\partial x}\right). \end{equation}

The current forms when the plume impacts the impermeable baffle, and has initial condition

(2.27)\begin{equation} h(x,0) = 0. \end{equation}

Subsequently the gravity current is fed by the plume and has boundary conditions,

(2.28a–c)$$\begin{gather} \left[\frac{kg'}{\nu}h\frac{\partial h}{\partial x}\right]_{x=0} ={-}(q - q_F), \quad \left[\frac{kg'}{\nu}h\frac{\partial h}{\partial x}\right]_{x=x_N} = 0, \quad h(x_N,t) = 0, \end{gather}$$

which describe the input flux at the origin (equal to the plume input flux $q$, minus the fault leakage flux $q_F$), a no flux condition through the nose of the current and zero thickness at the nose of the current. The current satisfies global conservation of mass given by,

(2.29)\begin{equation} \phi \int_0^{x_N} h\,{\textrm{d}\kern0.06em x} = q t - \frac{d_f \phi_f k_f g'}{\nu h_f} \int_0^{t} h_0\left( 1 + \frac{h_f}{h_0} \right){\textrm{d}}t, \end{equation}

where the final term comes from (2.25) and is equal to the total volume of fluid that has leaked through the fault. Note that time $t=0$ in (2.29) is the instance the plume first reaches the impermeable baffle at $z=H$.

2.2.1. Non-dimensionalisation

Based on (2.23a,b), (2.26) and (2.29) we define the following dimensionless variables,

(2.30a)\begin{gather} \tilde{x} = \frac{g' k_f^{2}d_f^{2}\phi_f^{2}}{kq\nu h_f^{2}}x = \frac{{\rm \pi} k_f^{2} d_f^{2} \phi_f^{2}}{8 k^{2}h_f^{2}\phi^{1/2}\alpha^{1/2} H^{1/2}(1+\theta)^{1/2}}x, \end{gather}

(2.30b)\begin{gather}\tilde{h} = \frac{g' k_f d_f\phi_f}{q\nu h_f}h = \frac{{\rm \pi} k_f d_f \phi_f}{8 k h_f \phi^{1/2}\alpha^{1/2} H^{1/2}(1+\theta)^{1/2}} h, \end{gather}

(2.30c)\begin{gather}\tilde{t} = \frac{ g'^{2} k_f^{3} d_f^{3} \phi_f^{3}}{k\phi q\nu^{2} h_f^{3}} t = \frac{{\rm \pi}^{3/2}k_f^{3} d_f^{3}\phi_f^{3}q_0^{1/2}g_0^{\prime 1/2}}{8^{3/2}k^{5/2}h_f^{3}\nu^{1/2} \phi^{7/4}\alpha^{3/4} H^{3/4}(1+\theta)^{3/4}} t. \end{gather}

By substituting (2.30a–c), (2.26)–(2.29) become,

(2.31)\begin{gather} \frac{\partial \tilde{h}}{\partial \tilde{t}} = \frac{\partial }{\partial \tilde{x}} \left(\tilde{h}\frac{\partial \tilde{h}}{\partial \tilde{x}}\right), \end{gather}

(2.32a,b)\begin{gather}\tilde{h}(\tilde{x},0) = 0 \quad \mathrm{and} \quad \tilde{h}(\widetilde{x_N},\tilde{t}) = 0, \end{gather}

(2.33)\begin{gather}\int_0^{\tilde{x}_N} \tilde{h}\,{\textrm{d}\kern0.06em }\tilde{x} = \tilde{t} -\int_0^{\tilde{t}} \tilde{h}_0 \left( 1+ \frac{h_f}{\tilde{h}_0}\frac{k_f d_f \phi_f g'}{q \nu h_f}\right) {\textrm{d}}\tilde{t}. \end{gather}

We introduce the dimensionless parameter

(2.34)\begin{equation} \lambda = \frac{g' k_f d_f \phi_f}{q\nu} = \frac{{\rm \pi} k_f d_f \phi_f}{8 k \phi^{1/2} \alpha^{1/2} H^{1/2} (1+\theta)^{1/2}}, \end{equation}

that characterises the strength of leakage through the fault due to the buoyancy of fluid in the fault so that (2.33) now has the form,

(2.35)\begin{equation} \int_0^{\tilde{x}_N} \tilde{h}\,{\textrm{d}}\kern0.06em \tilde{x} = (1-\lambda)\tilde{t} -\int_0^{\tilde{t}} \tilde{h}_0 \,{\textrm{d}}\tilde{t}. \end{equation}

Here, $\lambda =0$ describes the case where the leakage through the fault is only driven by the hydrostatic pressure within the underlying gravity current, and not by the buoyancy of fluid in the fault.

4.1. Thinning plume model

The scenario we consider is illustrated in figure 5(a). Note the redefined position of $z=0$. The secondary plume $z\geq 0$ is fed by a vertical inlet velocity $w_f$ which is slightly smaller than the buoyancy velocity $w_b=k{\rm \Delta} \rho g/\mu$, such that $w_f=w_b(1-\epsilon )$ for some small parameter $\epsilon \ll 1$. As a result, the width of the plume, which we denote $d(z)$, must reduce from an initial value $d_f$ to a far-field value $d_b=w_f d_f/w_b$. Note that it is possible for $\epsilon$ to be negative, in which case the width of the secondary plume would increase. We model the flow into the secondary plume using Darcy's law and mass conservation,

(4.1a,b)\begin{equation} \boldsymbol{u}={-}\frac{k}{\mu} \boldsymbol{\nabla} \left( p + \rho_1 g z \right),\quad \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u}=0, \end{equation}

along with boundary conditions corresponding to constant inflow across the fault, impermeability at the vertical left-hand wall and far-field velocity conditions

(4.2)\begin{gather} w=w_b(1-\epsilon),\quad \mathrm{at}\ z=0, \end{gather}

(4.3)\begin{gather}u=0,\quad \mathrm{at}\ x=0, \end{gather}

(4.4)\begin{gather}w\rightarrow w_b,\quad \mathrm{at}\ z\rightarrow \infty, \end{gather}

respectively, where $u$ and $w$ are the fluid velocities in the $x$ and $z$ directions. At the edge of the steady plume, we impose a kinematic boundary condition and we enforce that the pressure equals the ambient hydrostatic

(4.5a,b)\begin{equation} u=w\frac{\mathrm{d} d}{\mathrm{d}z},\quad \mathrm{and} \quad p=p_a-\rho_a g z, \quad \mathrm{at}\ x=d(z). \end{equation}

The pressure at the edge of the plume is hydrostatic (4.5), since it must match with the pressure in the adjacent static ambient fluid. We therefore consider an asymptotic solution of the form

(4.6)\begin{gather} u=\epsilon \hat{u}(x,z)+\cdots, \end{gather}

(4.7)\begin{gather}w=w_b+\epsilon \hat{w}(x,z)+\cdots, \end{gather}

(4.8)\begin{gather}p=p_a-\rho_a g z+\epsilon \hat{p}(x,z)+\cdots, \end{gather}

(4.9)\begin{gather}d=d_f+\epsilon \hat{d}(z)+\cdots, \end{gather}

for $\epsilon \ll 1$. Clearly in the limit $\epsilon \rightarrow 0$ the leading-order terms in (4.6)–(4.9) satisfy (4.1)–(4.5), as expected. At first order in $\epsilon$, the pressure must satisfy the linear system

(4.10)\begin{gather} \nabla^{2} \hat{p}=0, \end{gather}

(4.11)\begin{gather}\frac{\partial \hat{p}}{\partial z}={\rm \Delta} \rho g,\quad \mathrm{at}\ z=0, \end{gather}

(4.12)\begin{gather}\frac{\partial \hat{p}}{\partial x}=0,\quad \mathrm{at}\ x=0, \end{gather}

(4.13)\begin{gather}\frac{\partial \hat{p}}{\partial z}\rightarrow 0,\quad \mathrm{at}\ z\rightarrow \infty, \end{gather}

(4.14)\begin{gather}\hat{p}=0,\quad \mathrm{at}\ x=d_f, \end{gather}

(4.15)\begin{gather}\frac{\partial \hat{p}}{\partial x}={-}{\rm \Delta} \rho g \frac{\mathrm{d}\hat{d}}{\mathrm{d}z},\quad \mathrm{at}\ x=d_f. \end{gather}

Conservation of mass dictates that the pressure must satisfy Laplace's equation, which since it is linear results in (4.10). Likewise, the boundary conditions (4.2)–(4.4) and the hydrostatic condition in (4.5b) are all linear equations. Hence, they must apply to leading- and first-order pressure terms alike. The leading-order pressure solution (first part of (4.8)) satisfies all of these trivially. The only nonlinear equation is (4.5a), also known as the kinematic condition. In this case, one must expand out the variables, keeping only first-order terms, which gives (4.15). Conveniently, (4.10)–(4.14) can be solved independently of (4.15), indicating that the plume shape and the pressure are decoupled. The pressure solution is calculated by separation of variables, giving

(4.16)\begin{equation} \hat{p}={-}\frac{8{\rm \Delta} \rho g d_f}{{\rm \pi}^{2}}\sum_{n=0}^{\infty} \frac{({-}1)^{n}}{\left( 2n+1\right)^{2}} \cos\left[\left( 2n+1\right) \frac{{\rm \pi} x}{2d_f}\right]\exp({-\left( 2n+1\right){{\rm \pi} z}/{2d_f}}). \end{equation}

Hence, the pressure (including both leading- and first-order terms) at $x=z=0$, which we denote $p^{-}$, is given by

(4.17)\begin{equation} p^{-}= p_a- {\rm \Delta} \rho g d_f \left[ \frac{8G}{{\rm \pi}^{2}} \left( 1-\frac{w_f}{w_b} \right) \right], \end{equation}

where $G=\sum _{n=0}^{\infty } (-1)^{n}/(2n+1)^{2}\approx 0.9160$ is Catalan's constant. The plume shape is calculated by evaluating (4.15) and (4.16), which converges to the differential equation

(4.18)\begin{equation} \frac{\mathrm{d}\hat{d}}{\mathrm{d}z}={-}\frac{4}{\rm \pi}\tanh^{{-}1}\left[ \exp({-{\rm \pi} z/2 d_f}) \right].\end{equation}

We solve (4.18) with boundary condition $\hat {d}(0)=0$ and plot the solution in figure 6(a). Clearly, $\hat {d}\rightarrow -d_f$ as $z\rightarrow \infty$, such that the the total plume shape $d_f+\epsilon \hat {d}$ asymptotes to the far-field plume width $d_b=d_f w_f/w_b$, as required.

Figure 6.

(a) Perturbation of the width of the secondary plume (4.18). (b,c) Contour plots of the pressure above the baffle $p^{-}$ (4.17), normalised to give a dimensionless value $(p^{-}-p_a)/{\rm \Delta} \rho g d_f$, for the different parameters of the problem $\tilde {h}_0=h_0/d_f$, $\tilde {h}_f=h_f/d_f$ and $\tilde {k}=k/k_f$. Regions of the contour plots resulting in a widening plume $w_f/w_b>1$ and positive values of $p^{-}$ are omitted.

We note that both the pressure $p^{-}$ and the inlet velocity $w_f$ are unknown in (4.17). Hence, to close the system we require a second equation relating these two quantities. This is given by considering the flow through the fault $-h_f< z\leq 0$. In particular, as illustrated in figure 5, the flow through the fault is driven by a difference in pressure $p^{+}-p^{-}$, where $p^{+}$ is set by the thickness of the gravity current

(4.19)\begin{equation} p^{+}=p_a + {\rm \Delta} \rho g h_0 + \rho_a g h_f. \end{equation}

Hence, by approximating Darcy's law across the fault we arrive at a second relationship,

(4.20)\begin{equation} w_f={-}\frac{k_f}{\mu}\left[\frac{p^{-}-p^{+}}{h_f}-\rho_1 g\right].\end{equation}

Rearranging (4.17) and (4.20), we arrive at dimensionless equations for $w_f/w_b$ and $p^{-}$,

(4.21)\begin{gather} \frac{w_f}{w_b}=\frac{1+( \tilde{h}_f+\tilde{h}_0 ){\rm \pi}^{2}/8G}{1+\tilde{k}\tilde{h}_f{\rm \pi}^{2}/8G}, \end{gather}

(4.22)\begin{gather}p^{-}= p_a- {\rm \Delta} \rho g d_f \left[\frac{\tilde{h}_f(\tilde{k}-1)-\tilde{h}_0 }{1+\tilde{k}\tilde{h}_f{\rm \pi}^{2}/8G}\right], \end{gather}

where $\tilde {h}_0=h_0/d_f$, $\tilde {h}_f=h_f/d_f$ and $\tilde {k}=k/k_f$ are known parameters. From (4.21) and (4.22) it is clear that negative values of $p^{-} - p_a$ correspond to $w_f/w_b<1$, and positive values of $p^{-} - p_a$ correspond to $w_f/w_b>1$.

The pressure $p^{-}$ (written in dimensionless form $(p^{-}-p_a)/{\rm \Delta} \rho g d_f$) is illustrated in figures 6(b) and 6(c) using contour plots. Largest negative pressures are observed for small values of the gravity current thickness $\tilde {h}_0$, or large values of $\tilde {h}_f$ and $\tilde {k}$ (corresponding to small velocity ratios $w_f/w_b$). We note, however, that we do not expect our asymptotic solution to be valid for $w_f/w_b\ll 1$. In such situations a full numerical simulation is required to determine $p^{-}$. Nevertheless, for the parameter range of the current study ($w_f/w_b\in [0.5,1.2]$) we expect (4.22) to remain a good approximation.

5.1. Experimental set-up

The experiments were performed within a Perspex cell of length 40 cm, height 70 cm and internal thickness 1 cm (figure 7). The cell was filled with glass ballotini which formed a porous layer. The glass ballotini filling the majority of the cell had diameter $b_0 = 3.1 \pm 0.2$ mm except for the region adjacent to a central plastic spacer that lies across the cell where the ballotini diameter $b_f = 1.0 \pm 0.1$ mm. The porosity of 3 mm glass ballotini in a cell of width 1 cm was measured by Sahu & Neufeld (Reference Sahu and Neufeld2020) and found to be $\phi = 0.41 \pm 0.01$, which is slightly higher than the value of $\phi = 0.37$ for randomly close-packed beads as the width of the cell is comparable to the bead size. The permeability of the larger porous medium filling the majority of the cell was estimated using the Kozeny–Carman relation

(5.1)\begin{equation} k = \frac{b_0^{2}}{180}\frac{\phi^{3}}{(1-\phi)^{2}} \approx 1.06 \pm 0.15 \times 10^{{-}8}\,\text{m}^{2}. \end{equation}

The Kozeny–Carman relationship breaks down when applied to a small volume of beads, hence the permeability of the region with smaller ballotini was treated as an unknown parameter and fitted for (see § 5.3). The entire cell submerged in a large water tank with dimensions $90\,\textrm {cm} \times 80\,\textrm {cm}\times 80\,\textrm {cm}$ filled with fresh water of density $\rho _a = 0.998$ g cm$^{-3}$. The cell was closed on three sides, with one side open but lined with a permeable mesh so that this unconfined side created a hydrostatic pressure boundary condition at the right-hand side of the cell. A total of 15 experiments were conducted, using four different fault geometries (A–D), as summarised in table 1.

Figure 7.

Schematic of the experimental set-up.

Table 1.

Parameter values used in our experiments with associated value of $\lambda$ from (2.34). Typical uncertainties in these measurements are $d_f \pm 0.05$, $h_f \pm 0.05$, $H \pm 0.1$ cm, $q_0 \pm 0.001$ cm$^{3}$ s$^{-1}$ and $\rho _0 \pm 0.005$ g cm$^{-3}$.

During the experimental run, a peristaltic pump was used to inject aqueous solutions of sodium chloride (brine) dyed with red food colouring into the cell at a constant rate $q_0$ which ranged from 0.039–0.245 cm$^{3}$ s$^{-1}$, and was calculated by measuring the mass of the brine container over the course of the experiment. The injected brine solutions had an initial dye concentration of $c_0 = 3.00 \pm 0.05$ g l$^{-1}$ and an initial density $\rho _0 = 1.030$ or 1.070 g cm$^{-3}$, assuming a water temperature of 20 $^{\circ }$C (Green & Southard Reference Green and Southard2019). The kinematic viscosity of the water was 0.01 cm$^{2}$ s$^{-1}$ and the transverse dispersivity of the glass ballotini was assumed to be $\alpha =0.005$ cm (Delgado Reference Delgado2007). The experimental parameters were selected so experiments spanned a range of $\lambda$ values. We used a Nikon D5000 DSLR camera with a resolution of $4288 \times 2848$ pixels to capture images over the course of the experiment with a time gap which ranged from 4 to 10 s. The cell was backlit using a LED light panel with a Perspex diffuser to ensure uniform illumination.

5.2. Post-processing scheme

A set of photographs for experiment B2 is shown in figure 8, rotated by 180$^{\circ }$ for ready comparison between the analytical and experimental results. The theoretical predictions for the position of the gravity current interface $h(x,t)$ in the bottom portion of the cell is plotted for the thinning plume model, along with the interface for the zero leakage model (Huppert & Woods Reference Huppert and Woods1995). The time $t=0$ is defined as the point where the plume first makes contact with the bottom of the central spacer. The horizontal extent of the gravity current was measured as a function of time by picking the furthest front of the current above a threshold value. The error range of the measurements was set by calculating the extent for a range of threshold values. The results were verified by manually checking the front location from photographs and showed excellent agreement. The height of the current was more difficult to interpret due to dispersion from the plume making the top of the current uncertain.

Figure 8.

Images showing the evolution of experiment B2. The theoretical predictions for the position of the interface $h(x,t)$ are shown for the thinning plume model presented in § 4 and the zero leakage model (Huppert & Woods Reference Huppert and Woods1995). Images are rotated 180$^{\circ }$ to allow comparison with theoretical results (see § 5).

To measure the flux of fluid leaking through the fault, the concentration of the injected fluid throughout the cell must be determined, and so a series of calibration experiments were performed to determine the functional relationship between the image intensity and dye concentration. In each calibration experiment, the cell was uniformly saturated with a red dye solution with concentration $C_d$. The light intensity for the calibration image was calculated by subtracting the light intensity from a reference image containing no dye. A total of 20 different concentrations between 0 and $3.00 \pm 0.05$ g l$^{-1}$ were used to construct calibration curves for the green and blue colour channels (figure 9). The dye calibration is more sensitive to different colour channels at different dye concentrations, so a hybrid calibration curve was used which weighted contributions from the green and blue curves. The full function form of the hybrid calibration curve is given in Appendix B.

Figure 9.

Calibration curves of image intensity against dye concentration for the green and blue colour channels (black dashed lines). Also plotted is the effective concentration where both calibration curves have equal weighting (grey solid line) and the 99 % intervals above and below which the green and blue curves dominate the effective concentration (grey dashed lines). See Appendix B for functional forms of the curves.

The mass of dye across the cell is calculated by summing the average concentration within $0.5\,\text {cm} \times 0.5\,\text {cm}$ sized bins, with the leaked mass equal to the total mass below the top of the baffle. The total error in measurement is the sum of two errors. First, the leaked dye mass shows a linear behaviour at late times. A root-mean-square deviation from a regression line was calculated, with the deviations likely caused by small variations in light intensity between each photo, for example due to trapped air bubbles. Second, the total measured mass of dye across the cell was compared with the known input dye mass for each experiment. The measured mass increases linearly as a function of time but the post-processing recovers 83 %–97 % of the input dye mass across experiments. The partial measurement of the input dye is likely due to sensitivity of the calibration curve to smaller concentrations. For each experiment, the measured dye mass is scaled to the known input dye mass and the difference is defined as the error in measurement.

5.3. Experimental results and comparison with theory

Figure 10(a) shows the total leaked dye mass (black circles) and leakage flux (pink circles) as a function of time for experiment A3. The leaked flux is calculated by taking the gradient of a quadratic polynomial fitted over a seven element moving window. The errors bars represent a 95 % confidence interval for the regression coefficients.

Figure 10.

(a) Total leaked dye mass (green circles) and leakage flux (pink circles) plotted as a function of time for experiment A3 along with results for three leakage models (green and pink lines). Grey dashed lines display thinning plume model solutions for an uncertainty range for input $k_f \pm 20\,\%$. (b) Picture of the experiment at time $t=700$ s with width of fault $d_f$ and width of secondary plume $d_b$ indicated.

The experimental results are compared with theoretical results from three different models presented in §§ 2 and 4. The first model only considers the contribution from the hydrostatic pressure within the underlying current ($\lambda = 0$). The second model considers contributions from the underlying current and the buoyancy of fluid in the fault ($\lambda \neq 0$), and the third model includes both of these effects but also considers the contribution of flow-enhancing pressure deviations directly above the fault (thinning plume model). There is very little difference between the solutions for the second and third models, suggesting that the effects of pressure deviations above the baffle due to thinning of the secondary plume are minimal. This agrees with experimental observations (figure 10b), which show that the width of the secondary plume is the same as the width of the fault ($d_b/d_f \approx 1$).

Due to the difficulty in estimating the fault permeability $k_f$, this was calculated by minimising the misfit between the experimental data and the thinning plume model using a least-squares regression for one of the experiments. Experiment A3 (figure 10) was selected to fit the fault permeability due to the agreement between the thinning plume model and the $\lambda \neq 0$ model. A value of $k_f = k_0/3.9 \approx 2.7\pm 0.3\times 10^{-9}\,\text {m}^{2}$ was obtained and this value was used to calculate the theoretical results across the other experiments. It is possible that there is some variation in the fault permeability across different experiments, but this should be small as the faults are a similar size and were all packed using the same method. The sensitivity of the thinning plume model to changes in fault permeability is shown by the grey dashed lines in figure 10(a), in which the model solutions are displayed for $k_f \pm 20\,\%$.

Figure 11(a) shows the total leaked dye mass (black circles) and leakage flux (pink circles) as a function of time for experiment D4, along with theoretical results for the three models. There is a clear difference between the model which allows for a thinning plume above the baffle and the model only considering contributions from the underlying current and buoyancy in the fault, suggesting that acceleration and thinning of the secondary plume is leading to enhanced pressure gradients within the fault. Experimental observations are in agreement (figure 11b), where significant thinning of the secondary plume can be seen above the fault ($d_b/d_f < 1$).

Figure 11.

(a) Total leaked dye mass (green circles) and leakage flux (pink circles) plotted as a function of time for experiment D4 along with results for three leakage models (green and pink lines). (b) Picture of the experiment at time $t=660$ s with width of fault $d_f$ and width of secondary plume $d_b$ indicated.

Note that the model considering buoyancy in the fault assumes that the fault is initially full of the injected fluid. Furthermore, the thinning plume model assumes that a secondary plume has reached a quasi-steady state profile above the fault. In reality, at early times the fault remains filled with the ambient fluid so the only driving force is the hydrostatic pressure in the underlying current. To account for this, a breakthrough time is introduced, defined as the time it takes for the injected fluid to fill the fault and breakthrough into the upper reservoir, which is obtained from experimental observations. Prior to the breakthrough time, all three models only consider the contribution from the underlying current. After the breakthrough time, the other driving forces are taken into account. Further details on how the breakthrough time is calculated and the sensitivity of the model to the breakthrough time is given in the supplementary material available at https://doi.org/10.1017/jfm.2021.970. This modification to the theory results in a discontinuous leakage flux profile (pink dashes, figures 10(a), 11(a)), where the leakage flux increases rapidly at the breakthrough time. These predictions broadly agree with the behaviour seen in the experimental data, where initial leakage rates are slow as the fault fills with the injected fluid before increasing rapidly to a higher and constant rate. The thinning plume model shows good agreement with the experimental data in both total leaked dye mass and leakage flux.

Figure 12 shows the total leaked dye mass as a function of time across all experiments, with symbols listed in table 1. The solutions for the thinning plume model are plotted for each experiment and show good agreement across all experimental set-ups (a–d).

Figure 12.

Total leaked volume through the fault plotted as a function of time for the experimental set-ups (a–d), with symbols listed in table 1. Theoretical predictions of the leakage volume for each experiment calculated from the thinning plume model plotted as dashed lines.

The horizontal extent of the underlying gravity current along with solutions from the thinning plume model is plotted in figure 13 for experimental set-ups A–D. In general the theoretical model shows good agreement with the experimental data, although it tends to overpredict the horizontal extent of the current. A possible explanation is that the gravity current entrains ambient fluid as it moves into the reservoir which is not captured in the model. This decreases the reduced gravity of the current, and hence reduces the distance it travels. This effect would also cause the thickness of the current to increase, and this can be seen in figure 8, where the experimental current has a greater thickness than the thinning plume model predicts, but appears more dispersed towards its outer edge.

Figure 13.

Evolution of horizontal extent of the gravity current $x_N$ plotted for the experimental set-ups (a–d), with symbols listed in table 1. Theoretical predictions for the extent growth plotted as black dashed lines.