2. Model

The set-up is shown in figure 1. A similar geometry was studied by Nordbotten & Celia (Reference Nordbotten and Celia2006) and Guo et al. (Reference Guo, Zheng, Celia and Stone2016). Fluid of density $\rho$ and viscosity $\mu _i$ is injected into a porous medium of thickness $H_0$, horizontal permeability $k_x$, and porosity $\phi$. The medium is initially filled with a second fluid of greater density $\rho + \Delta \rho$ and different viscosity $\mu _a$. The top and bottom boundaries of the layer are impermeable. The input fluid is injected with volume flux $Q$ from a point source at the upper boundary, which we take to be the origin of our coordinate system. The $Z$ coordinate is measured downwards from the upper boundary at which $Z=0$. We use cylindrical polar coordinates, with $R$ denoting the distance from the $Z$ axis, and $\theta$ denoting the angular coordinate. Time is denoted by $T$. We neglect any intermingling of the fluids and assume that there is a sharp interface at $Z=H(R,\theta,T)$ between the fluids. We also assume that the two fluids are segregated by buoyancy so that $H$ is a single-valued function.

When the input fluid has spread far from the source, the flow is predominantly in the radial direction. We may then apply the shallow flow approximation and the pressure is approximately hydrostatic,

(2.1)\begin{gather} P= \bar{P} + \rho g Z \quad \text{for} \quad 0< Z< H, \end{gather}

(2.2)\begin{gather}P= \bar{P} -\Delta \rho\,g H + (\rho+\Delta \rho) g Z \quad \text{for} \quad H< Z< H_0, \end{gather}

in the input and ambient fluids, respectively, where $\bar {P}(R,\theta,T)$ is the pressure at the upper boundary, $Z=0$, which is to be determined. The Darcy velocities are given by

(2.3)\begin{gather} \boldsymbol{U}_i ={-}\frac{k_x}{\mu_i}\,\boldsymbol{\nabla} \bar{P}, \end{gather}

(2.4)\begin{gather}\boldsymbol{U}_a ={-}\frac{k_x}{\mu_a}\left( \boldsymbol{\nabla} \bar{P} - \Delta \rho\,g\,\boldsymbol{\nabla} H \right), \end{gather}

in the input and ambient fluids, respectively. The governing equations of the flow are

(2.5)\begin{equation} \phi\,\frac{\partial H}{\partial T} + \boldsymbol{\nabla}\boldsymbol{\cdot} \left( H \boldsymbol{U}_i \right) = 0 \end{equation}

and

(2.6)\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot} \left[ H \boldsymbol{U}_i + (H_0-H) \boldsymbol{U}_a \right] = 0, \end{equation}

which represent mass conservation integrated over the thickness of the input fluid, and mass conservation integrated over the thickness of the layer, respectively.

At $T=0$, the porous layer is filled entirely with the ambient fluid, and constant-flux injection of the input fluid begins. Global mass conservation of the input fluid takes the form

(2.7)\begin{equation} \int_0^{2{\rm \pi}} \int_0^{R_f(\theta,T)} \phi R H \, \mathrm{d} R\, \mathrm{d} \theta = QT, \end{equation}

where $R_f(\theta,T)$ denotes the frontal contact line where $H=0$ (see figure 1b).

Eventually, the fluid–fluid interface will touch the lower boundary (Nordbotten & Celia Reference Nordbotten and Celia2006). It is well-known that the predicted interface shape of an unconfined axisymmetric gravity current has an unphysical singularity at the origin, whilst in confined layers, the flow will always eventually encompass the entire layer near the source, regardless of the thickness of the layer (Lyle et al. Reference Lyle, Huppert, Hallworth, Bickle and Chadwick2005; Guo et al. Reference Guo, Zheng, Celia and Stone2016). Recently, Benham, Neufeld & Woods (Reference Benham, Neufeld and Woods2022) carried out a detailed analysis of the flow near the source, incorporating the pressure-driven vertical flow there. Their results corroborate the conclusion that a confined layer becomes filled with input fluid at sufficiently late times (see also Huppert & Pegler Reference Huppert and Pegler2022). After the fluid–fluid interface has touched the lower boundary, there is a second contact line at $R=R_b(\theta,T)$ where $H=H_0$. The flow then has three regions: there is the ‘fully flooded’ zone in which the input fluid fills the layer ($H=H_0$, $0\leqslant R < R_b(\theta,T)$); next is the partially filled zone in which the fluid–fluid interface lies between the top and bottom boundaries ($0< H< H_0$, $R_b(\theta,T)< R< R_f(\theta,T)$); finally, beyond this there is an uninvaded zone in which the ambient fluid fills the layer ($H=0$, $R > R_f(\theta,T)$); see figure 1(b).

When the input fluid has filled the layer near the source, the boundary conditions are as follows. At the source,

(2.8a,b)\begin{equation} H \to H_0, \quad \frac{\partial \bar{P}}{\partial R} \to \frac{-\mu_i Q}{2 {\rm \pi}k_x H_0 R} \quad \text{as} \ R \to 0, \end{equation}

where the latter arises from volume conservation. In the far field, the boundary condition is

(2.9a,b)\begin{equation} H \to 0, \quad \frac{\partial \bar{P}}{\partial R} \to \frac{-\mu_a Q}{2 {\rm \pi}k_x H_0 R} \quad \text{as} \ R \to \infty. \end{equation}

The dependent variables, $\bar {P}(R,\theta,T)$ and $H(R,\theta,T)$, are defined on $\theta \in [0,2 {\rm \pi})$, with periodic boundary conditions $\bar {P}(R,0,T)=\bar {P}(R,2 {\rm \pi},T)$ and $H(R,0,T)=H(R,2 {\rm \pi},T)$. We also require that $0 \leqslant H \leqslant H_0$. Equations (2.5) and (2.6) combined with global mass conservation (2.7) and the boundary conditions (2.8a,b), (2.9a,b) define a complete coupled system for determining the two dependent variables $\bar {P}(R,\theta,T)$ and $H(R,\theta,T)$.

Although the system derived above is sufficient to solve the problem numerically (see § 3), the stability investigation requires detailed treatment of the behaviour in the vicinity of the two contact lines where the interface touches the layer boundaries. Across these lines, the gradient of the interface and the gradient of the pressure at the top of the layer can be discontinuous. Indeed, it is well-known that the stability of the fluid–fluid interface for equally dense fluids is controlled by the jump in the pressure gradient at the interface (see § 1). Therefore, we derive boundary conditions for the interface shape and the pressure gradient across each contact line.

At the trailing contact line $R=R_b(\theta, T)$, continuity of the interface, continuity of the pressure and continuity of the flux of the input fluid take the forms

(2.10a–c)\begin{equation} H(R=R_b^{+})=H(R=R_b^{-})=H_0, \quad \left[\bar{P}\right]^{+}_- = 0, \quad \left[ \frac{\partial \bar{P}}{\partial R} - \frac{\partial R_b}{\partial \theta}\,\frac{1}{R}\, \frac{\partial \bar{P}}{\partial \theta} \right]^{+}_- = 0, \end{equation}

respectively, where $[\cdot ]^{+}_-$ denotes the discontinuity in the quantity in square brackets at the contact line $R=R_b(\theta, T)$. In addition, we have a kinematic boundary condition; at the trailing contact line, the normal velocity of the ambient fluid is equal to the normal velocity of the interface, which takes the form

(2.11)\begin{gather} -\frac{k_x}{\mu_a} \left.\left(\frac{\partial \bar{P}}{\partial R} - \Delta \rho\,g\, \frac{\partial H}{\partial R} \right)\right\vert_{R=R_b^{+}} = \left.\frac{\partial R_b}{\partial T} -\frac{k_x}{\mu_a}\,\frac{\partial R_b}{\partial \theta} \left(\frac{1}{R}\,\frac{\partial \bar{P}}{\partial \theta} - \frac{\Delta \rho\,g}{R}\,\frac{\partial H}{\partial \theta} \right)\right\vert_{R=R_b^{+}}. \end{gather}

At the leading contact line $R=R_f(\theta,T)$, the corresponding boundary conditions are

(2.12a,b)\begin{gather} H(R=R_f^{+})=H(R=R_f^{-})=0, \quad \left[\bar{P}\right]^{+}_- = 0, \end{gather}

