2.1. Steady plumes

To address the steady-state regime, we describe the shape of the plume $A(z)$ above the injection height ($z\geq 0$) in general terms that apply to both linear and circular sources. At the injection point, the vertical inflow velocity is

(2.1)\begin{equation} w_0=\frac{Q}{A_0}.\end{equation}

Likewise, the buoyancy velocity of the injected fluid, which is the equilibrium rise speed in the porous medium (i.e. buoyancy balancing viscous resistance), is given by

(2.2)\begin{equation} w_b=\frac{k\,{\rm \Delta}\rho\,g}{\mu},\end{equation}

where $k$ is the permeability of the medium, and ${\rm \Delta} \rho =\rho _2-\rho _1$. Due to mass conservation, the far field cross-section is given by $A_\infty =Q/w_b$. (Note that we use the term ‘far field’ here and throughout the paper to refer to the length scale over which the plume velocity approximately matches the buoyancy velocity. This is not to be confused with the even greater length scales over which the effects of mixing are important, since these are not studied here.) Hence the key dimensionless parameter in this study is the ratio between the inlet and buoyancy velocities,

(2.3)\begin{equation} W=\frac{w_0}{w_b}. \end{equation}

As the flow moves downstream (i.e. upwards), the plume must become thinner ($A_\infty < A_0$) for sub-buoyancy velocities $W<1$ and thicker ($A_\infty >A_0$) for super-buoyancy velocities $W>1$.

The flow within the injected fluid is governed by the Darcy equations

(2.4)\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}=0, \end{gather}

(2.5)\begin{gather}\boldsymbol{u}={-}\frac{k}{\mu}\left( \boldsymbol{\nabla} p + \rho_1 g \hat{\boldsymbol{k}}\right), \end{gather}

where $\boldsymbol {u}$ is the Darcy velocity vector, and $p$ is the pressure. The flow in the ambient fluid is coupled to the injected fluid only via the boundary conditions, which we discuss in the next subsection. Hence, for the purposes of this study, we omit further details of how to model the ambient flow outside the injected region, since the behaviour of the injected flow is of primary interest.

The Darcy equations (2.4)–(2.5) are accompanied by boundary conditions that take a different form depending on whether the flow is injected from a line source or a circular source. Hence we address the former and latter cases separately in §§ 2.2 and 2.4.

2.2. Plume shape: the case of a line source

In the case of the line source (with $A(z)=a(z)\,d$, where $d$ is the depth in the third dimension), we impose boundary conditions within the injected fluid region of the form

(2.6)\begin{gather} u=0,\quad x=0, \end{gather}

(2.7)\begin{gather}w=w_0,\quad z=0, \end{gather}

(2.8)\begin{gather}w\rightarrow w_b,\quad z\rightarrow \infty, \end{gather}

(2.9)\begin{gather}u=w\,a'(z),\quad x=a(z), \end{gather}

(2.10)\begin{gather}p=p_a-\rho_2 g z,\quad x=a(z), \end{gather}

where $p_a$ is the ambient hydrostatic pressure at $z=0$, and $a'(z)=\mathrm {d}a/\mathrm {d}z$. The above boundary conditions correspond with imposing symmetry on the $z$ axis (2.6), constant inflow at the source (2.7), matching with the far-field buoyancy velocity (2.8), and applying the kinematic (2.9) and dynamic conditions (2.10) at the sharp interface. For further details on the governing equations and boundary conditions for flow in porous media, see Bear (Reference Bear2013).

The above system (2.4)–(2.10) is a free boundary problem for both the flow and the shape of the interface $a(z)$. To proceed, we seek a solution of the form

(2.11)\begin{gather} p=p_a-\rho_2 g z +\hat{p}, \end{gather}

(2.12)\begin{gather}u=\hat{u}, \end{gather}

(2.13)\begin{gather}w=w_b+\hat{w}, \end{gather}

where the hatted Darcy velocities satisfy

(2.14)\begin{gather} \hat{u}={-}(k/\mu) \hat{p}_x, \end{gather}

(2.15)\begin{gather}\hat{w}={-}(k/\mu) \hat{p}_z, \end{gather}

in which subscripts denote partial derivatives. Hence the hatted pressure satisfies the new system of equations

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

(2.17)\begin{gather}\hat{p}_x=0,\quad x=0 , \end{gather}

(2.18)\begin{gather}\hat{p}_z={\rm \Delta}\rho\,g(1-W),\quad z=0, \end{gather}

(2.19)\begin{gather}\hat{p}_z\rightarrow 0,\quad z\rightarrow \infty, \end{gather}

