2.1. A note on anisotropic permeability

Before discussing the flow scenario, we first briefly discuss the geological context of anisotropy in porous media. In subsurface geological reservoirs it is common for the permeability of rocks to differ significantly depending on the direction of the flow (Corbett & Jensen Reference Corbett and Jensen1992; Martinius et al. Reference Martinius, Ringrose, Næss, Wen and Hentz1999). This may result from post-depositional compaction of the formation or from the deposition of successive layers of fine and coarse material. Hence, the permeability field is sometimes written as a three-dimensional diagonal matrix, $\boldsymbol {k}=\mathrm {Diag}(k_x,k_y,k_z)$, or equivalently in a different basis depending on the direction of compaction (e.g. see studies on cross-bedding Woods Reference Woods2015).

We account for such situations, but we restrict our attention to the case where the compaction is aligned with the vertical coordinate, and hence the permeability is reduced to horizontal and vertical variation only, $(k_x,k_y,k_z)=(k_H,k_H,k_V)$, where $k_H$ and $k_V$ are constants. Hence, the anisotropy is characterised by the single dimensionless parameter

(2.1)\begin{equation} \gamma=k_H/k_V. \end{equation}

This formulation is relevant either for flow through a single thick sedimentary layer that has been compacted vertically, resulting in a binary permeability field, or for flow through a system of many horizontal sedimentary layers, in which the values $k_H$, $k_V$, can be interpreted as effective permeabilities (Cardwell & Parsons Reference Cardwell and Parsons1945; Kumar et al. Reference Kumar, Farmer, Jerauld and Li1997; Woods Reference Woods2015). For example, in the latter scenario, if the layers vary between permeabilities $k_1$ and $k_2$ over alternating layer widths $d_1$ and $d_2$, then the effective permeability values are given by the arithmetic and harmonic mean values

(2.2a,b)\begin{equation} k_H\approx \frac{d_1 k_1 +d_2k_2}{d_1+d_2},\quad k_V\approx \frac{d_1+d_2}{d_1/k_1+d_2/ k_2}. \end{equation}

This tends to be a good approximation as long as the depth of the flow is much larger than the widths of the layers $d_1,d_2$. Hence, in this system of horizontal sedimentary layers the anisotropy is always such that $\gamma \geq 1$. In fact, for many geological systems the horizontal and vertical permeability have been observed to differ by several orders of magnitude (Martinius et al. Reference Martinius, Ringrose, Næss, Wen and Hentz1999), resulting in anisotropy values as large as $\gamma ={O}(10^4)$.

2.2. The flow scenario at early times and subsequent transition behaviour

Next, we move on to describe the flow scenario we consider, and describe how this evolves at early times, before the effects of gravity become significant. For the isotropic case, as discussed by Huppert & Pegler (Reference Huppert and Pegler2022), the constant input of fluid from a point source results in the current expanding in the shape of a self-similar hemisphere at early times. The radial and vertical extent of the current are equal and grow like the cube root of time. The flow transitions to a gravity-driven regime once the weight of the fluid dominates over the injection pressure. In the following analysis, we extend this early-time analysis to the case of injection into an anisotropic medium, demonstrating how the anisotropy $\gamma$ affects both the early-time dynamics and the transition behaviour.

We consider the constant axisymmetric injection of fluid (at flow rate $Q$) into a porous medium with anisotropy $\gamma$ and porosity $\phi$, bounded from below by an impermeable substrate, as illustrated in figure 1. The porous medium is initially saturated with an ambient phase which has a lighter density than the injected fluid $\rho _2<\rho _1$, whereas the viscosities of the two fluids are assumed to be equal $\mu _1=\mu _2=\mu$. For simplicity, we assume that there is no mixing between the two fluids, such that the interface between them remains sharp. We choose a cylindrical polar coordinate system $(r\geq 0,\ z\geq 0)$ and we denote the radial and vertical extent of the current, $R(t)$ and $H(t)$, which increase over time due to the constant injection rate $Q$ from a point at the origin.

Figure 1. Schematic diagram illustrating the transition from a pressure-driven flow at early times (a) to a gravity-driven flow with a pressure-driven boundary layer near the origin at late times (b) in an anisotropic porous medium. Contour lines are indicative of the typical pressure field and approximate streamlines are sketched: (a) $t\ll t^*$; (b) $t\gg t^*$.

At early times the pressure is dominated by the viscous resistance due to injection and gravity has a negligible effect. Hence, due to conservation of mass the pressure must satisfy the anisotropic version of Laplace's equation (i.e. with non-uniform coefficients), such that

(2.3)\begin{equation} k_H\frac{1}{r}\frac{\partial}{\partial r} \left(r\frac{\partial p}{\partial r}\right) +k_V\frac{\partial^2 p}{\partial z^2}=0. \end{equation}

Since $k_H$ and $k_V$ are constants, the dynamics can be easily mapped from an isotropic medium, and hence the lateral and vertical extents must satisfy

(2.4)\begin{equation} H^2/R^2=1/\gamma.\end{equation}

Likewise, conservation of mass (within a stretched hemisphere) dictates that