(2.13)\begin{gather} \left[ \frac{\partial (\bar{P}-\Delta \rho\,g H)}{\partial R} - \frac{\partial R_f}{\partial \theta}\,\frac{1}{R}\,\frac{\partial (\bar{P}-\Delta \rho\,g H)}{\partial \theta} \right]^{+}_- = 0, \end{gather}

(2.14)\begin{gather} \left.-\frac{k_x}{\mu_i}\,\frac{\partial \bar{P}}{\partial R}\right\vert_{R=R_f^{-}} = \left.\frac{\partial R_f}{\partial T} -\frac{k_x}{\mu_i R}\,\frac{\partial R_f}{\partial \theta}\, \frac{\partial \bar{P}}{\partial \theta}\right\vert_{R=R_f^{-}}, \end{gather}

where $[\cdot ]^{+}_-$ denotes the discontinuity in the quantity in square brackets at the contact line $R=R_f(\theta, T)$. Equation (2.14) relates the normal velocity of the interface to the normal velocity of the input fluid just behind the leading contact line (at $R=R_f^{-}$).

2.1. Transformed coordinate system

Mass conservation of the input fluid suggests that the radial extent of the interface grows with $R^{2} \sim QT/(\phi H_0)$. This motivates introducing the transformed dimensionless coordinates

(2.15a,b)\begin{equation} \eta = \left(\frac{2 {\rm \pi}\phi H_0}{Q T}\right)^{1/2} R, \quad {\tau=\log \left(T/T_{{ref}} \right)}, \end{equation}

where $T=T_{{ref}}>0$ is a reference time, which corresponds to $\tau =0$. We also scale the dependent variables to obtain the dimensionless quantities

(2.16a,b)\begin{equation} \mathcal{H}(\eta,\theta,\tau) = \frac{H(R,\theta,T)}{H_0}, \quad \mathcal{P}(\eta,\theta,\tau)=\frac{2 {\rm \pi}k_x H_0}{\mu_i Q}\, \bar{P}(R,\theta,T). \end{equation}

The two mass conservation equations (2.5) and (2.6) become

(2.17)\begin{gather} \frac{\partial \mathcal{H}}{\partial \tau}-\frac{\eta}{2}\,\frac{\partial \mathcal{H}}{\partial \eta} - \tilde{\boldsymbol{\nabla}} \boldsymbol{\cdot} \left( \mathcal{H}\,\tilde{\boldsymbol{\nabla}} \mathcal{P} \right) = 0, \end{gather}

(2.18)\begin{gather}\tilde{\boldsymbol{\nabla}} \boldsymbol{\cdot} \left[\left(M+(1-M)\mathcal{H}\right) \tilde{\boldsymbol{\nabla}} \mathcal{P} - G M (1-\mathcal{H})\, \tilde{\boldsymbol{\nabla}} \mathcal{H} \right] = 0, \end{gather}

where

(2.19)\begin{equation} \tilde{\boldsymbol{\nabla}}=\left(\frac{\partial}{\partial (\eta \cos \theta)}, \, \frac{\partial}{\partial (\eta \sin \theta)} \right) \end{equation}

is the gradient operator with respect to the transformed coordinates, and we have introduced the two dimensionless groups

(2.20a,b)\begin{equation} M=\frac{\mu_i}{\mu_a}, \quad G= \frac{2 {\rm \pi}k_x\,\Delta \rho\,g H_0^{2}}{\mu_i Q}, \end{equation}

which are the viscosity ratio and gravity number, respectively. The latter represents the magnitude of pressure gradients associated with buoyancy relative to pressure gradients associated with injection of fluid. In the case of equally dense fluids, the gravity number is $G=0$. The boundary conditions at the source and in the far field are

(2.21a,b)\begin{equation} \mathcal{H} \to 1, \quad \frac{\partial \mathcal{P}}{\partial \eta} \to \frac{-1}{\eta} \quad \text{as} \ \eta \to 0, \end{equation}

and

(2.22a,b)\begin{equation} \mathcal{H} \to 0, \quad \frac{\partial \mathcal{P}}{\partial \eta} \to \frac{-M^{{-}1}}{\eta} \quad \text{as} \ \eta \to \infty. \end{equation}

The trailing contact line in the transformed coordinates is denoted by $\eta =\eta _b(\theta,\tau )$, and the boundary conditions there ((2.10a–c) and (2.11)) are recast as

(2.23a–c)\begin{gather} \mathcal{H}(\eta=\eta_b^{+})=\mathcal{H}(\eta=\eta_b^{-})=1, \quad \left[\mathcal{P} \right]^{+}_- =0, \quad \left[\frac{\partial \mathcal{P}}{\partial \eta} - \frac{\partial \eta_b}{\partial \theta}\, \frac{1}{\eta}\,\frac{\partial \mathcal{P}}{\partial \theta} \right]^{+}_- =0, \end{gather}

(2.24)\begin{gather} -M\left.\left( \frac{\partial \mathcal{P}}{\partial \eta} - G\,\frac{\partial \mathcal{H}}{\partial \eta} \right)\right\vert_{\eta=\eta_b^{+}}=\left.\frac{\partial \eta_b}{\partial \tau} +\frac{1}{2}\,\eta_b -M\,\frac{\partial \eta_b}{\partial \theta} \left( \frac{1}{\eta}\, \frac{\partial \mathcal{P}}{\partial \theta} -\frac{G}{\eta}\,\frac{\partial \mathcal{H}}{\partial \theta}\right)\right\vert_{\eta=\eta_b^{+}}, \end{gather}

whilst at the leading contact line, $\eta =\eta _f(\theta,\tau )$,

(2.25a,b)\begin{gather} \mathcal{H}(\eta=\eta_f^{+})=\mathcal{H}(\eta=\eta_f^{-})=0, \quad \left[\mathcal{P} \right]^{+}_- =0, \end{gather}

(2.26)\begin{gather}\left[ \frac{\partial (\mathcal{P}-G \mathcal{H})}{\partial \eta} - \frac{\partial \eta_f}{\partial \theta}\, \frac{1}{\eta}\,\frac{\partial (\mathcal{P}-G\mathcal{H})}{\partial \theta} \right]^{+}_- =0, \end{gather}

(2.27)\begin{gather}\left.-\frac{\partial \mathcal{P}}{\partial \eta}\right\vert_{\eta=\eta_f^{-}} =\left.\frac{\partial \eta_f}{\partial \tau} +\frac{1}{2}\,\eta_f -\frac{1}{\eta}\,\frac{\partial \eta_b}{\partial \theta}\,\frac{\partial \mathcal{P}}{\partial \theta}\right\vert_{\eta=\eta_f^{-}}. \end{gather}

Mass conservation (2.7) takes the form

(2.28)\begin{equation} {\frac{1}{2 {\rm \pi}} \int_0^{2 {\rm \pi}} \int_0^{\eta_f} \eta \mathcal{H} \, \mathrm{d} \eta \, \mathrm{d} \theta = 1}. \end{equation}

2.2. Self-similar axisymmetric solution

The flow has a self-similar axisymmetric solution with (Nordbotten & Celia Reference Nordbotten and Celia2006; Guo et al. Reference Guo, Zheng, Celia and Stone2016)

(2.29a,b)\begin{equation} \mathcal{H}(\eta,\theta,\tau)=\mathcal{H}_0(\eta), \quad \mathcal{P}(\eta,\theta,\tau) = \mathcal{P}_0(\eta), \end{equation}

and

(2.30a,b)\begin{equation} \eta_b(\theta,\tau) = {\eta_{b_0}}, \quad \eta_f(\theta,\tau) = {\eta_{f_0}}, \end{equation}

so that the contact lines are circles, fixed in similarity space. We are interested in analysing the stability of and convergence to these solutions in §§ 3, 4 and 5. In the present section, we derive a single ordinary differential equation that governs the axisymmetric self-similar interface shape, and we recall an important analytical solution of this equation from Nordbotten & Celia (Reference Nordbotten and Celia2006) and Guo et al. (Reference Guo, Zheng, Celia and Stone2016). For axisymmetric flow, the pressure at the top boundary, $\mathcal {P}_0(\eta )$, can be eliminated from the problem as follows. We integrate (2.18) with respect to $\eta$, and apply the source flux boundary condition (2.21a,b) to obtain