(2.20)\begin{gather}\hat{p}_x=(\hat{p}_z-{\rm \Delta} \rho\,g)\,a'(z),\quad x=a(z), \end{gather}

(2.21)\begin{gather}\hat{p}=0,\quad x=a(z). \end{gather}

The governing equation (2.16) is obtained by inserting (2.12)–(2.15) into the continuity equation (2.4), whereas (2.17)–(2.21) are obtained from each of the boundary conditions (2.6)–(2.10), respectively.

In the case where the inlet velocity is close to the buoyancy velocity ($W\approx 1$), the solution to (2.16)–(2.21) was calculated by Gilmore et al. (Reference Gilmore, Sahu, Benham, Neufeld and Bickle2022) using the method of separation of variables, in which case

(2.22) \begin{equation} \hat{p}\approx{-}(1-W){\frac{8\,{\rm \Delta} \rho\,g a_0}{{\rm \pi}^2}}\sum_{n=0}^\infty {\frac{({-}1)^n}{(2n+1)^2}} \cos\left[{(}2n+1{)}{{\rm \pi} {x}}/{2a_0}\right]\exp\bigl(-(2n+1){{\rm \pi} z}/{2a_0}\bigr),\end{equation}

where $a_0=a(0)$. Likewise, the approximate plume shape is given by the solution to

(2.23)\begin{equation} a'(z)\approx{-}(1-W)\,\frac{4}{\rm \pi}\tanh^{{-}1}\left[\exp(-{\rm \pi} z/2 a_0)\right].\end{equation}