(2.5)\begin{equation} 2{\rm \pi} \phi HR^2/3 = Q t.\end{equation}

This indicates that the early-time dynamics is described by

(2.6)$$\begin{gather} H=\gamma^{{-}1/3}(3Qt/2{\rm \pi}\phi)^{1/3}, \end{gather}$$

(2.7)$$\begin{gather}R=\gamma^{1/6}(3Qt/2{\rm \pi}\phi)^{1/3}. \end{gather}$$

Hence, at early times, the effect of anisotropy is to stretch/squash the flow in the more/less permeable direction. We note that it is also possible to extend (2.6)–(2.7) to the case of a time-varying injection rate $Q(t)$, although we do not include this analysis here.

From this early-time analysis it follows that the thickness of the current, $z=h(r,t)$, is given by

(2.8)\begin{equation} h/H=[1-(r/R)^2]^{1/2},\end{equation}

at early times. This is plotted in figure 2(a,c) for both isotropic ($\gamma =1$) and anisotropic ($\gamma =10$) porous media. Numerical data are also plotted (at both early (a,c) and late times (b,d)), which are discussed in more detail in § 2.5. It is surprising to note that (within this early-time regime) a strongly anisotropic medium results in a current which is very long and thin, but in which the pressure is not hydrostatic. This is contrary to the common assumption that long and thin currents automatically imply a hydrostatic pressure.

Figure 2. Evolution of the shape of the injected flow $z=h(r,t)$ in isotropic (a,b) and anisotropic (c,d) porous media. Analytical (approximate) solutions (2.8), (2.29) are shown at early times in red (a,c) and at late times in blue (b,d). Numerical solutions (black dotted lines) are shown at logarithmically spaced time intervals with $t/t^*=6\times (10^{-7}\unicode{x2013}10^{3})$ (other data for larger times do not fit within the axis limits). The position of the inner radius $r=R_p(t)$ is shown at late times with red points. The vertical scale of (b,d) is stretched by $\times 10.5$ for illustration purposes (note, (a,c) are magnified regions of (b,d)).

After a significant amount of time the current grows to a thickness where the weight of fluid begins to dominate over the injection pressures. This is equivalent to the moment when the input flow rate $Q$ is matched by the flow rate due to hydrostatic pressure gradients. For a hydrostatic current, we have that $\partial p/\partial r\approx \Delta \rho g \partial h/\partial r$, so that the integrated gravity-driven flux scales like

(2.9)\begin{equation} Q_g\sim{-}\frac{k_H\Delta \rho g}{\mu} r h\frac{\partial h}{\partial r}\sim Q.\end{equation}

By inserting approximate scalings, $h\sim H$, $r\sim R$, $\partial h/\partial r\sim -H/R$, this provides the balance

(2.10)\begin{equation} Q_g\sim u_b H^2 \sim Q,\end{equation}

where the buoyancy velocity $u_b=k_H\Delta \rho g/\mu$. This immediately provides an expression for the thickness of the current at the transition between pressure- to buoyancy-driven flow,

(2.11)\begin{equation} H^*=\left( \frac{Q}{u_b}\right)^{1/2},\end{equation}

and by equating this to the early-time behaviour (2.6) we find that the transition time is

(2.12)\begin{equation} t^*= \frac{2\gamma\phi{\rm \pi}}{3}\left(\frac{Q}{u_b^3}\right)^{1/2}.\end{equation}

Therefore, anisotropy causes the transition to occur at later times, as expected. Indeed, the significance of this delayed transition is immediately appreciable when one considers that some geological formations have anisotropy values as large as $\gamma ={O}(10^4)$ (Martinius et al. Reference Martinius, Ringrose, Næss, Wen and Hentz1999; Bergmo et al. Reference Bergmo, Emmel, Anthonsen, Aagaard, Mortensen and Sundal2017).

The critical time $t^*$ distinguishes two distinct regimes. At early times $t\ll t^*$, or when $H\ll H^*$, the effects of gravity are negligible and injection pressures dominate the flow, whereas at late times $t\gg t^*$, or when $H\gg H^*$, gravity plays a significant role. However, clearly the flow must remain pressure-driven very close to the source, even at late times. This pressure-driven region is responsible for distributing the flow across the vertical extent of the current, whereupon it diverts to the remaining gravity-driven region far away from the origin (see figure 1b). Whilst various authors have alluded to the existence of this pressure-driven boundary layer, it is not known how its lateral and vertical extents evolve over time, nor how it couples with the remaining gravity current. Furthermore, it is unknown how anisotropy affects the flow after transition occurs, which is the subject of the following sections.

2.3. Late-time dynamics: classical analysis

At much later times $t\gg t^*$, once the flow has grown to a thickness $H\gg H^*$, it retains much of the character of a classical porous gravity current. There is a long history of studies on self-similar and gravity-driven flows in porous media, with the earliest work developed by Pattle (Reference Pattle1959) and Barenblatt (Reference Barenblatt1952). Later, Lyle et al. (Reference Lyle, Huppert, Hallworth, Bickle and Chadwick2005) described the late-time dynamics in an isotropic porous medium by assuming a hydrostatic pressure profile everywhere within the current. This results in the corresponding radially symmetric thin-film equation