(2.31)\begin{equation} \frac{\mathrm{d} \mathcal{P}_0}{\mathrm{d} \eta} = \frac{-1}{\eta\left(M+(1-M)\mathcal{H}_0\right)} + \frac{GM (1-\mathcal{H}_0)}{\left(M+(1-M)\mathcal{H}_0\right)}\,\frac{\mathrm{d} \mathcal{H}_0}{\mathrm{d} \eta}. \end{equation}

Upon substituting (2.31) into (2.17), we obtain the following ordinary differential equation for $\mathcal {H}_0(\eta )$:

(2.32)\begin{equation} -\frac{\eta}{2}\,\frac{\mathrm{d} \mathcal{H}_0}{\mathrm{d} \eta} + \frac{1}{\eta}\, \frac{\mathrm{d}}{\mathrm{d} \eta} \left(\frac{\mathcal{H}_0}{M+(1-M) \mathcal{H}_0} \right) = \frac{1}{\eta}\,\frac{\mathrm{d}}{\mathrm{d} \eta} \left(\eta\, \frac{G M \mathcal{H}_0(1-\mathcal{H}_0)}{M+(1-M) \mathcal{H}_0}\,\frac{\mathrm{d} \mathcal{H}_0}{\mathrm{d} \eta} \right). \end{equation}

This axisymmetric system can be solved numerically for $\mathcal {H}_0(\eta )$ as described in Guo et al. (Reference Guo, Zheng, Celia and Stone2016) for a wide range of $M$ and $G$. For $G\gg 1$ (buoyancy-dominated flow), the input fluid forms a thin current near the top of the layer (except in a small neighbourhood of the source, where the layer is fully flooded). In this case, the displacement of the ambient fluid is unimportant over most of the radial extent of the flow. The flow is effectively ‘unconfined’ and single-phase. Hence we do not expect the large $G$ solution to be unstable to viscous fingering.

Instead, we focus on the case of relatively larger input flux ($G \leqslant {O}(1$)) and a less viscous input fluid ($M<1$) for which stability is not well understood. For $G \ll 1$ and $M<1$, (2.32) has an approximate analytic solution. We neglect the right-hand side of (2.32). Although this removes the highest derivative of $\mathcal {H}_0$, the boundary conditions (2.24) and (2.27) become independent of $\mathcal {H}_0$ for $G\ll 1$ and simply give the locations of the contact points, so the problem remains complete. Integrating the left-hand side of (2.32) furnishes the interface shape (Guo et al. Reference Guo, Zheng, Celia and Stone2016)

(2.33)\begin{gather} \mathcal{H}_0=1, \quad 0<\eta< {\eta_{b_0}}, \end{gather}

(2.34)\begin{gather}\mathcal{H}_0 = \frac{(2M)^{1/2}}{(1-M)\eta} -\frac{M}{1-M}, \quad {\eta_{b_0}}<\eta< {\eta_{f_0}}, \end{gather}

(2.35)\begin{gather}\mathcal{H}_0=0, \quad {\eta_{f_0}}<\eta, \end{gather}

where the contact lines are given by

(2.36a,b)\begin{equation} {\eta_{b_0}=\sqrt{2M}} \quad \text{and} \quad {\eta_{f_0} =\sqrt{2/M}.} \end{equation}

Note that these provide inner bounds on the locations of the contact lines for general $G> 0$ provided that $M<1$:

(2.37a,b)\begin{equation} {\eta_{b_0}} \leqslant \sqrt{2M}, \quad {\eta_{f_0}} \geqslant \sqrt{2/M}, \end{equation}

with equality as $G \to 0^{+}$. These inequalities arise because buoyancy acts to extend the interface (its effect is proportional to $G$).

The corresponding pressure at the top of the layer takes the form

(2.38)\begin{gather} \mathcal{P}_0={-}1-\log\left(\frac{\eta}{(2M)^{1/2}}\right), \quad 0<\eta< {\eta_{b_0}}, \end{gather}

(2.39)\begin{gather}\mathcal{P}_0={-}\frac{\eta}{(2M)^{1/2}}, \quad {\eta_{b_0}}<\eta< {\eta_{f_0}}, \end{gather}

(2.40)\begin{gather}\mathcal{P}_0={-}M^{{-}1} -M^{{-}1} \log\left(\frac{\eta}{(2/M)^{1/2}}\right), \quad {\eta_{f_0}}<\eta, \end{gather}

where $\mathcal {P}_0(\eta )$ is defined up to addition of an arbitrary constant.

The numerical solution of (2.32) for $\mathcal {H}_0$ is plotted as a blue line in figure 2(a) for $G=0.1$ and $M=0.2$. The small-$G$ analytical solution (2.34) is shown as a black dashed line, and there is good agreement with the numerical result. Figure 2(b) shows the analytical solutions for various values of $M<1$. When the input fluid is of relatively lower viscosity (smaller $M$), it forms a finger that intrudes further into the ambient fluid at the top of the layer.

Figure 2. Self-similar axisymmetric flow. (a) Interface shape for viscosity ratio $M=0.2$ and $G=0.1$ (blue line) and the small-$G$ approximate solution ((2.34), black dashed line). Stability is controlled by the pressure gradients at points (i)–(iv); see § 4. (b) Small-$G$ approximate interface shapes (2.34) for four values of $M$.

It is worth considering the behaviour in the limit $G \to 0$, corresponding to equally dense fluids and no buoyancy force. As $G \to 0$, the axisymmetric analytic solution (2.34) becomes a more accurate approximation of the numerical solution, with the input fluid preferentially flowing near the upper boundary (figure 2). However, $G=0$ corresponds to no buoyancy segregation ($\Delta \rho =0$), so the input fluid should have no preference for the top of the layer. Indeed, in the case of equally dense fluids ($G=0$) and $M<1$, one would observe classical Saffman–Taylor fingering (Paterson Reference Paterson1981). The derivation of the axisymmetric analytical solution (2.34) assumes implicitly that the input and ambient fluid have been segregated by buoyancy, which requires $G>0$ and takes a time proportional to $G^{-1}$. Hence the case $G=0$ has qualitatively different late-time behaviour to the case of small but non-zero $G$, and the limit $G \to 0$ is singular. In this paper, we analyse the stability when $G$ is small but positive; we exclude the case $G=0$, for which buoyancy has no effect.

For $G$ small but non-zero, the interfaces in figure 2(b) occur provided that (i) buoyancy segregation has occurred and (ii) the radial extent of the flow is much greater than the layer thickness so that the shallow approximation applies. For $M<1$, in dimensional terms, the former requires

(2.41)\begin{equation} T \gg \frac{\mu_a H_0}{\Delta \rho\,g k_z}, \end{equation}

where $k_z$ is the vertical permeability in the layer. Equation (2.41) corresponds to the time for vertical flow across the layer. The shallow approximation requires

(2.42)\begin{equation} T \gg \frac{2 {\rm \pi}\phi H_0^{3}}{Q}. \end{equation}

Prior to buoyancy segregation, Saffman–Taylor fingering can occur. The present paper is concerned with the post-segregation stability of the axisymmetric solutions of (2.32) to angular-dependent perturbations with $\Delta \rho >0$. For $G\ll 1$, the first condition (2.41) implies the second (2.42). Notice that for small density differences $\Delta \rho$, the time for buoyancy segregation is very large. It is also important to note that the segregation timescale (2.41) depends only on the vertical permeability $k_z$, whereas the gravity number (and hence the self-similar solutions) depends only on horizontal permeability $k_x$ (cf. Benham et al. Reference Benham, Neufeld and Woods2022). In an anisotropic layer, the horizontal permeability $k_x$ is different to the vertical permeability $k_z$. For an isotropic medium, $k_x=k_z$. The results that we derive in this paper concerning flow stability extend to porous layers that have different vertical and horizontal permeabilities (as is common in many aquifers). The case of more complicated anisotropy, such as different horizontal permeabilities in different directions, is beyond the scope of the present work. (Vertically varying permeability is analysed in § 6.)

For simplicity, we have modelled the case where the source of input fluid is located at a point on the upper boundary. However, the analysis and results of this paper apply for any location of the injection source (for example, at the lower boundary, or over a vertical line within the layer). This is because the input fluid will always eventually fill the thickness of the layer near the source. Subsequently, the flow becomes predominantly radial. There will, of course, be different early-time transient behaviour for different source locations, but once (2.41) and (2.42) apply, the exact source location becomes unimportant.