This is obtained from the linearised version of (2.20) (i.e. $\hat {p}_x=-{\rm \Delta} \rho \,g\,a'(z)$ at $x=a_0$) after inserting the formula for the pressure (2.22). Then the infinite sum converges to the inverse hyperbolic tangent function in (2.23).

In the case where the inlet velocity and the buoyancy velocity are not similar ($W\not \approx 1$), a numerical method must be used to calculate the solution. By converting to a set of scaled dimensionless variables

(2.24a–d)\begin{equation} X=x/a(z),\quad Z=z/a_0,\quad P(X,Z)=\hat{p}(x,z)/{\rm \Delta} \rho\,g a_0,\quad \alpha(Z)= a(z)/a_0,\end{equation}

Laplace's equation (2.16) becomes

(2.25)\begin{equation} \left[\alpha^{{-}2}\,\partial_{XX}+\left(\partial_Z-X\alpha'\alpha^{{-}1}\,\partial_X\right)^2\right]P=0,\end{equation}

whilst the remaining boundary conditions become

(2.26)\begin{gather} P_X=0,\quad X=0, \end{gather}

(2.27)\begin{gather}P_Z-X\alpha'\alpha^{{-}1} P_X=1-W,\quad Z=0, \end{gather}

(2.28)\begin{gather}P_Z-X\alpha'\alpha^{{-}1} P_X\rightarrow 0,\quad Z\rightarrow \infty, \end{gather}

(2.29)\begin{gather}P_X={-}\alpha\alpha'- {\alpha'}^2P_X,\quad X=1, \end{gather}

(2.30)\begin{gather}P=0:\quad X=1. \end{gather}

The nonlinear system (2.25)–(2.30) is solved using Newton's method in combination with a finite difference scheme. The domain is discretised using a rectangular grid, with $X\in [0,1]$ and $Z\in [0,H]$, where the boundary condition (2.28) is approximated at a large but finite value $H=10$. We calculate the solution for a variety of values of $W$ (the only dimensionless parameter of the problem), using an 8th-order finite difference scheme with a grid of $20\times 200$ points in the $X,Z$ directions (Note that in some cases – such as for $W>1$ – the discretisation was increased to $40\times 400$ points to resolve accurately sharp gradients in the plume shape near the origin.) By employing the method of continuation using incremental changes in $W$, Newton's method converges in approximately three steps (for each increment).

We plot examples of the shape $\alpha (Z)$ in figure 2(a) for $W=0.3,0.6,0.9,1.1$. Thinning plumes are observed for $W<1$, whereas thickening plumes are observed for $W>1$, as expected. One salient feature of the analysis is the distance over which the plume approaches its far-field width $a_\infty$ (or $\alpha \rightarrow W$ in dimensionless terms). We define the $99\,\%$ boundary layer distance $\delta$ (dimensionless) as

(2.31)\begin{equation} \left|\frac{\alpha(\delta)-W}{1-W}\right|=0.01,\end{equation}

and this is plotted in figure 2(b) for different values of $W$. We find that $\delta (W)$ is monotone increasing for $W\in [0.12,1.14]$. Outside this range, our numerical method fails to converge to a real-valued solution (which we discuss shortly), so no data are plotted. We also compare these values of $\delta$ with the analytical value in the case where $W$ is close to $1$ (e.g. via the solution to (2.23)). In this case, the boundary layer distance is independent of $W$ at leading order, and is given by the approximate value $\delta \approx 2.83$ (see dotted lines in figure 2b). We note that despite the fact that the boundary condition (2.28) is imposed at $Z\rightarrow \infty$, no more than a few plume-width scales are required for the plume to adjust to its far-field shape ($99\,\%$ of the way, more specifically), as seen by the $\delta$ values in figure 2(b).

Figure 2. Numerical and analytical results for a thinning/thickening plume resulting from a line source. (a) Plume shape for different velocity ratios $W=w_0/w_b$, and (b) $99\,\%$ boundary layer distance $\delta$ as defined in (2.31). Critical values $W_0$, $W^*$ and $W_n$ are related to the sign of the discriminant $\varDelta$ (see (2.32)) and the stability/existence of a steady solution. Dotted lines indicate the approximate solution when $W\approx 1$, for which the plume shape is given by the solution to (2.23).

Due to the kinematic boundary condition (2.29), for the solution to remain real-valued we require a non-negative discriminant

(2.32)\begin{equation} \Delta(Z):=\alpha^2-4\left.P_X\right|_{X=1}^2\geq 0,\end{equation}

for all values of $Z$. By writing the pressure gradient as a dimensionless velocity $P_X=-U$, we see that (2.32) can be interpreted as a balance criterion between the width of the plume $\alpha$ and the horizontal velocity $U$ required to sustain that width. For example, a thinning plume ($U<0$) with a shape that tapers to smaller than thickness $\alpha <-2U$ is not permitted by (2.32). In general, whenever (2.32) cannot be satisfied, this indicates that a smooth and continuous steady plume shape is not possible.

The behaviour of the discriminant function (2.32) depends on the value of the velocity ratio $W$. There are four solution regimes defined by three values of the velocity ratio, $W^*\approx 0.5$, $W_n\approx 1.14$ and $W_0\approx 0.12$. For velocity ratios $W^*< W< W_n$, the discriminant is strictly positive, $\varDelta >0$, for all values of $Z$. For velocity ratios in the range $W_0< W< W^*$, the discriminant is non-negative, $\varDelta \geq 0$, but equals zero at some critical distance $Z=Z^*$ downstream of the inflow. For $0< W< W_0$ or $W>W_n$, Newton's method fails to converge to a real-valued solution, indicating that a steady solution in the form of a smooth continuous shape $\alpha (Z)$ may not exist. The four different solution regimes are illustrated with shading in figure 2(b). It should be noted that the sign of the discriminant is closely linked with the stability criteria for the plume, and we discuss this later, in § 3.

Before discussing the comparison with experiments, it is interesting to note the possible effects of the ambient fluid that we have so far neglected in this study. Within the ambient region, the flow satisfies a set of Darcy equations similar to (2.4)–(2.5) except with $\rho =\rho _2$. Appropriate boundary conditions consist of the impermeability condition at the bottom boundary ($w=0$ on $z=0$), far-field hydrostatic conditions ($p \rightarrow p_a-\rho _2 g z$, as $x,z,\rightarrow \infty$), and the kinematic and dynamic boundary conditions, (2.9) and (2.10), at the interface $x=a(z)$. One can see immediately that a hydrostatic pressure profile $p=p_a-\rho _2 g z$ satisfies all of these conditions. Hence the flow in the ambient region is simply zero everywhere in the steady state, and therefore affects the injected fluid region only via the density difference.

It is also interesting to note that if we change the viscosity of the ambient fluid (i.e. such that $\mu _2\neq \mu _1$), then this has no effect on the preceding argument. Therefore, it follows that the steady plume shape is independent of the viscosity ratio between the injected and ambient fluids, $M=\mu _1/\mu _2$. However, as we discuss later, in § 3, this is no longer true in the unsteady case.

2.3. Comparison with experiments

In this subsection, we compare our results for the steady plume shape (in the case of a line source) to the porous bead experiments of Gilmore et al. (Reference Gilmore, Sahu, Benham, Neufeld and Bickle2022). These experiments were conducted in a thin rectangular tank of dimensions $40\times 70\,\mathrm {cm}^2$ in the $x,z$ directions and $1\,\mathrm {cm}$ thick in the transverse ($y$) direction. The tank was filled with $3\,\mathrm {mm}$ Ballotini beads and initially saturated with fresh water. Salty water dyed with red food colouring was injected into the top of the tank, using different salt concentrations to modify the density contrast. Since salty water is heavier than fresh water, the experiments resulted in a falling plume rather than the rising plume studied at present. Therefore, we have inverted their experimental photos for comparison with our model. The inverted system behaves in approximately the same way as the current system due to the Boussinesq approximation (Soltanian et al. Reference Soltanian, Amooie, Dai, Cole and Moortgat2016; Amooie et al. Reference Amooie, Soltanian and Moortgat2018).

The focus of the Gilmore et al. (Reference Gilmore, Sahu, Benham, Neufeld and Bickle2022) study was on the leakage of salty water through a gap in an impermeable division midway down the tank. However, for comparison with the present study, we focus on the flow below this division only, and we ignore all the flow details above this. Therefore, we restrict our attention to the lower $40\times 40\times 1\,\mathrm {cm}^3$ of their tank. In this way, the leakage rates of salty water into this lower section of the tank (which were calculated in their study) correspond to the injection flow rate $Q$ in our theoretical model. As described by Gilmore et al. (Reference Gilmore, Sahu, Benham, Neufeld and Bickle2022), both the leakage flux and the near-field plume shape within this lower section of the tank were approximately steady (after an initial transient). Hence, for comparison with our model, a constant inflow $Q$ and a steady shape $\alpha (Z)$ can be assumed to good approximation.

Examples of steady plumes from the study of Gilmore et al. (Reference Gilmore, Sahu, Benham, Neufeld and Bickle2022) are shown in experimental photographs in figures 3(a,b). The plume width at the inlet for each case is $a_0=3,2.5$ cm, respectively. To calculate the velocity ratios $W$ for each case, we estimate the far-field plume width $a_\infty$ downstream, noting that $W=a_\infty /a_0$. This results in $W=0.58$ and $W=0.74$ for figures 3(a) and 3(b), respectively. We have also tried estimating $W$ using values of $w_0$ and $w_b=k\,{\rm \Delta} \rho \,g/\mu$ given by Gilmore et al. (i.e. using quoted values for the permeability, density, viscosity and leakage flux), which produces similar calculations.

Figure 3. Experimental photos (taken from the study of Gilmore et al. Reference Gilmore, Sahu, Benham, Neufeld and Bickle2022) of thinning plumes with (a) $W=0.58$ and (b) $W=0.74$, compared with numerical (solid lines) and analytical (dotted lines) solutions for the steady plume shape in the case of a line source. The photos are partially obscured by a clamp (part of the apparatus), which is labelled for clarity.

The steady plume shape predicted by our numerical model is compared to each of these photos with solid lines. The analytical approximation, given by integrating (2.23), is also shown with dotted lines. Overall, good agreement is observed between the numerical model, the analytical approximation and the experiments. Dispersion causes the plume shape to slightly diffuse downstream of the inlet, which is not captured by the sharp interface in our model.

One advantage of our simple model is that it depends on only a single dimensionless parameter $W$, which is easily calculated by estimating the plume width at two locations (i.e. $a_0$ and $a_\infty$). A more complicated model that accounts for dispersion, for example, would require further parameter values of the fluid-medium properties, such as the diffusion and dispersion coefficients of the salt/dye.

2.4. Plume shape: the case of a circular source

In the case of a circular source, the boundary conditions (2.6)–(2.10) are replaced by corresponding conditions in cylindrical radial coordinates (e.g. with $x$ replaced by $r$, and $u$ replaced by $u_r$, the radial velocity). In this case, the radius of the plume (measured from the $z$ axis) is given by $a(z)=(A(z)/{\rm \pi} )^{1/2}$. As before, we seek a solution of the form (2.11)–(2.13) (with $u$ replaced by $u_r$). The pressure $\hat {p}$ satisfies a system of equations similar to (2.16)–(2.21), with $x$ replaced by $r$. In the case where $W$ is close to unity, the solution is calculated by separation of variables (see Appendix A), giving

(2.33)\begin{equation} \hat{p}\approx{-}2(1-W)\,{\rm \Delta} \rho\,g a_0\sum_{n=1}^\infty \frac{J_0(j_{0,n}r/a_0)}{j_{0,n}^2\,J_1(j_{0,n})} \exp({-j_{0,n}z/a_0}),\end{equation}

where $J_0$ and $J_1$ are the 0th- and 1st-order Bessel functions of the first kind, and $j_{0,n}$ is the $n$th zero of $J_0$. Likewise, the plume shape is given by the solution to

(2.34)\begin{equation} a'(z)\approx{-}2(1-W)\sum_{n=1}^\infty \frac{\exp({-j_{0,n}z/a_0})}{j_{0,n}}.\end{equation}

This is obtained from the linearised version of (2.20) in radial coordinates (i.e. $\hat {p}_r=-{\rm \Delta} \rho \,g\,a'(z)$ at $r=a_0$) after inserting the formula for the pressure (2.33). By using the identity $J_0'(r)=-J_1(r)$, this simplifies to (2.34).

In the case where $W$ is not close to unity, we calculate the solution via the numerical method described earlier. After introducing a scaled dimensionless radial coordinate

(2.35)\begin{equation} R=r/a(z),\end{equation}

Laplace's equation (2.16) becomes

(2.36)\begin{equation} \left[\alpha^{{-}2}R^{{-}1}\,\partial_{R}(R\,\partial_R)+\left(\partial_Z-R\alpha'\alpha^{{-}1}\, \partial_R\right)^2\right]P=0,\end{equation}

whilst the remaining boundary conditions stay the same as (2.26)–(2.30) except with $X$ switched to $R$. The system of equations is solved using the same finite difference scheme as in § 2.2 (except with the radius truncated at a small but finite value $R=0.01$ to avoid a singular Laplacian). Likewise, a similar discriminant function is defined as

(2.37)\begin{equation} \Delta(Z):=\alpha^2-4\left.P_R\right|_{R=1}^2,\end{equation}

which indicates whether or not a real solution exists.

In figure 4(a), we display the steady plume shapes calculated for $W=0.3,0.6,0.9,1.1$. Likewise, the $99\,\%$ boundary layer distance $\delta$ (as defined in (2.31)) is plotted in figure 4(b). Overall, the behaviour is similar to the case of a line source, except that the plume adjusts over a shorter vertical length scale, resulting in smaller values of $\delta$. The approximate solution calculated in the case where $W\approx 1$ (i.e. via integration of (2.34)) results in a boundary layer distance $\delta \approx 1.76$, as shown with dotted lines in figure 4(b). Critical values of the velocity ratio (see earlier discussion in § 2.2) are $W_0=0.1$, $W^*=0.5$ and $W_n=1.17$, which are very similar to the case of a line source.

Figure 4. Numerical and analytical results for a thinning/thickening axisymmetric plume resulting from a circular source. (a) Plume shape for different velocity ratios $W$, and (b) $99\,\%$ boundary layer distance $\delta$ as defined in (2.31). Critical values $W_0$, $W^*$ and $W_n$ are related to the sign of the discriminant $\varDelta$ in (2.37) and the stability/existence of a steady solution. Dotted lines indicate the approximate solution when $W\approx 1$, for which the plume shape is given by the solution to (2.34).

3.1. Small velocity ratios $W_0< W< W^*$

As described earlier, it is expected that the flow may become unstable for small velocity ratios, such that an unsteady model is required. In the unsteady case, the only equation that requires modification is the kinematic boundary condition (2.9), which becomes

(3.1)\begin{equation} u=\phi a_t + w a_z,\quad x=a(z,t),\end{equation}

where $\phi$ is the porosity. Written in terms of the stretched dimensionless coordinates (2.24a–d), this becomes

(3.2)\begin{equation} -\alpha^{{-}1}P_X=\alpha_T + \left[ 1+ \alpha_Z \alpha^{{-}1}P_X \right]\alpha_Z,\quad X=1,\end{equation}

where $T=t w_b/a_0\phi$ is the dimensionless time. The rest of the governing equations and boundary conditions remain the same as in the steady case, i.e. (2.25)–(2.28) and (2.30).

We consider a small perturbation applied to the plume shape and pressure of the form

(3.3)\begin{gather} \alpha=\bar{\alpha}(Z)+\epsilon\,\tilde{\alpha}(Z,T), \end{gather}

(3.4)\begin{gather}P=\bar{P}(X,Z)+\epsilon\,\tilde{P}(X,Z,T), \end{gather}

where $\epsilon \ll 1$ is a small parameter, and $\bar {\alpha },\bar {P}$ solve the steady problem. Inserting (3.3) and (3.4) into (3.2), and linearising, we get

(3.5)\begin{equation} \tilde{\alpha}_T+(1+2\bar{\beta}\bar{P}_X)\tilde{\alpha}_Z -\bar{P}_X(\bar{\beta}^2+\bar{\alpha}^{{-}2}) \tilde{\alpha} +(\bar{\beta}^2\bar{\alpha}+\bar{\alpha}^{{-}1})\tilde{P}_X=0:\quad X=1,\end{equation}

where we have introduced the notation $\bar {\beta }=\bar {\alpha }'/\bar {\alpha }$. In addition to (3.5), we require a set of equations and boundary conditions for the perturbed pressure $\tilde {P}$ to complete the system. These are placed in Appendix B for convenience.

The stability of the system is elucidated by considering the discriminant function (2.32) at leading order, $\bar {\varDelta }=\bar {\alpha }^2-4\bar {P}_X^2$. As described earlier, for small values of the velocity ratio $W_0< W< W^*$, the discriminant function becomes zero at a critical point $Z=Z^*$ downstream of the inlet. To illustrate this, we have plotted the discriminant function in figure 5(a) for $W=0.2,0.3,0.4,0.5$, and the corresponding critical points $Z^*$ are plotted alongside in figure 5(b). Clearly, $\bar {\Delta }(Z)$ is a non-monotone function that touches zero just once, with critical values in the range $Z^*\in [0,0.37]$ and $Z^*=0$ at $W=W^*$ (also note the corresponding values for the case of a circular source, shown as a dashed line).

Figure 5. (a) Discriminant function $\bar {\varDelta }$ (steady state), and (b) vertical position of the critical point $Z^*$, where the discriminant equals zero, in the case of a line source. The values of $Z^*$ calculated in the case of a circular source are shown in (b) with a dashed line.

By definition, at the critical point ($\bar {\Delta }(Z^*)=0$) the shape function $\bar {\alpha }$ and its derivative satisfy

(3.6)\begin{gather} \bar{\alpha}(Z^*)= 2\bar{P}_X(1,Z^*), \end{gather}

(3.7)\begin{gather}\bar{\beta}(Z^*)={-}\bar{\alpha}(Z^*)^{{-}1}. \end{gather}

Hence, inserting (3.6) and (3.7) into the linearised kinematic condition (3.5) at the critical point $Z=Z^*$ gives the ordinary differential equation

(3.8)\begin{equation} \bar{\alpha}^*\dot{\tilde{\alpha}}^*=\tilde{\alpha}^*-2\left.\tilde{P}_X^*\right|_{X=1},\end{equation}

where a dot indicates differentiation with respect to time, and $\ast$ superscripts indicate evaluation at the critical point (e.g. $\tilde {\alpha }^*(T)=\tilde {\alpha }(Z^*,T)$). The stability of the perturbation at the critical point therefore depends on the right-hand side of (3.8), which incidentally is proportional to the perturbed horizontal velocity at the edge of the plume ($\tilde {\alpha }^*-2\tilde {P}_X^*=2\bar {\alpha }^*\tilde {\mathcal {U}}^*|_{X=1}$), as shown in Appendix C. Hence (3.8) is rewritten as

(3.9)\begin{equation} \dot{\tilde{\alpha}}^*=2\left.\tilde{\mathcal{U}}^*\right|_{X=1}.\end{equation}

To assess the stability, we consider how the sign of the perturbed velocity $\tilde {\mathcal {U}}|_{X=1}$ relates to the sign of the perturbation $\tilde {\alpha }$ (i.e. whether the shape is deformed inwards or outwards at the critical point).

For positive perturbations $\tilde {\alpha }>0$, the expanded plume shape needs to be filled with fluid, so we expect a net positive velocity perturbation in the vicinity of the critical point, $\tilde {\mathcal {U}}|_{X=1}>0$. Likewise, for negative perturbations $\tilde {\alpha }<0$, a net negative velocity perturbation is required, $\tilde {\mathcal {U}}|_{X=1}<0$, such that the plume can shrink inwards. Hence $\tilde {\mathcal {U}}|_{X=1}$ is positively correlated with $\tilde {\alpha }$, indicating that (3.9) produces unstable solutions that grow unbounded with time when the perturbation is applied close to the critical point. Note, however, that when the perturbation (and therefore the flow) extends far away from the critical point, it is not obvious how $\tilde {\alpha }$ and $\tilde {\mathcal {U}}|_{X=1}$ are correlated.

Next, we demonstrate the existence of a small localised perturbation that grows unbounded over time. To do so we choose a simple Gaussian function for the initial perturbation, which is of the form

(3.10)\begin{equation} \tilde{\alpha}(Z^*,0)=\exp\left[ -{(Z-Z^*)^2}/{2\sigma^2} \right],\end{equation}

where $\sigma$ is the standard deviation. Whilst there are many possible local perturbations that cause instability to occur, we use this one since it is simple and demonstrates the point effectively.

We apply the perturbation (3.10) (using $\sigma =0.1$) to the case where $W=0.3$ (for which $Z^*=0.35$), and we plot the results in figure 6. Vector fields for the perturbed velocity, $(\tilde {\mathcal {U}},\tilde {\mathcal {W}})$ (see Appendix C for full expressions), are plotted in figures 6(a,b) for the cases of both a positive perturbation $\tilde {\alpha }>0$ and a negative perturbation $\tilde {\alpha }<0$. As shown by the velocity vectors, positive/negative perturbations result in a local flow outwards/inwards from the steady plume shape. Hence the perturbed horizontal velocity $\tilde {\mathcal {U}}|_{X=1}$, which is plotted in figure 6(c), is positive/negative (in the vicinity of the critical point) for positive/negative perturbations, indicating that the solution is unstable.

Figure 6. Vector fields for the perturbed velocity, $(\tilde {\mathcal {U}},\tilde {\mathcal {W}})$ (see Appendix C), in the cases of both (a) a negative perturbation $\tilde {\alpha }<0$, and (b) a positive perturbation $\tilde {\alpha }>0$. The stability is determined by the sign of the perturbed horizontal velocity at the edge of the plume, $\tilde {\mathcal {U}}|_{X=1}$, which is plotted in (c) for each case. The steady-state plume shape is indicated with a solid blue line in (a,b). The velocity ratio for this case is $W=0.3$, which has a critical point at $Z^*=0.35$.

Whilst this is just a specific case, it nevertheless demonstrates that an infinitesimal solution can be constructed that becomes unstable. Although we do not include the results here, we have also developed a time-dependent implicit numerical scheme that solves the perturbed equations for the pressure and plume shape at first order. These time-dependent numerical simulations confirm that localised perturbations of the form (3.10) become unstable when applied near the critical point. A more detailed eigenvalue analysis could explore the fastest growing perturbation (eigenfunction) and corresponding eigenvalue exponent. However, this lies outside the scope of the current study.

3.2. Moderate velocity ratios $W>W^* \wedge W\approx 1$

Now that we have addressed the case of small velocity ratios, for which the discriminant becomes zero at a critical point, we next address the case of moderate velocity ratios for which $\varDelta$ is always positive. We restrict our attention to velocity ratios that are larger than $W^*$ but which are still ${O}(1)$ (i.e. ignoring $W\gg 1$). Hence, taking $W\approx 1$, the pressure and the plume shape are well approximated by the expressions (2.22) and (2.23) (which are both ${O}(1-W)$). Meanwhile, the time-dependent kinematic condition (3.1) approximates to

(3.11)\begin{equation} -({k}/{\mu})\hat{p}_x\approx \phi a_t + w_b a_z,\quad x\approx a_0.\end{equation}

As before, we consider a small perturbation applied to the plume shape and pressure, but now we avoid converting to dimensionless coordinates for simplicity. Hence we consider a perturbation of the form

(3.12)\begin{gather} a=\bar{a}(z)+\epsilon\,\tilde{a}(z,t), \end{gather}

(3.13)\begin{gather}\hat{p}=\bar{p}(x,z)+\epsilon\,\tilde{p}(x,z,t). \end{gather}

Attention is required when performing this decomposition, since there are now two small parameters in the problem, namely $\epsilon >0$ and $\varepsilon =1-W>0$ (we consider $W<1$ without loss of generality). Hence in the following analysis, it is assumed that the perturbation to the shape is relatively much smaller than the perturbation to uniform flow, such that an asymptotic hierarchy $0<\epsilon \ll \varepsilon \ll 1$ is maintained.

After inserting (3.12) and (3.13) into (2.16)–(2.19), (2.21) and (3.11), and expanding in powers of $\epsilon$, it is clear that the leading-order expressions for the pressure and shape, $\bar {p}$ and $\bar {a}$, are precisely (2.22) and (2.23). Meanwhile, the unsteady terms satisfy

(3.14)\begin{equation} -({k}/{\mu})\tilde{p}_x=\phi\tilde{a}_t+w_b\tilde{a}_z,\quad x\approx a_0,\end{equation}

and the pressure perturbation $\tilde {p}$ satisfies homogeneous versions of the governing equations and boundary conditions (2.16)–(2.19), (2.21) (i.e. with zero on the right-hand sides of all the equations). Hence the pressure perturbation is trivial, $\tilde {p}=0$, and consequently (3.14) also becomes homogeneous. In this way, perturbations to the plume shape $\tilde {a}$ are simply advected downstream with dimensional velocity $w_b/\phi$, and consequently the system is stable.

3.3. A note on the case of large velocity ratios $W>1$

It is worth discussing briefly the plume behaviour in the case where $W$ is larger than 1, since this has so far been neglected. As discussed earlier, velocity ratios $W>1$ are associated with an expanding plume shape. This expansion requires large horizontal velocities, and hence large values of $-P_X|_{X=1}$, suggesting that the discriminant function (2.32) may become zero or negative for large values of $W$.

It is found that the numerical solution of the steady plume shape for $W>1$ has a sharp gradient near the origin, followed by a slow tapering off. This sharp expansion causes the discriminant to become negative ($\varDelta <0$) near the origin for velocity ratios larger than $W_n\approx 1.14$ in the case of a line source, and $W_n\approx 1.19$ in the case of a circular source. For $W$ larger than these values, a numerical solution satisfying the steady kinematic condition (2.29) cannot be found. This suggests that a steady plume shape may not be possible, or at least may not take the form of a smooth continuous curve $\alpha (Z)$, as we have prescribed. For example, it is possible that the shape may jump discontinuously from $\alpha (0)=1$ to a larger value near the origin. However, the current modelling approach cannot be used to investigate such scenarios since we require a smooth continuous transformation of the form (2.24a–d).

Whilst further analytical treatment of this case lies outside the scope of this study, it is worth commenting on the validity of the modelling approach for $W\gg 1$. In such cases, large velocities near the origin may violate the assumptions of Darcy's law. In particular, the Reynolds number of the flow near the origin is given by

(3.15)\begin{equation} Re=\frac{\rho_1 w_0 d_p}{\mu}=\frac{W \rho_1 w_b d_p}{\mu}, \end{equation}

where $d_p$ is the pore size. Since $Re \propto W$, it is expected that inertial effects may become important for $W\gg 1$ (e.g. see Sahu & Flynn Reference Sahu and Flynn2015). Likewise, the Péclet number is also proportional to $W$, indicating that mixing due to dispersion cannot be neglected for $W\gg 1$. Hence, to study the shape and stability of the plume in such scenarios, a more sophisticated model that accounts for laminar flow and dispersive mixing might be more appropriate than the linear stability analysis based on Darcy's law used here.

3.4. Further cases and considerations

We have therefore shown that plume shapes are stable in the regime $W\approx 1$, and unstable in the regime $W_0< W< W^*$. It is not known whether the plume shapes for $W>W^*$ but $W\not \approx 1$ are stable or unstable (e.g. one could argue that $W=0.6\not \approx 1$), and for this a full eigenvalue analysis is required. Such an analysis could predict the largest value of $W$ that onsets instability (i.e. producing a positive real-valued eigenvalue). However, for the purposes of this study, and since our experimental data suggest a stable shape for velocity ratios as small as $W=0.58$ (see figure 3a), this analysis covers the relevant and interesting range of cases.

To go beyond the current analysis and resolve the nonlinear instability at very small velocity ratios $W< W_0$, a time-dependent numerical simulation is required in either two or three dimensions (e.g. see Hewitt et al. Reference Hewitt, Neufeld and Lister2013b; Hewitt & Lister Reference Hewitt and Lister2017). In particular, it is not clear exactly how the flow evolves over time, whether it forms fingers, filaments or disconnected regions. Moreover, in the presence of dispersion, it is likely that thin disconnected regions may fuse together. It is possible that the stability diagram in figures 2(b) and 4(b) may not be complete, but in fact could contain other stable, unstable or periodic branches. Hence such a numerical simulation could be used to explore the full range of solution branches, and classify any bifurcation points, such as at $W=W_0$, that may be a fold (saddle node) or a Hopf bifurcation.

As described earlier, it is worth commenting briefly on the effects of the ambient fluid in the unsteady case. In particular, in the above analysis we have assumed implicitly (i.e. via the dynamic boundary condition (2.10)) that the injected fluid and the ambient fluid are decoupled. However, in the unsteady case (unlike the steady case) the ambient fluid is not quiescent. In fact, as the plume expands or contracts over time, the ambient fluid must be displaced accordingly. However, as long as the perturbation to the plume shape has a sufficiently large aspect ratio, the rearrangement of the ambient fluid should not affect the plume growth to good approximation. This is similar to the shallow-water approximation used for gravity currents in porous media (Huppert & Woods Reference Huppert and Woods1995), for which the injected fluid is not coupled to the ambient fluid so long as the aspect ratio of the current remains large.

This argument does not extend to the case where the viscosity of the injected fluid is different to that of the ambient fluid. In this case, the viscosity contrast affects the pressure gradients on either side of the perturbation to the interface as it grows. If $M=\mu _1/\mu _2$ denotes the viscosity contrast between injected and ambient fluids, then the interface becomes destabilised for $M<1$ and stabilised for $M>1$, as described in the famous study by Saffman & Taylor (Reference Saffman and Taylor1958). For this reason, we do not expect the results of § 3 to apply in the case where $M\neq 1$.