(2.13)\begin{equation} \phi\frac{\partial h}{\partial t}=u_b\frac{1}{r}\frac{\partial }{\partial r} \left[r h \frac{\partial h}{\partial r}\right],\end{equation}

for the evolution of the current thickness $z=h(r,t)$. The above equation is accompanied by boundary conditions imposing zero thickness at the nose, constant injection flux at the origin and zero flux through the nose, such that

(2.14)$$\begin{gather} h=0,\quad \mathrm{at}\ r=R(t), \end{gather}$$

(2.15)$$\begin{gather}-2{\rm \pi} r h \frac{\partial h}{\partial r}\rightarrow Q,\quad \mathrm{as}\ r\rightarrow 0, \end{gather}$$

(2.16)$$\begin{gather}-2{\rm \pi} r h \frac{\partial h}{\partial r}\rightarrow 0,\quad \mathrm{as} r\rightarrow R(t). \end{gather}$$

It should be noted that, whilst these equations were formulated for the case of an isotropic porous medium, the analysis can be easily transposed to the case of anisotropic permeability, $k_H$, $k_V$. However, the vertical permeability $k_V$ vanishes from the analysis since the injected flow and the ambient flow are decoupled (an assumption which we revisit later), resulting in an identical system of equations.

The system of (2.13)–(2.16) admits the similarity variables

(2.17a,b)\begin{equation} h=\left(\frac{Q}{u_b}\right)^{1/2} f(\eta/\eta_N),\quad \eta =\frac{r}{(Q u_b)^{1/4}(t/\phi)^{1/2}},\end{equation}

where the dimensionless shape $f$ and prefactor $\eta _N$ satisfy

(2.18)$$\begin{gather} -\frac{1}{2}\eta f'=\frac{1}{\eta}\frac{\mathrm{d}}{\mathrm{d}\eta} \left[\eta f \frac{\mathrm{d} f}{\mathrm{d}\eta}\right], \end{gather}$$

(2.19)$$\begin{gather}f=0, \quad \mathrm{at}\ \eta = \eta_N, \end{gather}$$

(2.20)$$\begin{gather}-2{\rm \pi} \eta f \frac{\mathrm{d} f}{\mathrm{d}\eta} \rightarrow 1,\quad \mathrm{as}\ \eta \rightarrow 0, \end{gather}$$

(2.21)$$\begin{gather}-2{\rm \pi} \eta f \frac{\mathrm{d} f}{\mathrm{d}\eta} \rightarrow 0,\quad \mathrm{as}\ \eta \rightarrow \eta_N. \end{gather}$$

The solution can be calculated numerically, giving a prefactor value $\eta _N=1.155$, and the resultant shape is plotted in figure 8 in Appendix B with green lines. The influx condition (2.20) results in singular behaviour of the thickness $h$ near the origin, which is a consequence of the hydrostatic pressure assumption. This is unphysical and can be addressed by considering the pressure-driven flow near the origin. In particular, the flux within this region is not buoyancy driven, and so the inflow boundary condition (2.20) (resulting in singular behaviour of the thickness) no longer applies. Instead, the inflow condition requires a non-hydrostatic flux component within the pressure-driven zone, and as we will see, this maintains a dominant contribution to the flow (even at late times) such that the thin-film approximation (2.13) cannot be used near the origin.

2.4. Back to the future: late-time dynamics revisited

Within the inner pressure-driven region strong vertical velocities deliver the injected flow all the way to the uppermost extent of the current. The inner pressure-driven radius and the thickness scaling are related (for the same reasons as (2.4)) according to

(2.22)\begin{equation} H^2/R_p^2=1/\gamma.\end{equation}

Our isotropic porous bead experiments, discussed later in § 2.6, show that the thickness of the current near the origin increases slowly over time (e.g. see figure 6). Therefore, we attempt to describe the solution within this inner region using a set of rescaled variables where we allow for the possibility of an evolving thickness at the origin, as shown in figure 3. Using a weak power law value $0< a<1$, the variables are rescaled according to

(2.23a–c)\begin{align} \xi=r/R_p(t),\quad h=\gamma^{{-}1/2}\mathcal{G}(\xi,t)R_p(t),\quad R_p(t)=\xi_0 u_b^{{(3a-1)}/{2}}Q^{{(1-a)}/{2}} (t/\phi)^a,\end{align}

where $\mathcal {G}$ is the dimensionless shape function, $\xi _0$ is a dimensionless constant that we determine later, and the powers of $u_b$ and $Q$ in (2.23a–c) are chosen to be dimensionally consistent. Note that we have kept the time dependence in $\mathcal {G}$ since it is not clear that the shape is self-similar under the proposed scalings with a growing thickness. We later show that this dependence can be neglected to good approximation at very late times.

Figure 3. Schematic diagram showing the different variables of the late-time regime. The analysis in § 2.4 shows that the power law values are $a=1/7,\ b=3/7$. Note the aspect ratio is exaggerated for illustration purposes, but in fact $H$ may be smaller than $R_p$ due to anisotropy (2.22).