3. Numerical results

In order to investigate the response of the flow to perturbations with angular dependence, we develop a numerical method to solve the coupled system for the evolution of the fluid–fluid interface $\mathcal {H}(\eta, \theta, \tau )$, and the pressure at the upper boundary $\mathcal {P}(\eta, \theta, \tau )$. This consists of (2.17) and (2.18) with boundary conditions (2.21a,b) and (2.22a,b), and an appropriate initial condition. We solve this system numerically using a finite-difference scheme, details of which are given in Appendix A.

The numerical method is applied to a wide range of axisymmetric and non-axisymmetric initial conditions over many values of the parameters $G$ and $M$. In all cases, the solution converges to the axisymmetric similarity solution. This demonstrates both the veracity of the numerical method and that the similarity solution provides the intermediate asymptotics for many initial conditions. Examples are shown in figures 3 and 4 illustrating how mode six and mode three fingers are suppressed even though the input fluid is of lesser viscosity than the ambient. We find that any perturbation is suppressed, including radial perturbations, provided that the fluids remain buoyancy segregated. If the interface becomes non-monotonic as a function of vertical coordinate, then the model breaks down. Stability in this case is not well understood and is beyond the scope of the present work.

Figure 3. Suppression of mode six fingers for $G=0.1$ and $M=0.2$. (a) Initial condition at the reference time ($\tau =0$) for the flow thickness $\mathcal {H}$ in transformed coordinates. (b,c) The flow thicknesses at $\tau =1.25$ and $\tau =4$. (d) The evolution of the contact line, where $\mathcal {H}=0$.

Figure 4. Suppression of mode three fingers for $G=0.5$ and $M=0.5$. (a) Initial condition ($\tau =0$) for the interface shape. (b) Initial condition ($\tau =0$) for the pressure gradient (log of its magnitude is shown). (c,d) Corresponding plots at $\tau =1$. Radial cross-sections of these results are shown in figure 5.

Figure 5. Suppression of mode three fingers for $G=0.5$ and $M=0.5$ (same parameters and initial condition as figure 4). Radial cross-sections at $\theta =0$ (black lines) and $\theta ={\rm \pi} /3$ (blue lines) of the flow thickness $\mathcal {H}$, and radial pressure gradient at the upper boundary $\partial \mathcal {P}/\partial \eta$, at three times: (a,d) $\tau =0.25$, (b,e) $\tau =1$, (c,f) $\tau =4$. The red dashed lines denote the locations of the frontal contact line $\eta _f$ at $\theta =0$ and $\theta ={\rm \pi} /3$. The discontinuity in the pressure gradient at the frontal contact line is slightly smoothed at earlier times owing to non-axisymmetry.

The numerical results also demonstrate that there is a jump in the radial pressure gradient at the upper boundary across the leading contact line (see figure 5). This pressure gradient jump is associated with the density difference between the two fluids, and the discontinuity vanishes as $G \to 0^{+}$. The discontinuity has a stabilising effect because the magnitude of the pressure gradient is smaller in the ambient fluid (see § 4). However, this discontinuity in the hydrostatic pressure is not required for the interface to be stable. Indeed, in § 5 we show that the interface is stable even in the limit $G \to 0^{+}$.

4. Stability mechanism

In this section, we describe the mechanism by which buoyancy segregation suppresses viscous fingering. In the classical case of radial flow with no cross-layer gravity (discussed in § 1), the pressure gradient either side of the fluid–fluid interface is proportional to the fluid viscosity owing to continuity of the radial flux, which drives the well-known instability.

We now analyse the more complicated three-dimensional case in which the fluids have different densities and are segregated by buoyancy. Between the contact lines, the pressure gradient in the input fluid is $\mathrm {d} \mathcal {P}_0/\mathrm {d} \eta$, and the pressure gradient in the ambient fluid is $\mathrm {d} (\mathcal {P}_0-G\mathcal {H}_0)/\mathrm {d} \eta$, where the $-G \mathcal {H}_0$ contribution arises from the density difference between the fluids. Since both $\mathrm {d} \mathcal {P}_0/\mathrm {d} \eta$ and $\mathrm {d} \mathcal {H}_0/\mathrm {d} \eta$ are negative, the magnitude of the pressure gradient is larger in the input fluid (in contrast to the classical case, which is independent of gravity), which suggests that between the contact lines, the interface is stable for $G>0$. The stability at the contact lines must be treated separately because there is a discontinuity in $\mathrm {d} \mathcal {P}_0/\mathrm {d} \eta$ and $\mathrm {d} \mathcal {H}_0/\mathrm {d} \eta$ between the input and ambient fluids, as shown in figure 5. We consider the driving pressure gradients in the input and ambient fluids either side of the contact lines; these locations are labelled (i)–(iv) in figure 2a.

At location (i), just ahead of the leading contact line at $\eta ={\eta _{f_0}}^{+}$, the interface gradient and pressure gradient at the top boundary are

(4.1a,b)\begin{equation} \frac{\mathrm{d} \mathcal{H}_0}{\mathrm{d} \eta} =0, \quad \frac{\mathrm{d} \mathcal{P}_0}{\mathrm{d} \eta} ={-}\frac{1}{M {\eta_{f_0}}}. \end{equation}

The latter is calculated from (2.31). The pressure gradient driving the ambient fluid at (i) is

(4.2)\begin{equation} \frac{\mathrm{d} (\mathcal{P}_0-G \mathcal{H}_0)}{\mathrm{d} \eta}={-}\frac{1}{M {\eta_{f_0}}}. \end{equation}

At location (ii), just behind the leading contact line at $\eta ={\eta _{f_0}}^{-}$, the interface and pressure gradients are

(4.3a,b)\begin{equation} \frac{\mathrm{d} \mathcal{H}_0}{\mathrm{d} \eta} =\frac{1}{{\eta_{f_0}} M G} -\frac{{\eta_{f_0}}}{2G}, \quad \frac{\mathrm{d} \mathcal{P}_0}{\mathrm{d} \eta} ={-}\frac{{\eta_{f_0}}}{2}. \end{equation}

The former is calculated by taking $\mathcal {H} \to 0$ in (2.32). The pressure gradient driving the input fluid at (ii) is given by (4.3b), and its magnitude is greater than (or equal to) the magnitude of the pressure gradient in the ambient fluid (4.2) because ${\eta _{f_0}} \geqslant \sqrt {2/M}$ (see (2.37a,b)). Hence small perturbations to the contact line decay because they experience unfavourable pressure gradients, which suggests interfacial stability provided that $G>0$. The driving pressure gradients either side of the contact line are equal in the limit $G \to 0^{+}$; this case is discussed in more detail in § 5.

Inspection of equations (4.1b) and (4.3b) demonstrates that there is a discontinuity in the pressure gradient at the top boundary, $\mathrm {d} \mathcal {P}_0/\mathrm {d} \eta$, across $\eta ={\eta _{f_0}}$, which can be observed in the numerical results in figure 5. This discontinuity vanishes as $G \to 0$ even for fluids with different viscosities. In other words, the destabilising jump in the pressure gradient associated with the viscosity contrast is eliminated by the radial extension of the interface. Indeed, the jump in the pressure gradient is stabilising for $G>0$.

Similar analysis applies at the trailing contact line. At location (iii), just downstream of the trailing contact line $\eta ={\eta _{b_0}}^{+}$, the interface gradient, the pressure gradient at the upper boundary and the pressure gradient driving the ambient fluid are given by

(4.4a–c)\begin{equation} \frac{\mathrm{d} \mathcal{H}_0}{\mathrm{d} \eta} =\frac{{\eta_{b_0}}}{2 M G} - \frac{1}{G {\eta_{b_0}}}, \quad \frac{\mathrm{d} \mathcal{P}_0}{\mathrm{d} \eta} ={-}\frac{1}{{\eta_{b_0}}}, \quad \frac{\mathrm{d} (\mathcal{P}_0-G \mathcal{H}_0)}{\mathrm{d} \eta}={-} \frac{{\eta_{b_0}}}{2M}, \end{equation}

respectively. At location (iv), just behind the trailing contact line $\eta ={\eta _{b_0}}^{-}$, the interface gradient and the pressure gradient at the upper boundary (which drives the input fluid) are given by

(4.5a,b)\begin{equation} \frac{\mathrm{d} \mathcal{H}_0}{\mathrm{d} \eta} =0, \quad \frac{\mathrm{d} \mathcal{P}_0}{\mathrm{d} \eta} ={-}\frac{1}{{\eta_{b_0}}}. \end{equation}

The magnitude of the pressure gradient driving the input fluid, $1/{\eta _{b_0}}$, is greater than (or equal to) the magnitude of the pressure gradient driving the ambient fluid, ${\eta _{b_0}}/(2M)$, since ${\eta _{b_0}} \leqslant \sqrt {2M}$, with equality when $G \to 0^{+}$. This suggests that the trailing contact line is stable to small perturbations for $G>0$.

At the leading contact line (where $\mathcal {H}_0=0$), with $G \ll 1$, the pressure gradient in the ambient fluid (4.1b) is proportional to the relative viscosity of the ambient fluid divided by the radial distance. Similarly at $\mathcal {H}_0=1$, the pressure gradient in the input fluid (4.5b) is proportional to the relative viscosity of the input fluid divided by the radial distance. This dependence on the relative viscosity and radial location is as in the classical case (see (1.1)). However, the key difference is that in contrast to the classical instability, the contact lines are separated so that the radial coordinate is different at $\mathcal {H}_0=0$ and $\mathcal {H}_0=1$. The change in the pressure gradient between location (iv) and location (i) (figure 2a) is controlled by the viscosity contrast and how far apart the contact lines are, which reduces the magnitude of the pressure gradient downstream. These two effects act in opposition. If the intrusion of input fluid has a relatively short radial extent, then the flow is unstable and the intrusion grows radially until the distance between the contact lines acts to stabilise the interface. Hence the self-similar axisymmetric solutions have larger radial extent at smaller viscosity ratios. Stability is ensured because the radially extensive intrusion along the upper boundary associated with buoyancy segregation counteracts the destabilising discontinuity in the pressure gradient associated with the viscosity contrast. The vertically segregated radial intrusion could be thought of as a single axisymmetric viscous finger, which is stable to angular-dependent fingers.

It is not gravity that acts against the pressure gradient to stabilise the contact lines. The role of gravity is simply to segregate the fluids, and stability occurs after the segregation has occurred. Indeed, the radial extent of the interface is insensitive to $G$ at small values of $G$. In § 5, we show that provided that buoyancy has segregated the fluids, the small-$G$ axisymmetric self-similar solutions are linearly stable. This formally confirms the results in the present section.

5. Linear stability for small density difference ($G \ll 1$)

In this section, we demonstrate the linear stability of the $G \ll 1$ axisymmetric self-similar solutions (given by (2.34)). Note that $G>0$ as some buoyancy is required to segregate the fluids. Confirming the stability of the $G\ll 1$ solution will demonstrate stability for $G > 0$ as discussed in § 4. The linear stability analysis also gives the rate of decay of each mode and its dependence on the viscosity ratio.

We consider $\theta$-dependent perturbations to the axisymmetric self-similar solutions of the form

(5.1)\begin{gather} \mathcal{H}=\mathcal{H}_0(\eta) + \epsilon\,\mathcal{H}_1(\eta) \exp({\sigma \tau + {\rm i} n \theta}), \end{gather}

(5.2)\begin{gather}\eta_f ={\eta_{f_0}} + \epsilon {\eta_{f_1}} \exp({\sigma \tau + {\rm i} n \theta}), \end{gather}

(5.3)\begin{gather}\eta_b ={\eta_{b_0}} + \epsilon {\eta_{b_1}} \exp({\sigma \tau + {\rm i} n \theta}), \end{gather}

(5.4)\begin{gather}\mathcal{P} = \mathcal{P}_0(\eta) +\epsilon\,\mathcal{P}_1(\eta) \exp({\sigma \tau + {\rm i} n \theta}), \end{gather}

where $\epsilon \ll 1$ and $n$ is an integer. We seek to determine the stability of these perturbations. Note that the perturbation corresponds to ${\rm e}^{\sigma \tau } \sim T^{\sigma }$ in terms of real time $T$ (see (2.15a,b)).

Often, linear stability analyses of viscous gravity currents require rescaling the spatial domain owing to the singular behaviour of the interface at the contact lines (e.g. Mathunjwa & Hogg Reference Mathunjwa and Hogg2006; Kowal & Worster Reference Kowal and Worster2019). For the present porous gravity current, the interface is linear at the contact lines so such a transformation is not required. Instead, we linearise the boundary conditions about the leading-order locations of the contact lines, $\eta _f={\eta _{f_0}}$ and $\eta _b={\eta _{b_0}}$. We consider perturbations with $\theta$-dependence ($n \geqslant 1$) first as the stability in this case has not yet been investigated. For $n \geqslant 1$, global mass conservation of the input fluid (2.28) is identically satisfied by the form of the perturbations. The case of axisymmetric perturbations ($n=0$) is included at the end of this section for completeness.

In the single-phase regions (upstream of the trailing contact line and downstream of the leading contact line), the pressure $\mathcal {P}$ satisfies Laplace's equation, hence the pressure perturbation $\mathcal {P}_1$ is given by the solution of

(5.5)\begin{equation} \frac{1}{\eta}\,\frac{\mathrm{d}}{\mathrm{d} \eta}\left( \eta\, \frac{\mathrm{d} \mathcal{P}_1}{\mathrm{d} \eta} \right) - \frac{n^{2}}{\eta^{2}}\, \mathcal{P}_1 = 0. \end{equation}

Since the pressure perturbation remains finite as $\eta \to \infty$ and as $\eta \to 0$, the solution in the single-phase regions is

(5.6)\begin{gather} \mathcal{P}_1=c_n \eta^{n} \quad \text{for} \ \eta< {\eta_{b_0}}, \end{gather}

(5.7)\begin{gather}\mathcal{P}_1=d_n \eta^{{-}n} \quad \text{for} \ \eta > {\eta_{f_0}}, \end{gather}

where $c_n$ and $d_n$ are constants.

In the interface region ($\eta _{b_0}<\eta <\eta _{f_0}$), the linearised governing equations (2.17) and (2.18) become

(5.8)\begin{gather} (M-1)\,\frac{\mathrm{d}}{\mathrm{d} \eta} (\eta \mathcal{H}_1 ) + 2M \left( \frac{\mathrm{d}^{2} \mathcal{P}_1}{\mathrm{d} \eta^{2}} - \frac{n^{2}}{\eta^{2}}\, \mathcal{P}_1 \right) =0, \end{gather}

(5.9)\begin{gather}(1-M)(\sigma +1/2) \mathcal{H}_1 ={-}\frac{M}{\eta}\,\frac{\mathrm{d} \mathcal{P}_1}{\mathrm{d} \eta}. \end{gather}

We can then eliminate $\mathcal {H}_1$ and obtain a single equation for $\mathcal {P}_1$ in the interface region:

(5.10)\begin{equation} {\left(\frac{\sigma +1}{\sigma +1/2}\right) \frac{\mathrm{d}^{2} \mathcal{P}_1}{\mathrm{d} \eta^{2}} - \frac{n^{2}}{\eta^{2}}\,\mathcal{P}_1 = 0.} \end{equation}

To determine $\mathcal {P}_1$, we solve (5.10) with three boundary conditions at each contact line. These arise from the linearised version of continuity of the pressure, continuity of the flux, and the kinematic boundary condition (see (2.23a–c)–(2.27)). At $\eta ={\eta _{b_0}}^{+}$, these boundary conditions take the form

(5.11)\begin{gather} \mathcal{P}_1\rvert_{{\eta_{b_0}}^{+}} = c_n {\eta_{b_0}}^{n}, \end{gather}

(5.12)\begin{gather}\left.\frac{\mathrm{d} \mathcal{P}_1}{\mathrm{d} \eta} \right\vert_{{\eta_{b_0}}^{+}} = \frac{{\eta_{b_1}}}{{\eta_{b_0}}^{2}} + n c_n {\eta_{b_0}}^{n-1}, \end{gather}