The consequence of having a growing inner region is that the radial extent of the corresponding outer region does not grow according to the classical $R\sim t^{1/2}$ power law described earlier. In particular, conservation of mass (in conjunction with (2.22)) dictates that

(2.24)\begin{equation} {\rm \pi}\phi H R^2= {\rm \pi}\phi \gamma^{{-}1/2}R_p R^2\approx Q t.\end{equation}

If $a>0$, it follows that $R$ grows as ${\sim }t^b$, where $b<1/2$. Hence, we introduce a corresponding set of rescaled outer variables

(2.25a–c)\begin{equation} \zeta=r/R(t),\quad h=\gamma^{{-}1/2}\mathcal{F}(\zeta,t)R_p(t),\quad R(t)=\zeta_N u_b^{{(3b-1)}/{2}}Q^{{(1-b)}/{2}} (t/\phi)^b,\end{equation}

where $\mathcal {F}$ is the dimensionless shape profile, the nose position $\zeta _N$ is a dimensionless constant to be determined and $h$ is rescaled to match with the inner region. Mass conservation (2.24) indicates that

(2.26)\begin{equation} a+2b=1.\end{equation}

Given these rescaled variables (2.25a–c), the thin-film equation (2.13) in the outer (hydrostatic) region becomes

(2.27)\begin{equation} \frac{R^2}{R_p}\frac{\partial \mathcal{F}}{\partial t}-\frac{R R'}{R_p}\zeta \frac{\partial \mathcal{F}}{\partial \zeta}+\frac{R^2 R_p'}{R_p^2}\mathcal{F}= \frac{u_b}{ \gamma^{1/2}\phi\zeta}\frac{\partial }{\partial \zeta}\left[\zeta \mathcal{F} \frac{\partial \mathcal{F}}{\partial \zeta}\right].\end{equation}

The second two terms on the left-hand side of (2.27) are of the order ${O}(t^{-2a})$ and so become vanishingly small at large times. Assuming the time derivative of the scaled shape decreases as ${\partial \mathcal {F}}/{\partial t}\sim 1/t$ for large times, then the first term on the left-hand side is also of the order ${O}(t^{-2a})$. Hence, at sufficiently late times (2.27) yields the result

(2.28)\begin{equation} {-}2{\rm \pi} \zeta \mathcal{F} \frac{\mathrm{d} \mathcal{F}}{\mathrm{d} \zeta}\approx \mathcal{Q},\end{equation}

where $\mathcal {Q}$ is the scaled flux of injected fluid, which is uniform in magnitude across the extent of the gravity current, and which is yet unknown. Upon further integration and applying the boundary condition (2.14), we arrive at an expression for the outer solution, which is

(2.29)\begin{equation} \mathcal{F}\approx\left(-\frac{\mathcal{Q}}{\rm \pi} \log \zeta\right)^{1/2}.\end{equation}

We note that the integrated gravity-driven flux is given by

(2.30)\begin{equation} Q_g={-}2{\rm \pi} u_b {r} {h} \frac{\partial {h}}{\partial {r}},\end{equation}

which grows according to

(2.31)\begin{equation} Q_g=\gamma^{{-}1}\mathcal{Q}u_b R_p(t)^2.\end{equation}

At the transition thickness, $H=\gamma ^{-1/2}R_p=H^*=(Q/u_b)^{1/2}$, the gravity-driven flux (2.31) equals the input flux $Q$, indicating that $\mathcal {Q}=1$. We also note that, after the transition time $t>t^*$, the gravity-driven flux (2.31) continues to grow unbounded over time. Hence, this can satisfy neither the boundary condition at the origin (2.15), nor that at the nose (2.16). Losing the ability to impose the boundary condition at the nose indicates that a shock profile will form, over which the flux discontinuously drops to zero (this is a feature of such equations, e.g. see Whitham Reference Whitham2011) Hence, in terms of the classification of partial differential equations, the problem is hyperbolic, unlike the case of a two-dimensional gravity current due to a line source (Huppert & Woods Reference Huppert and Woods1995) which is parabolic.

Next, we consider conservation of mass to determine the model parameters. This is given by the volume integral

(2.32)\begin{equation} 2{\rm \pi}\phi \int_0^R r h\,\mathrm{d}r =Q t.\end{equation}

Hence, by inserting (2.29) into (2.32), we arrive at a relationship involving $\zeta _N$ and $\xi _0$, which is

(2.33)\begin{equation} \frac{2{\rm \pi} \zeta_N^2\xi_0}{({\rm \pi} \gamma)^{1/2}} \int_{\epsilon(t)}^{1} y (-\log y)^{1/2}\,\mathrm{d}y = 1,\end{equation}

where $\epsilon =R_p(t)/R(t)={O}(t^{-(b-a)})$. Noting that $\epsilon$ becomes vanishingly small at late times (since $b>a$), and noting that $\int _0^1 y (-\log y)^{1/2}\,\mathrm {d}y =({\rm \pi} /32)^{1/2}$, we can rearrange (2.33) to give