(5.13)\begin{gather}(\sigma+1/2) {\eta_{b_1}} =\left. - M\,\frac{\mathrm{d} \mathcal{P}_1}{\mathrm{d} \eta} \right\vert_{{\eta_{b_0}}^{+}} ={-}M \left(\frac{{\eta_{b_1}}}{{\eta_{b_0}}^{2}} + n c_n {\eta_{b_0}}^{n-1}\right), \end{gather}

where we have used the upstream behaviour (5.6). At $\eta ={\eta _{f_0}}^{-}$, the analogous boundary conditions take the form

(5.14)\begin{gather} \mathcal{P}_1\rvert_{{\eta_{f_0}}^{-}} = d_n {\eta_{f_0}}^{{-}n}, \end{gather}

(5.15)\begin{gather}\left.\frac{\mathrm{d} \mathcal{P}_1}{\mathrm{d} \eta} \right\vert_{{\eta_{f_0}}^{-}}= \frac{{\eta_{f_1}}}{2} - n d_n {\eta_{f_0}}^{{-}n-1}, \end{gather}

(5.16)\begin{gather}(\sigma+1/2) {\eta_{f_1}} =\left. -\frac{\mathrm{d} \mathcal{P}_1}{\mathrm{d} \eta} \right\vert_{{\eta_{f_0}}^{-}}={-}\left(\frac{{\eta_{f_1}}}{2} - n d_n {\eta_{f_0}}^{{-}n-1}\right), \end{gather}

where we have used the downstream behaviour (5.7). The system for $\mathcal {P}_1$ between the contact lines is linear and has six boundary conditions with six unknown constants: $c_n$, $d_n$, ${\eta _{b_1}}$, ${\eta _{f_1}}$ and two constants arising from solving the ordinary differential equation (5.10). The system is an eigenvalue problem for the growth rate $\sigma$. The general solution of (5.10) is characterised by the value of $\sigma$ relative to the critical value

(5.17)\begin{equation} \sigma_c = \frac{-2n^{2} -1}{4n^{2} +1}. \end{equation}

This critical value of $\sigma$ corresponds to a repeated root in the auxiliary equation for the ordinary differential equation (5.10). There are no non-trivial solutions that satisfy all the boundary conditions for $\sigma \geqslant \sigma _c$. For $-1<\sigma <\sigma _c$, the solution takes the form

(5.18)\begin{equation} \mathcal{P}_1 = \eta^{1/2} \left[ a_1 \cos( \alpha \log \eta) + a_2 \sin( \alpha \log \eta)\right], \end{equation}

where

(5.19)\begin{equation} \alpha= \sqrt{-\frac{1}{4}-\left(\frac{\sigma+1/2}{\sigma +1} \right) n^{2}}, \end{equation}

and $a_1$ and $a_2$ are constants. We apply the six boundary conditions to obtain the following dispersion relation for $\sigma$:

(5.20)\begin{equation} {\mathcal{F}(\sigma,n,M) =0,} \end{equation}

where

(5.21)\begin{align} &\mathcal{F}(\sigma,n,M)\nonumber\\ &\quad = \left[\left(\sigma+{\tfrac{1}{2}}\right)n s_b -(\sigma+1)\left({\tfrac{1}{2}}s_b+\alpha c_b\right) \right] \left[ \left(\sigma+{\tfrac{1}{2}}\right)n c_f+(\sigma+1)\left({\tfrac{1}{2}}c_f-\alpha s_f\right) \right]\nonumber\\ &\qquad - \left[\left(\sigma+{\tfrac{1}{2}} \right)n s_f +(\sigma+1)\left({\tfrac{1}{2}}s_f+\alpha c_f\right)\right] \left[ \left(\sigma+{\tfrac{1}{2}}\right)n c_b-(\sigma+1)\left({\tfrac{1}{2}}c_b-\alpha s_b\right)\right] \end{align}

and

(5.22a,b)\begin{gather} s_b=\sin(\alpha \log {\eta_{b_0}}), \quad s_f=\sin(\alpha \log {\eta_{f_0}}), \end{gather}

(5.23a,b)\begin{gather}c_b=\cos(\alpha \log {\eta_{b_0}}), \quad c_f=\cos(\alpha \log {\eta_{f_0}}). \end{gather}

The locations of the contact lines $\eta _{b_0}$ and $\eta _{f_0}$ are functions of the viscosity ratio $M$, so the dispersion equation (5.20) also depends on $M$.

The dispersion equation (5.20) can be solved to obtain the growth rate $\sigma$ for any $n \geqslant 1$ and $0< M<1$. The corresponding eigenfunctions are given by (5.18). Since $\sigma < \sigma _c<0$, the growth rate is always negative, hence the axisymmetric similarity solutions are linearly stable.

The dispersion equation (5.20) can have multiple solutions. The function $\mathcal {F}(\sigma,n,M)$ is plotted against $\sigma$ in figure 6 for $M=0.5$ and four values of $n$. For each $n$, the largest zero for $\sigma$ is $\sigma =\sigma _c$ (these solutions are indicated by stars in figure 6). However, the case $\sigma =\sigma _c$ corresponds to a trivial solution for $\mathcal {P}_1$, which is not admissible. The next largest solution $\sigma$ of (5.20) corresponds to the slowest-decaying non-trivial solution, and these values are indicated by circles in figure 6. We expect this to be the decay rate that is observed, and we take this solution as the correct value of $\sigma$.

Figure 6. The dispersion relation function $\mathcal {F}(\sigma,n,M)$ (5.21) plotted versus the growth rate $\sigma$ for $M=0.5$. It has been normalised by its maximum value $\mathcal {F}_{{max}}$ over the range of $\sigma$. The zeros of the curves correspond to admissible decay rates $\sigma$ (5.20) with the exception of the largest zero (shown as stars). The next largest zero (shown as circles) will be the decay rate observed (further details in the text).

The eigenfunctions (5.18) that correspond to this correct solution $\sigma$ are defined up to a multiplicative constant. The form of the interfacial perturbation $\mathcal {H}_1$ for each eigenfunction (5.18) can be calculated from (5.9). These interfacial perturbation eigenfunctions are shown in figure 7 for different values of the viscosity ratio $M$ and the mode $n$. The shapes are more oscillatory at higher values of $n$.

Figure 7. Eigenfunctions for the interface perturbation $\mathcal {H}_1$, where we have set $\mathcal {H}_1(\eta =1)=1$: (a) $M=0.2$, and (b) $M=0.5$. The eigenfunctions are calculated from (5.9) and (5.18). The horizontal axis is $\eta _{b_0}< \eta < \eta _{f_0}$ for (a) $\eta _{b_0}=0.63$, $\eta _{f_0}=3.16$, and (b) $\eta _{b_0}=1$, $\eta _{f_0}=2$.

The decay rate $\sigma$ calculated for different values of $M$ and $n$ using the dispersion relation (5.20) is plotted using continuous lines in figure 8. The rate of decay is slower for larger values of $n$ and smaller viscosity ratios $M$. The crosses in figure 8 denote the decay rate calculated from the numerical method for $M=0.5$ and $M=0.8$ for three values of $n$ with $G=0.01$; details of this calculation are given in § A.1. There is excellent agreement between the numerically derived prediction for $\sigma$ and the values obtained from the linear stability analysis.

Figure 8. The decay rate $\sigma$ for $M=0.1$, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9 and $0 \leqslant n \leqslant 10$. The crosses show the decay rate observed from our numerical results (see § A.1).

Finally, we consider the case $n=0$. For $\eta <{\eta _{b_0}}$, the pressure perturbation is constant, $\mathcal {P}_1=c_0$, whilst for $\eta >{\eta _{f_0}}$, $\mathcal {P}_1=d_0$. In the interface region, the pressure perturbation is linear: $\mathcal {P}_1=a+b\eta$. The kinematic boundary conditions (5.13) and (5.16) furnish $\sigma =-1$. The thickness perturbation takes the form (see (5.9))

(5.24)\begin{equation} \mathcal{H}_1 = \frac{2 M b}{1-M}\,\eta^{{-}1}. \end{equation}

The numerical results for axisymmetric perturbations to the self-similar solutions corroborated that $\sigma =-1$ (shown as a cross at $n=0$ in figure 8). Axisymmetric perturbations decay faster than $\theta$-dependent perturbations, with the error decaying proportional to $T^{-1}$ as has been found previously for unconfined gravity currents (Grundy & McLaughlin Reference Grundy and McLaughlin1982; Mathunjwa & Hogg Reference Mathunjwa and Hogg2006). For a different stability argument in the case of axisymmetric perturbations, see the Appendix of Nordbotten & Celia (Reference Nordbotten and Celia2006).

6. Vertically varying permeability

The results obtained in previous sections generalise to layers in which the permeability varies vertically (i.e. as a function of $Z$). The flow is axisymmetric and self-similar after buoyancy segregation. Here, we find the analytical solutions that arise in the small-$G$ regime with a less viscous input fluid ($M<1$). One key difference from the uniform case is that shock-like regions in which the interface is relatively steep can occur even for $M<1$ (see § 6.1). Despite this, we show that at late times, the axisymmetric solutions for any continuous permeability profile are stable.

We define the dimensionless horizontal permeability as

(6.1)\begin{equation} \mathcal{K}_x(Z/H_0) = \frac{k_x(Z)}{\dfrac{1}{H_0} \displaystyle\int_0^{H_0} k_x(Z)\,\mathrm{d}Z}, \end{equation}

where the denominator is the mean horizontal permeability. We also define the dimensionless depth-integrated horizontal permeability (Hinton & Woods Reference Hinton and Woods2018)

(6.2)\begin{equation} \varPhi(\mathcal{H}) = \int_0^{\mathcal{H}} \mathcal{K}_x(s) \, \mathrm{d} s. \end{equation}

Note that $\varPhi (0)=0$ and $\varPhi (1)=1$. For a uniform layer, $\mathcal {K}_x \equiv 1$ and $\varPhi (\mathcal {H})=\mathcal {H}$. The two mass conservation equations (2.17) and (2.18) generalise to

(6.3)\begin{gather} \frac{\partial \mathcal{H}}{\partial \tau}-\frac{\eta}{2}\,\frac{\partial \mathcal{H}}{\partial \eta} - \tilde{\boldsymbol{\nabla}} \boldsymbol{\cdot} \left( \varPhi(\mathcal{H})\,\tilde{\boldsymbol{\nabla}} \mathcal{P} \right) =0, \end{gather}

(6.4)\begin{gather}\tilde{\boldsymbol{\nabla}} \boldsymbol{\cdot} \left[\left(M+(1-M)\,\varPhi(\mathcal{H})\right) \tilde{\boldsymbol{\nabla}} \mathcal{P} - G M (1-\varPhi(\mathcal{H}))\,\tilde{\boldsymbol{\nabla}} \mathcal{H} \right] = 0, \end{gather}

where the horizontal permeability in the definition of $G$ (2.20a,b) is now given by the mean horizontal permeability. This model applies provided that the fluids have become vertically segregated by buoyancy. In a uniform layer, this applies at times given by (2.41). In a layer with vertically varying permeability and less viscous input fluid ($M<1$), the dimensional time for buoyancy segregation is adjusted to

(6.5)\begin{equation} {T \gg \frac{\mu_a}{\Delta \rho\,g} \int_0^{H_0} \frac{1}{k_z(Z)} \, \mathrm{d} Z,} \end{equation}

where $k_z(Z)$ is the vertical permeability. In an isotropic vertically heterogeneous layer, $k_z(Z)=k_x(Z)$. The time (6.5) becomes singular if the vertical permeability is zero at any height in the layer (and in this case our model does not apply). The time for the shallow approximation to apply is given by (2.42).

As for a uniform layer, there are self-similar axisymmetric solutions $\mathcal {H}_0=\mathcal {H}_0(\eta )$, $\mathcal {P}=\mathcal {P}_0(\eta )$ to (6.3) and (6.4), with the pressure gradient at the top boundary given by

(6.6)\begin{equation} \frac{\mathrm{d} \mathcal{P}_0}{\mathrm{d} \eta} = \frac{-1}{\eta\left(M+(1-M)\,\varPhi(\mathcal{H}_0)\right)} + \frac{GM (1-\varPhi(\mathcal{H}_0))}{M+(1-M)\,\varPhi(\mathcal{H}_0)}\,\frac{\mathrm{d} \mathcal{H}_0}{\mathrm{d} \eta}, \end{equation}

and $\mathcal {H}_0$ satisfies the ordinary differential equation

(6.7)\begin{align} &{-\frac{\eta}{2}}\,\frac{\mathrm{d} \mathcal{H}_0}{\mathrm{d} \eta} + \frac{1}{\eta}\,\frac{\mathrm{d}}{\mathrm{d} \eta} \left(\frac{\varPhi(\mathcal{H}_0)}{M+(1-M)\,\varPhi(\mathcal{H}_0)} \right) \nonumber\\ &\qquad=\frac{1}{\eta}\,\frac{\mathrm{d}}{\mathrm{d} \eta} \left(\eta\,\frac{G M\, \varPhi(\mathcal{H}_0)\,(1-\varPhi(\mathcal{H}_0))}{M+(1-M)\,\varPhi(\mathcal{H}_0)}\,\frac{\mathrm{d} \mathcal{H}_0}{\mathrm{d} \eta} \right). \end{align}

The small-$G$ analytic solutions are given by neglecting the right-hand side of (6.7), from which we obtain the implicit solution

(6.8)\begin{gather} \mathcal{H}_0={1,} \quad 0<\eta< {\eta_{b_0}}, \end{gather}

(6.9)\begin{gather}\eta=\frac{\sqrt{2M\,\mathcal{K}_x(\mathcal{H}_0)}}{M+(1-M)\,\varPhi(\mathcal{H}_0)}, \quad 0<\mathcal{H}_0< 1, \end{gather}

(6.10)\begin{gather}\mathcal{H}_0={0,} \quad {\eta_{f_0}}<\eta, \end{gather}

where the contact lines are

(6.11a,b)\begin{equation} {\eta_{f_0}} = \sqrt{\frac{2\,\mathcal{K}_x(0)}{M}}, \quad {\eta_{b_0}} =\sqrt{2 M\,\mathcal{K}_x(1)}. \end{equation}

The associated pressure gradient in the interface region is given by

(6.12)\begin{equation} \frac{\mathrm{d} \mathcal{P}_0}{\mathrm{d} \eta} ={-}\frac{1}{\sqrt{2M\, \mathcal{K}_x(\mathcal{H}_0)}}. \end{equation}

In a uniform layer, the interface (6.9) is monotonic whenever $M<1$, but in a layer with a vertical permeability variation, turning points can arise even when the input fluid is less viscous than the ambient (Hinton & Woods Reference Hinton and Woods2018). The solution (6.9) is valid only when the interface is monotonic. Otherwise, the interface has a turning point with buoyant input fluid lying below denser ambient fluid. This would invalidate the model, which assumed that the fluids have segregated owing to buoyancy. For example, for a linear permeability structure with dimensionless variation $\Delta k$,

(6.13)\begin{equation} \mathcal{K}_x(s) = 1+ \Delta k\,(s-1/2), \end{equation}

the interface (6.9) is monotonic if and only if (Hinton & Woods Reference Hinton and Woods2018)

(6.14)\begin{equation} M < \frac{(2-\Delta k)^{2}}{4 - 2\,\Delta k + (\Delta k)^{2}}. \end{equation}

Provided that (6.9) is monotonic, the interface forms an axisymmetric self-similar intrusion, which extends along the upper boundary in a qualitatively identical fashion to the case of a uniform layer (see figure 9a). In the case that (6.9) has a turning point, a shock must be introduced; this regime is studied in § 6.1.

Figure 9. Self-similar axisymmetric interface shapes in a layer with vertically varying permeability. (a) Small-$G$ analytic solutions (6.9) for different linear permeability variations (6.13). (b) Small-$G$ analytic solutions in the case that there is a shock-like region.

We apply our numerical method to layers with vertical variations in permeability for a wide range of initial conditions, and find that the self-similar axisymmetric solutions are stable to both $\theta$-dependent and axisymmetric perturbations for any $G > 0$. Next, we generalise § 4 to show how the pressure gradients suppress instabilities in a layer with vertically varying permeability (but no shock-like regions; for that case, see § 6.1). Between the contact lines, the magnitude of the pressure gradient in the ambient fluid is less than or equal to that in the input fluid owing to the contribution of the term $-G\,\partial \mathcal {H}_0/\partial \eta$ (as in § 4). The interface gradient at the contact lines is discontinuous, and the stability there is treated separately.