(2.34)\begin{equation} \zeta_N^2\xi_0 \approx {{(8\gamma)^{1/2}}}/{\rm \pi}.\end{equation}

A further final condition is required to fully determine the coefficients $\zeta _N$, $\xi _0$, and this is provided by matching the flux of the current between the inner and outer regions. Moving from the outer region to the origin the thickness tends to a finite but growing value $h=H(t)$ (with zero slope) and a balance is sustained at the matching radius $r=R_p$, so the gravity-driven flux becomes

(2.35)\begin{equation} -2{\rm \pi} u_b {r} {h} \frac{\partial {h}}{\partial {r}}= 2{\rm \pi} u_b \frac{R_p^3}{\gamma R} = Q,\quad\mathrm{at}\ r=R_p(t).\end{equation}

Hence, (2.35) suggests that

(2.36)\begin{equation} b=3a,\end{equation}

which together with (2.26), gives the power law values

(2.37a,b)\begin{equation} a=1/7,\quad b=3/7.\end{equation}

Furthermore, (2.35) implies that the coefficients satisfy

(2.38)\begin{equation} \xi_0^3/\zeta_N = \gamma/2{\rm \pi}.\end{equation}

Hence, the two equations (2.34) and (2.38), can be solved simultaneously to give

(2.39)$$\begin{gather} \xi_0=(2^{1/2}{\rm \pi}^3)^{{-}1/7}\gamma^{5/14}\approx 0.58\cdot \gamma^{5/14}, \end{gather}$$

(2.40)$$\begin{gather}\zeta_N=(2^{11/2}{\rm \pi}^{{-}2})^{1/7}\gamma^{1/14}\approx 1.24\cdot \gamma^{1/14}. \end{gather}$$

In addition, the thickness and radial scalings are given by

(2.41)$$\begin{gather} H=\frac{1}{\gamma^{1/7}}\left(\frac{Q^3}{2^{1/2}u_b^{2}{\rm \pi}^{3}}\right)^{1/7} \left(\frac{t}{\phi}\right)^{1/7}, \end{gather}$$

(2.42)$$\begin{gather}R=\gamma^{1/14}\left(\frac{2^{11/2}Q^2u_b}{{\rm \pi}^{2}} \right)^{1/7} \left(\frac{t}{\phi}\right)^{3/7}, \end{gather}$$

(2.43)$$\begin{gather}R_p=\gamma^{5/14}\left(\frac{Q^3}{2^{1/2}u_b^{2}{\rm \pi}^{3}}\right)^{1/7} \left(\frac{t}{\phi}\right)^{1/7}, \end{gather}$$

which now all reflect a dependence on the anisotropy $\gamma$. We note one further detail of this analysis. Whilst $H(t)$ (2.41) captures the overall thickness scaling of the gravity current, it does not reflect the precise value at the origin. In particular, if we insert the matching radius $r=R_p$ into the thickness profile in the outer region (2.29), we find that the thickness near the origin actually grows like

(2.44)\begin{equation} h(R_p(t),t)= \frac{R_p}{\gamma^{1/2}}\left(-\frac{1}{\rm \pi} \log \frac{R_p}{R}\right)^{1/2}\sim \left(\frac{2}{7{\rm \pi} \gamma} \log t\right)^{1/2} t^{1/7}. \end{equation}

However, since the relative growth of the logarithmic component is very slow, this only becomes appreciable at time scales $t/t^*$ larger than around $\exp ({7{\rm \pi} \gamma /2})\approx (5.96\times 10^{4})^\gamma$. To summarise, the scaling for $H(t)$ (2.41) provides the general behaviour of the thickness everywhere in the current except near the origin, and (2.44) must be used to capture the near-origin behaviour at very late times.

2.5. Comparison with numerical simulations

We compare our analytical predictions for the early-time, late-time and transition behaviour, with finite difference numerical simulations of axisymmetric Darcy flow in an anisotropic porous medium. Before discussing these results, we first briefly outline the numerical method. We model the flow using the Darcy equations (in cylindrical polar coordinates), which are

(2.45a,b)\begin{equation} \boldsymbol{u}={-}\frac{1}{\mu}\underline{\underline{\boldsymbol{k}}}\boldsymbol{\cdot} \boldsymbol{\nabla} \left( p + \rho_i g z \right),\quad \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u}=0,\quad \mathrm{in}\ \varOmega_i,\quad i=1,2,\end{equation}

where $\rho =\rho _1$ inside the injected fluid domain (which we denote $\varOmega _1$) and $\rho =\rho _2$ outside the injected fluid domain (which we denote $\varOmega _2$). The anisotropic permeability matrix $\boldsymbol {k}$ is accounted for by imposing different values, $k_H$ and $k_V$, in the radial and vertical directions.