First, we obtain bounds on the location of the contact lines of the self-similar axisymmetric solution. In the case that the small-$G$ analytic solution (6.9) is monotonic, the contact lines for $G > 0$ satisfy

(6.15a,b)\begin{equation} {\eta_{f_0}} \geqslant \sqrt{\frac{2\,\mathcal{K}_x(0)}{M}}, \quad {\eta_{b_0}} \leqslant \sqrt{2 M\,\mathcal{K}_x(1)}, \end{equation}

with equality as $G\to 0^{+}$. These inequalities reflect the action of buoyancy to extend the interface.

The interface gradient and pressure gradient in the ambient fluid just ahead of the leading contact line (at $\eta =\eta _{f_0}^{+}$) are given by

(6.16a,b)\begin{equation} \frac{\mathrm{d} \mathcal{H}_0}{\partial \eta} =0, \quad \frac{\mathrm{d}(\mathcal{P}_0-G \mathcal{H}_0)}{\mathrm{d} \eta} ={-} \frac{1}{M {\eta_{f_0}}}. \end{equation}

The interface gradient and pressure gradient just behind the leading contact line (at $\eta =\eta _{f_0}^{-}$) are

(6.17a,b)\begin{equation} \frac{\mathrm{d} \mathcal{H}_0}{\partial \eta} = \frac{1}{{\eta_{f_0}} M G} - \frac{\eta_{f_0}}{2 G\,\mathcal{K}_x(0)}, \quad \frac{\mathrm{d} \mathcal{P}_0}{\partial \eta} = \frac{-{\eta_{f_0}}}{2\,\mathcal{K}_x(0)}, \end{equation}

respectively. The former is obtained by taking $\mathcal {H}_0 \to 0^{+}$ in (6.7). The inequality (6.15a) ensures that the magnitude of the pressure gradient is larger in the input fluid, which stabilises the interface, with equality as $G\to 0^{+}$.

Just ahead of the trailing contact line (at $\eta =\eta _{b_0}^{+}$),

(6.18a,b)\begin{equation} \frac{\mathrm{d} \mathcal{H}_0}{\mathrm{d} \eta} = \frac{{\eta_{b_0}}}{2 M G\, \mathcal{K}_x(1)} - \frac{1}{{\eta_{b_0}} G}, \quad {\frac{\mathrm{d} \mathcal{P}_0}{\mathrm{d} \eta} = \frac{-1}{{\eta_{b_0}}}.} \end{equation}

Thus the pressure gradient in the ambient fluid just ahead of the trailing contact line is given by

(6.19)\begin{equation} \frac{\mathrm{d}(\mathcal{P}_0-G \mathcal{H}_0)}{\mathrm{d} \eta} ={-} \frac{{\eta_{b_0}}}{2 M\,\mathcal{K}_x(1)}, \end{equation}

whilst in the input fluid just upstream of the trailing contact line (at $\eta =\eta _{b_0}^{-}$), the pressure gradient is

(6.20)\begin{equation} \frac{\mathrm{d} \mathcal{P}_0}{\mathrm{d} \eta} ={-} \frac{1}{{\eta_{b_0}}}, \end{equation}

and the inequality (6.15b) ensures that the magnitude of the pressure gradient is larger in the input fluid, which stabilises the interface.

6.1. Shock-like regions

In the case in which the small-$G$ interface (6.9) has a turning point, a shock must be introduced to vertically segregate the fluids. This shock may occupy part or all of the layer (see figure 9b). By way of an example, we consider a linear variation in permeability with $\Delta k=1.5$ and viscosity ratio $M=0.3$. The axisymmetric self-similar profile has a shock-like region near the top of the layer, even though $M<1$ (red line in figure 9b). The shock at the top of the layer is at location $\eta =\eta _s$ and is of thickness $\mathcal {H}_s$, which satisfy

(6.21)\begin{equation} \eta_s=\sqrt{2\,F'(\mathcal{H}_s)}, \quad \eta_s^{2} \mathcal{H}_s = 2\,F(\mathcal{H}_s), \quad \text{where} \ F(\mathcal{H}) = \frac{\varPhi(\mathcal{H})}{M+(1-M)\,\varPhi(\mathcal{H})}. \end{equation}

These arise from continuity with the solution (6.9), which is valid for $\mathcal {H}>\mathcal {H}_s$, and mass conservation, respectively. The shock-like region may also extend across the entire layer (e.g. black line in figure 9b). For more complicated permeability profiles, there can be multiple shock-like regions separated by rarefaction-like regions (Hinton & Woods Reference Hinton and Woods2018).

We apply our numerical method to layers with shock-like regions and find that the self-similar axisymmetric solutions are stable to both $\theta$-dependent and axisymmetric perturbations. An example is shown in figure 10 with a linear variation in permeability with $\Delta k=1.5$, viscosity ratio $M=0.3$, and $G=0.05$, demonstrating that mode five fingers are suppressed. The small-$G$ analytic solution is shown as a dashed line in figure 10(c); there is a shock-like region in the top part of the layer. For $G>0$, the shock-like regions are smoothed by buoyancy so that they are not vertical (Hinton & Woods Reference Hinton and Woods2018). The smoothed shock-like region has radial extent proportional to $G$ in the similarity coordinate $\eta$. Hence, in dimensional terms, the aspect ratio of the shock evolves with

(6.22)\begin{equation} {\frac{L}{H_0} \sim G \left(\frac{QT}{2 {\rm \pi}\phi H_0^{3}} \right)^{1/2},} \end{equation}

and eventually the shock-like region of the flow is long and thin (when the right-hand side becomes large) and the shallow flow approximation applies at sufficiently late times. This smoothing of the shock-like region also occurs when the shock encompasses the entire thickness of the layer (e.g. figure 9b). Thus the fluid–fluid interface is not cylindrical in this case. Instead, due to buoyancy, the interface becomes gradually shallower over time, which suppresses the classical instability.

Figure 10. Shutdown of mode five fingers in a layer with vertically varying permeability. We use a linear variation with $\Delta k=1.5$, viscosity ratio $M=0.3$, and $G=0.05$. Flow thickness is shown at three radial cross-sections at times (a) $\tau =0.5$, (b) $\tau =1.5$, and (c) $\tau =20$ (results calculated numerically; see § 3). In (c), the black dashed line shows the axisymmetric small-$G$ solution with a shock-like region in the upper portion of the layer. (d) Location of the trailing contact line $\eta =\eta _b(\theta,\tau )$ at unit time intervals between $\tau =0$ and $\tau =20$ (lighter colours correspond to later times).

The stability analysis generalises to the case in which there are shock-like regions. At the top of the layer, the leading contact line for $G > 0$ satisfies

(6.23)\begin{equation} {\eta_{f_0}}\geqslant \eta_s \geqslant \sqrt{\frac{2\,\mathcal{K}_x(0)}{M}}. \end{equation}

The first inequality arises from buoyancy smoothing the interface beyond the shock location (see figure 10c), whilst the second arises because the profile (6.9) has a turning point in the shock-like region, and the shock conserves mass so it must be ahead of this profile at $Z=0$. Hence (6.15a) applies even when there is a shock-like region at the top of the layer. A similar argument applies when there is a shock-like region at the bottom of the layer and (6.15b) remains correct.

The gradient of the interface and the pressure at the upper boundary at the contact lines are given by (6.18a,b) and (6.19a,b), even when there are shock-like regions. Hence the stability argument of § 6 applies.

In summary, the self-similar axisymmetric solutions in a layer with vertically varying permeability are stable owing to the same mechanism as in a uniform layer. The less viscous input fluid intrudes into the ambient fluid sufficiently far that the destabilising increase in the magnitude of the pressure gradient associated with the viscosity contrast is smoothed over the extent of the interface. Buoyancy segregation corresponds to a monotonic interface, and the input fluid may be forced into the low-permeability region (see also Debbabi et al. Reference Debbabi, Jackson, Hampson and Salinas2018). We note that there may be significant intermingling of the fluids in the transient evolution to the buoyancy-segregated flow (e.g. Huppert, Neufeld & Strandkvist Reference Huppert, Neufeld and Strandkvist2013). There may also be a competition between viscous fingering and heterogeneity-driven fingering prior to buoyancy segregation (De Wit & Homsy Reference De Wit and Homsy1997).