The numerical solution at $r=0$ is not available due to the singular nature of the divergence in the continuity equation $\boldsymbol {\nabla }\boldsymbol {\cdot } \boldsymbol {u}=0$. Instead, the domain begins at a small but finite value of $r_0=0.05H^*$ close to the origin (we note that the solution is insensitive to the precise value $0.05$). The governing equations (2.45a,b) are accompanied by boundary conditions corresponding to impermeable/symmetric walls at $r=r_0$ and $z=0$, except for a constant input flux $Q$ injected over a small region near the origin, and far-field (Dirichlet) pressure conditions imposed at the boundaries of a large but finite domain. Further details on the boundary conditions and the dimensionless problem are given in Appendix A.

The interface between the domains $\varOmega _1$ and $\varOmega _2$ is given by a set of points $\partial \varOmega =\boldsymbol {r}(t)$, which is treated as a passive scalar and is therefore advected at the fluid velocity, such that

(2.46)\begin{equation} \frac{\mathrm{d}\boldsymbol{r}}{\mathrm{d}t}=\frac{1}{\phi}\boldsymbol{u}(\boldsymbol{r},t),\end{equation}

with initial conditions given by a small surface $\boldsymbol {r}(0)=\boldsymbol {r}_0$ close to the origin. Following the Lagrangian approach, the interface is discretised into points which are advected along streamlines. Therefore, at every time step the domains $\varOmega _1$ and $\varOmega _2$ are updated via interpolation over a gridded mesh of $150\times 150$ points. Likewise, the grid is stretched over time to account for the fact that the fluid region grows across several orders of magnitude. The dynamic equation for the fluid–fluid interface (2.46) is solved using an explicit Euler scheme with an adaptive time step, where at each time step the Darcy equations (2.45a,b), (with updated domains $\varOmega _i$) are solved using a second-order finite difference scheme. A small amount of Gaussian smoothing is applied to the interface at every time step to aid stability, but this is done under the constraint of mass conservation.

In figures 2 and 4 the evolution of the gravity current shape $h(r,t)$ is shown together with streamlines at different times. Both numerical and analytical (approximate) solutions are shown, and the evolution of the inner region is illustrated in red. At early times $t\ll t^*$, the shape of the current is well approximated by (2.8). Later, when $t\approx t^*$, the flow enters a transition regime in which gravity becomes appreciable in the majority of the current, except in the near-origin pressure-driven region (whose relative radial extent $R_p(t)/R(t)$ diminishes in size). Much later, when $t\gg t^*$, the shape of the gravity current converges to the approximate solution (2.29), and the streamlines away from the origin become horizontal. Meanwhile, streamlines within the inner region remain sloped, indicating that significant vertical velocities are present close to the origin.

Figure 4. (a–c) Numerical and analytical profiles for the gravity current shape at three late times $t/t^*=6\times (10^{1},10^{2},10^{3})$, corresponding to the last three curves in figure 2(b). Streamlines in cyan are calculated by integrating the numerical data using initial values close to the origin. The vertical scale is stretched by ${\times }10.5$ for illustration purposes.

We also plot the horizontal and vertical extents $R$, $H$, in figure 5. At early ($t/t^*<10^{-2}$) and late ($t/t^*>10^{2}$) times, the numerical solution matches the analytical scalings (2.6)–(2.7), (2.41)–(2.42) closely. At very late times, the slow logarithmic growth (2.44) can be observed as a slight deviation in the vertical extent $H$ away from the late-time scaling (2.41), which is consistent with our predictions. These plots are also displayed with linear scales in figure 9 in Appendix B for comparison.

Figure 5. Analytical, numerical and experimental data for the radial (a) and vertical (b) extents of the current, $R$, $H$. Numerical results are shown for both an isotropic medium $\gamma =1$ and an anisotropic medium $\gamma =10^4$, whereas experimental results are only for the isotropic case. Analytical early-and late-time scalings (2.6)–(2.7), (2.41)–(2.42), are shown with solid black lines and alternate scalings derived by Lyle et al. (Reference Lyle, Huppert, Hallworth, Bickle and Chadwick2005) ($R\sim t^{1/2},H\sim t^0$) are shown with dotted green lines. Best fit power law $H/H^*\propto (t/t^*)^{0.1537\pm 0.0176}$ for the late-time experimental thickness data is shown as an inset in (b).

To understand the role of the inner pressure-driven region, we have also calculated the maximum absolute value of the pressure deviation from hydrostatic $|p-p_{hyd}|_\infty$, normalised by the weight of the injected fluid $\Delta \rho g H$. This deviation does not fall to zero over time, but instead appears to tend to a constant value of around ${|p-p_{hyd}|_\infty /\Delta \rho g H\approx 3}$. This indicates that the flux is not well approximated by the gravity-driven component (2.30) everywhere within the current. As a result, the classical thin-film equation (2.13) and the similarity variables (2.17a,b) cannot be accurate in the late-time regime. In particular, the boundary condition at the origin feeding the gravity current with flux $Q$ cannot be supplied with a gravity-driven term of the form (2.30) as this would cause an infinite thickness. Instead, the only way for the thickness to remain finite is with a pressure-driven input flux at the origin, which explains why the pressure deviation from hydrostatic may never drop to zero over time.

We have also performed the same numerical calculation in the case of a two-dimensional injection from a line source (e.g. see figure 10 in Appendix B) and in this case the pressure deviation $|p-p_{hyd}|_\infty /\Delta \rho g H$ does drop to zero over time, indicating that (by contrast) the two-dimensional flux is well approximated by the corresponding gravity-driven component ($-u_b h \partial h/\partial x$) everywhere within the current at late times.

Since the power laws $t^{3/7}$ and $t^{1/2}$ are not hugely different, we have also tried comparing the numerical solution with the classical self-similar scalings (2.17a,b) proposed by Lyle et al. (Reference Lyle, Huppert, Hallworth, Bickle and Chadwick2005), as shown with dotted green lines in figure 5 and with solid green lines in figure 8 in Appendix B. Since the disparity between the different proposed scalings grows larger over time like ${\sim }t^{1/2-3/7}\sim t^{1/14}$ (in the case of radius) and like $\sim t^{0-1/7}\sim t^{-1/7}$ (in the case of thickness), after a time scale $t/t^*=10^6$ this corresponds to disparaging factors of $10^{6/14}\approx 2.68$ and $10^{-6/7}\approx 0.14$, respectively. Since these discrepancies lie outside the 5 % error observed in figure 5, this indicates that the scalings (2.41)–(2.42) are more accurate than the self-similar scalings (2.17a,b).

2.6. Comparison with experimental data

Next, we compare our results with experiments conducted in a porous medium tank. We use the same apparatus for our experiments as Lyle et al. (Reference Lyle, Huppert, Hallworth, Bickle and Chadwick2005) except that we use a peristaltic pump instead of a fluid reservoir to control the flow rate (for improved accuracy). A square Perspex tank with a $61\,{\rm cm}\times 61\,{\rm cm}$ base is filled with 0.3 cm glass Ballotini beads saturated with fresh water. Due to the difficulties associated with creating an anisotropic permeability field from these beads, our experiments are restricted to the isotropic case. Hence, the dependence on the anisotropy parameter $\gamma$ in the analytical scalings (2.6)–(2.7), (2.41)–(2.43) is checked against our numerical simulations, whilst the dependence on the remaining parameters (e.g. dimensional scalings and power laws) is checked against both numerical simulations and experiments (in the isotropic case).

Two cameras are positioned above and alongside the tank, pointing vertically downwards (plan view) and horizontally (side view), to capture both the radius and the thickness of the current. Salty water dyed with red food colouring is injected into the corner of the tank at different salt concentrations leading to a variety of density contrasts. By injecting into one corner of the tank the experiment represents one quarter of the full axisymmetric scenario whilst allowing for a larger range of radial data given the tank size.

A full list of the different input flow rates and reduced gravity values $g'=\Delta \rho g/\rho _w$ (where $\rho _w=1\,{\rm g}\,{\rm cm}^{-3}$ is the density of water) is given in table 1. The experiments at high flow rates $Q=17.36\unicode{x2013}31.62\,{\rm cm}^3\,{\rm s}^{-1}$ result in a current with an order ${O}(1)$ aspect ratio, with most of the injected fluid located in the interior region of the tank (away from the walls). Hence, the porosity and permeability are estimated using the values calculated for randomly packed beads, $\phi =0.37$ and $k=6.8 \times 10^{-9}\,{\rm m}^2$ (Lyle et al. Reference Lyle, Huppert, Hallworth, Bickle and Chadwick2005). However, the experiments conducted at low flow rates $Q=0.020\unicode{x2013}0.209\,{\rm cm}^3\,{\rm s}^{-1}$ result in a gravity current no thicker than 1–3 bead diameters ($z=d-3d$) and hence the flow is subject to wall effects. In these cases the porosity is estimated using the empirical formula derived by Ribeiro, Neto & Pinho (Reference Ribeiro, Neto and Pinho2010) accounting for wall effects, which is

(2.47)\begin{equation} \phi=0.373+0.917\exp({-0.824z/d})\approx0.55,\end{equation}

where we have used an average flow depth value $z=2d$, and the permeability is estimated using the Kozeny–Carman relationship

(2.48)\begin{equation} k=\frac{d^2\phi^3}{180(1-\phi)^2}\approx 4.36 \times10^{{-}8}\,{\rm m}^2. \end{equation}

It is worth noting how variations in the permeability (and hence the velocity) due to wall effects (2.47) may modify the flow. As discussed by Hinton & Woods (Reference Hinton and Woods2018), whilst such shear effects may alter the coefficients of the spreading rates, they do not modify the power laws, which motivates the use of different constant values here.

Table 1. List of experiments and corresponding parameter values.

It is also worth briefly discussing the appropriate size of a representative elementary volume (REV) for this medium, and consequently the possible limitations on modelling the flow using a homogenised set of equations such as Darcy's law. The form of (2.47) suggests that an appropriate length scale for the REV is simply $d$ (the bead diameter), which is the $e$-folding length scale for porosity changes. Hence, it is expected that the porosity and permeability may become more uncertain when the flow becomes smaller than one bead diameter in thickness (e.g. towards the outer radius of the gravity current at very small flow rates). However, in such cases upper and lower bound estimates can be calculated using (2.47), for example.

By varying the flow rate and reduced gravity, a variety of different transition time scalings are achieved, producing data in the range $t/t^*=10^{-3}\unicode{x2013}10^{5}$. We also complement our data with the experiments of Lyle et al. (Reference Lyle, Huppert, Hallworth, Bickle and Chadwick2005) which are in the intermediate range $t/t^*=0.08\unicode{x2013}22$. To reach such late-time values we inject using a very slow flow rate of $Q=0.02\,{\rm cm}^3\,{\rm s}^{-1}$ over one hour with photos taken at logarithmically spaced time intervals starting at 30 s (see figure 6). After post-processing the images, the radial and vertical extents of the current are extracted by dividing the pixels into binary values according to a threshold value (enabled by the colour contrast due to the dye).

Figure 6. (a,b) Close-up snapshots of the gravity current thickness near the origin during two injection experiments with different flow rates. Photos are taken over 30 min (a) and 1 h (b), and spaced logarithmically in time. (a) $Q = 0.209\,{\rm cm}^3\,{\rm s}^{-1}$, $H^* = 0.36$ cm; (b) $Q = 0.020\,{\rm cm}^3\,{\rm s}^{-1}$, $H^* = 0.11$ cm.

Data for the thickness and radius of the current are plotted together with numerical and analytical results in figure 5. The radius data $R(t)$ approximately follow the $3/7$ power law (2.42), however, one could argue that the $1/2$ power law (2.17a,b) is just as accurate within the error margins of the experimental data. More clarity is achieved by observing the variation of the thickness near the origin $H(t)$, as shown by the photos in figure 6. In particular, the thickness is clearly observed to increase beyond the transition value by as much as $H\approx 5.5 H^*$ over the course of the experiments, albeit with a very small growth rate. Indeed, without running the experiments for such a long time it would appear that the thickness has ceased growing entirely. By contrast, when observed at logarithmically spaced time intervals the thickness shows no signs of abating. We estimate that the error associated with measuring the thickness of the current is approximately one bead diameter, so we have added error bars to illustrate this.

Post-processing of the images to extract the interface position reveals that the thickness follows the 1/7 power law (2.41) to good approximation, as shown in figure 5. We have also tried using a least-squares minimisation that fits an arbitrary power law $H= A t^B$ to the combined experimental dataset (over four different flow rates) at late times $t>100 t^*$, as shown in the figure inset. This produces a power law of $B=0.1537\pm 0.0176$ which is close to $1/7\approx 0.1429$. Moreover, this slightly larger numerical value for the power law might hint towards the faster logarithmic growth predicted earlier.

It is not possible to discern the logarithmic correction to the thickness (2.44) quantitatively since this would require much longer time scales than we could achieve experimentally. Nor is it possible to accurately discern the convergence of the experimentally measured shape to the approximate solution (2.29) over long times, since the aspect ratio of the current is incredibly slender (around 200:1). To extend the experiments to larger values of $t/t^*$ is challenging in practice, since smaller flow rates than $Q=0.02\,{\rm cm}^3\,{\rm s}^{-1}$ would result in current thicknesses much smaller than a single bead size ($H^*\ll 0.3\,{\rm cm}$) that would be difficult to observe. Alternately one could run the experiments for a longer time, but this requires a much larger tank than we have available. For example to extend $t/t^*$ by a factor 100 would require a tank of dimensions $100^{3/7}\approx 7.2$ times larger (i.e. with a $439\,{\rm cm}\times 439\,{\rm cm}$ base). Nevertheless, our experiments clearly show an increasing thickness near the origin with the approximate 1/7 power law behaviour (2.41). Hence, by conservation of mass (2.24) it follows that the radius must grow according to a 3/7 power law, and the shape profile (2.29) must become a good approximation over long times (e.g. see analysis in § 2.4 and numerical solution in figures 2, 4).

It is important to note the possible effects of diffusion/dispersion over such long time scales. The length scales for molecular diffusion, fluid dispersion, and Taylor dispersion are approximately $\ell \sim (D t)^{1/2}$, $\ell \sim (d U t)^{1/2}$ and $\ell \sim (d^2 U^2 t / 48D)^{1/2}$, respectively, where $D$ is the molecular diffusivity of salt or dye (both similar values), $U$ is a characteristic velocity scale and $d=0.3\,{\rm cm}$ is the diameter of the Ballotini beads. As an estimate for the diffusivity we take $D=10^{-10}\,{\rm m}^2\,{\rm s}^{-1}$, and for the velocity we take the value at the top of the current $U\approx \mathrm {d}H/\mathrm {d}t$. Inputting these, we calculate diffusion/dispersion scales of $\ell \sim 0.06\,{\rm cm}$, $\ell \sim 0.13\,{\rm cm}$ and $\ell \sim 0.04\,{\rm cm}$, which are all considered small within the experimental context. It is clear from the images in figure 6 that the interface near the origin remains relatively sharp (although there are some light scattering effects due to the glass beads), indicating that the effects of dispersion are negligible there. Hence, this validates our approach of post-processing the images to extract the current thickness at late times.