1. Introduction
Among possible substitutes for hydrocarbon fuels, hydrogen has a high energy density and can be converted into heat or electricity without emitting 
$\mathrm {CO_2}$ (Andrews & Shabani Reference Andrews and Shabani2012). Although there is an obvious incentive to generate hydrogen from renewables, it is difficult to do so on a temporally consistent basis given the variability of e.g. wind forcing and solar radiation. As such, options for 
$\mathrm {H}_2$ storage must be pursued so that 
$\mathrm {H}_2$ produced in excess may be stored, then used when demand outstrips supply (Sainz-Garcia et al. Reference Sainz-Garcia, Abarca, Rubi and Grandia2017). An appealing option for 
$\mathrm {H}_2$ storage, particularly when considering large volumes of 
$\mathrm {H}_2$, may be underground storage (Panfilov Reference Panfilov2016). The technical and economic feasibility of underground 
$\mathrm {H}_2$ storage (UHS) has been studied in various locales e.g. Bulgaria, United States, United Kingdom, Poland, Spain and Turkey (Flesch et al. Reference Flesch, Pudlo, Albrecht, Jacob and Enzmann2018). Here, we consider UHS in the context of depleted reservoirs. In particular, we focus on depleted natural gas reservoirs, which avoid a possible contamination of 
$\mathrm {H}_2$ by the longer chain organic molecules present in depleted oil reservoirs.
Underground 
$\mathrm {H}_2$ storage in depleted natural gas reservoirs requires cushion gas, gas stored permanently in formation to maintain pressure for optimum injection and withdrawal. Although the cushion gas may be 
$\mathrm {H}_2$, it is more typically a dissimilar (heavier) gas such as 
$\mathrm {CO}_2$ or 
$\mathrm {N}_2$ (Feldmann et al. Reference Feldmann, Hagemann, Ganzer and Panfilov2016). Thus 
$\mathrm {H}_2$ injection into cushion gas may be associated with significant mixing, whether due to diffusion or dispersion. Mixing may be exacerbated by buoyancy effects, which result from the small size of the 
$\mathrm {H}_2$ molecule. In turn, 
$\mathrm {H}_2$ has a high mobility in formation and may therefore travel long lateral distances or else leak into adjoining stratigraphic layers (Lubon & Tarkowski Reference Lubon and Tarkowski2021). Leakage often arises from local fault(s), which act as pathways through otherwise impermeable layers – see e.g. Flett, Gurton & Taggart (Reference Flett, Gurton and Taggart2005). Complicating matters are the facts that (i) local faults may prove difficult to detect in surveys and (ii) monitoring injectate migration in UHS operations is non-trivial and expensive. There is a need, therefore, for tractable conceptual models of buoyancy-driven flow, drainage and dispersion that may inform key processes important to the techno-economic evaluation of UHS projects. Of course, such conceptual models might additionally consider effects such as viscous fingering, capillary effects, bio-geochemical reactions and, when residual liquid is present, capillarity, dissolution and chemical reaction. However, and in the interests of simplicity, we do not examine such additional effects here.
Another possible application of our work is to 
$\mathrm {CO}_2$ sequestration. In this related (and better-studied) problem, one likewise considers the eventual fate of a fluid that is injected at high pressure into a porous medium. Similar to UHS, the success of 
$\mathrm {CO}_2$ sequestration relies, in part, on considerations of the mixing (e.g. by dissolution cf. MacMinn et al. Reference MacMinn, Neufeld, Hesse and Huppert2012; Khan, Bharath & Flynn Reference Khan, Bharath and Flynn2022) between the injectate with the ambient fluid (i.e. brine) that occupies the pore space. However, 
$\mathrm {CO}_2$ sequestration flows are complicated by surface tension effects and the possibility of capillary trapping, which arise because of the relative immiscibility of 
$\mathrm {CO}_2$ and brine.
Buoyancy-driven flows in porous media are bookended by two canonical scenarios: a vertically ascending or descending plume and a horizontally propagating gravity current. Wide attention has been devoted to gravity current flow in porous media since the seminal work of Huppert & Woods (Reference Huppert and Woods1995), who studied the evolution of finite releases of fluid propagating in an expansive rectilinear porous medium. They considered that (dense) fluid moves under gravity along either horizontal or inclined boundaries and through a medium whose permeability is either uniform or else changes normal to the lower (impermeable) boundary. (For analytical convenience, many previous studies assume that the gravity current density is larger than that of the surrounding ambient. To be consistent with this earlier body of research, a similar assumption shall be adopted here. In this respect, the gravity currents to be described quantitatively in e.g. section 2 are ‘upside down’ relative to those expected in real UHS operations.) Huppert and Woods's analysis was extrapolated to axisymmetric geometries by Lyle et al. (Reference Lyle, Huppert, Hallworth, Bickle and Chadwick2005). Further extensions to Huppert & Woods's work have considered gravity current flow in shallow porous media (MacMinn et al. Reference MacMinn, Neufeld, Hesse and Huppert2012; Pegler, Huppert & Neufeld Reference Pegler, Huppert and Neufeld2014), through horizontally heterogeneous media (Zheng, Christov & Stone Reference Zheng, Christov and Stone2014) and through layered porous media that are either horizontal (Pritchard, Woods & Hogg Reference Pritchard, Woods and Hogg2001; Goda & Sato Reference Goda and Sato2011; Sahu & Flynn Reference Sahu and Flynn2017) or else make an angle to the horizontal (Bharath, Sahu & Flynn Reference Bharath, Sahu and Flynn2020). Studies have further considered non-Newtonian gravity currents (Ciriello et al. Reference Ciriello, Longo, Chiapponi and Federico2016) and gravity currents consisting of fluid having spatially variable densities (Pegler, Huppert & Neufeld Reference Pegler, Huppert and Neufeld2016).
Notable in the above summary of previous research are investigations involving a leakage of fluid from the gravity current underside. This leakage may be either localized (e.g. through a discrete fissure) or else distributed (e.g. into a lower layer of comparatively small permeability). Pritchard (Reference Pritchard2007) modelled gravity current flow in a porous medium with a series of line fissures in which drainage is due to the hydrostatic pressure exerted by the column of gravity current fluid situated directly above a particular fissure. Neufeld, Vella & Huppert (Reference Neufeld, Vella and Huppert2009) expanded this analysis by additionally considering the weight of the dense fluid inside and below the fissure. More recently, Gilmore et al. (Reference Gilmore, Sahu, Benham, Neufeld and Bickle2021) combined Neufeld et al.'s description with the plume solution of Sahu & Flynn (Reference Sahu and Flynn2015) to study flow in faults cross-cutting multiple aquifers. Meanwhile Avci (Reference Avci1994) studied local drainage in separated confined aquifers taking into account the effect of the injection pressure and the natural contrast of hydraulic head in two separated aquifers. Nordbotten, Celia & Bachu (Reference Nordbotten, Celia and Bachu2004) extended Avci's work for multiple abandoned leaky wells. The above studies mostly invoke a sharp interface approximation and so ignore mass transfer between gravity currents and the surrounding ambient fluid e.g. through dissolution, diffusion or dispersion. However, in miscible flow e.g. gas reservoir storage of 
$\mathrm {H}_2$, the numerical simulations of Feldmann et al. (Reference Feldmann, Hagemann, Ganzer and Panfilov2016) indicate that mixing between the injected and ambient fluids may be non-trivial. As we shall see, this feature becomes more prominent in the presence of draining.
Mixing in porous media involves dispersion and diffusion. Diffusion is a process driven at the molecular scale by concentration differences while dispersion is advection driven and is related to the macro-scale flow phenomena. The mixing that occurs between miscible fluids in porous medium flow depends on the Péclet number, which characterizes the importance of advection to diffusion. Mixing is due to diffusion for 
$Pe \ll 1$ and due to dispersion (transverse and longitudinal) for 
$Pe \gg 1$ (Delgado Reference Delgado2007). Studies that explore mixing in the context of buoyancy-driven porous medium flow include Szulczewski & Juanes (Reference Szulczewski and Juanes2013). They examined mixing due to diffusion for lock exchange flows consisting of two fluids in vertically confined permeable rock. Szulczewski & Juanes (Reference Szulczewski and Juanes2013) showed that, if a constant volume of dense fluid is released into light fluid, there is an evolution through the following regimes: (i) diffusion-dominated flow, (ii) slumping in which the interface between the two fluids is sharp and tilts in an S-shaped curve, (iii) slumping in which the interface remains sharp but changes from an S-shaped curve to a straight line, (iv) Taylor slumping where mixing increases due to Taylor dispersion at the aquifer scale and decelerates the flow and (v) late diffusion where, similar to (i), transport occurs primarily by diffusion. Hinton & Woods (Reference Hinton and Woods2018) modelled longitudinal shear dispersion due to a vertical gradient of permeability. They demonstrated that the pattern of longitudinal dispersion depends on a number of factors including (i) the viscosity (or mobility) ratio and (ii) the severity of the vertical permeability gradient. Huyakorn et al. (Reference Huyakorn, Andersen, Mercer and White1987) considered interfacial mixing associated with sea water intrusions into coastal aquifers. The study of Paster & Dagan (Reference Paster and Dagan2007) also applied boundary layer approximations and the von Kármán integral method to solve for the velocity-dependent transverse dispersion in sea water intrusions having a non-uniform flux field. Mixing in miscible gravity currents is studied directly by Sahu & Neufeld (Reference Sahu and Neufeld2020). In their study, the authors used the depth-averaged mass and buoyancy conservation equations to provide a theoretical model for porous medium gravity currents experiencing transverse dispersion only. Theoretically speaking, they determined that the gravity current buoyancy flux can be described by a self-similar solution. However, in contrast to the sharp interface case, the gravity current height and concentration are not self-similar.
A limitation of the study by Sahu & Neufeld (Reference Sahu and Neufeld2020) is that it ignored longitudinal dispersion. Furthermore, they considered a dispersed interface only and did not define the bulk interface separately from the dispersed interface. The gravity current nose position is therefore identical to that predicted by a sharp interface formulation. Sahu & Neufeld's model works well when the lower boundary is impermeable, however, when draining is allowed to occur (either locally or in a distributed fashion), experimental evidence from Sahu & Flynn (Reference Sahu and Flynn2015) (e.g. their figure 3e) and from Bharath et al. (Reference Bharath, Sahu and Flynn2020) shows that a significant separation may arise between the fronts for the bulk and dispersed phases of the gravity current.
Due to these shortcomings in the literature, we seek to provide a theoretical model for porous medium gravity current flow where the bulk and dispersed phases are accounted for separately using the depth-averaged mass and buoyancy conservation equations for each phase. Also important is to develop a complementary numerical model (e.g. using COMSOL Multiphysics) to validate our theoretical model. In § 2, we derive a theoretical model for gravity currents experiencing dispersion and local drainage. Section 3 describes the numerical simulations meant to corroborate model output. In § 4, we discuss results and compare the predictions of the theoretical model with those due to the numerical simulations. Section 5 illustrates the application of our theoretical model to UHS in depleted reservoirs. Finally, current work is summarized and ideas for future research are outlined in § 6.
2. Theoretical model
2.1. Governing equations
We consider gravity current flow due to a dense fluid injection of density 
$\rho _s$ along a punctured boundary that makes an angle 
$\theta$ with the horizontal, as depicted in figure 1. Simplifying assumptions are as follows: (i) initially, the porous medium is saturated with ambient fluid of uniform density 
$\rho _0$; (ii) the source fluid and ambient fluid are incompressible and also miscible i.e. capillary effects can be ignored both in the medium as well as in the fissure; (iii) the gravity current consists of a bulk phase and a dispersed phase both of which remain long and thin such that the gravity current flow is everywhere hydrostatic, i.e. the Dupuit approximation is applicable, and in addition, the depth of the ambient is much larger than the gravity current depth; (iv) consistent with the Boussinesq approximation, the dynamic viscosity, 
$\mu$, is independent of concentration and therefore the viscosity of the bulk and dispersed phases are assumed equal; (v) at least until the location of the isolated fissure, the leading edge of the dispersed phase remains close to that of the bulk phase such that negligible drainage of dispersed fluid occurs; (vi) 
$Pe \gg 1$, such that diffusion is ignored; and (vii) broadly analogous to the dissolution study of MacMinn et al. (Reference MacMinn, Neufeld, Hesse and Huppert2012), whatever mixing occurs along the bulk-dispersed boundary leaves, within the bulk phase, a core of fluid whose density remains 
$\rho _s$. This core of bulk fluid remains upstream of the dispersed phase and, in time, must extend downstream of the fissure.
Figure 1. Schematic of a leaky gravity current propagating along an inclined boundary with local drainage through an isolated fissure. The gravity current consists of bulk and dispersed phases. Variables 
$h_1$ (bulk phase height), 
$h_2$ (overall height), 
$u_1$ (bulk phase velocity), 
$u_2$ (dispersed phase velocity), 
$w_{e1}$ (entrainment velocity from bulk phase), 
$w_{e2}$ (entrainment velocity from ambient) and 
$\overline{c}_2$ (average concentration in dispersed phase) are functions of 
$x$ and 
$t$. Meanwhile, variables 
$x_{N_b}$ (bulk phase nose position) and 
$x_{N_d}$ (dispersed phase nose position) are functions of 
$t$ only.
Following Vella & Huppert (Reference Vella and Huppert2006), the coordinate system 
$(x,z)$ associated with the along- and cross-slope directions is obtained by a clockwise rotation of the natural coordinates 
$(X,Z)$ through an angle 
$\theta$. The origin of both coordinate systems is coincident with the isolated source, which is indicated by the red dot in figure 1. In the analysis to follow, we restrict our attention to 
$x \geq 0$.
If the gravity current experiences local drainage through a fissure situated at 
$x=x_f$ and having width 
$\xi$, the continuity equation as applied to the bulk phase reads
\begin{equation} \phi\frac{\partial h_1}{\partial t} +\frac{\partial }{\partial x}\int_0^{h_1} u_1\,\mathrm{d}z ={-} w_{e1} -w_{d} F(x,x_f,\xi), \end{equation}
where 
$\phi$ is the porosity, 
$h_1$ is the height of the bulk phase and 
$w_{e1}$ and 
$w_{d}$ are the Darcy velocities respectively accounting for entrainment from the bulk to the dispersed phase and drainage through the fissure. Meanwhile 
$F(x,x_f,\xi )$ is a boxcar function centred on the fissure, which is zero everywhere except within the interval 
$x_f-{\xi }/{2}< x< x_f+{\xi }/{2}$. Because the pressure is hydrostatic, the bulk phase velocity 
$u_1$ does not change significantly in a direction perpendicular to the bottom boundary. Thus, 
$u_1$ can be considered independent of 
$z$ in (2.1) (Happel & Brenner Reference Happel and Brenner1991). Accordingly, (2.1) may be simplified to
\begin{equation} \phi\frac{\partial h_1}{\partial t} + \frac{\partial }{\partial x}(u_1 h_1) ={-} w_{e1} -w_{d} F(x,x_f,\xi). \end{equation}
The solute concentration in the bulk phase is assumed to be equal to the source concentration, 
$c_s$; consequently, it is unnecessary to apply a solute conservation equation in the bulk phase.
The continuity equation for the dispersed phase is
\begin{equation} \phi\frac{\partial (h_2-h_1) }{\partial t} +\frac{\partial }{\partial x}\int_{h_1}^{h_2} u_2 \,\mathrm{d}z = w_{e1} + w_{e2},\end{equation}
in which 
$h_2 - h_1$ and 
$u_2$ are, respectively, the thickness and speed of the dispersed phase. (Note that, consistent with the caption to figure 1, 
$u_2$ is assumed independent of 
$z$). Finally, 
$w_{e2}$ is the entrainment velocity from the ambient to the dispersed phase. By simplifying and exploiting (2.2), the above result can be rewritten
\begin{equation} \phi \frac{\partial h_2 }{\partial t} + \frac{\partial }{\partial x} [u_2(h_2-h_1)] ={-}\frac{\partial }{\partial x}(u_1 h_1) + w_{e2}-w_{d} F(x,x_f,\xi) . \end{equation}Finally, solute conservation in the dispersed phase provides
\begin{equation} \phi\frac{\partial }{\partial t} \int_{h_1}^{h_2} c_2\,\mathrm{d}z+\frac{\partial }{\partial x}\int_{h_1}^{h_2} u_2 c_2\,\mathrm{d}z = w_{e1}c_s . \end{equation}
Here, 
$w_{e1}$ is defined only over the interval 
$0 \leq x \leq x_{N_b}$ (see figure 1). By defining the dispersed phase buoyancy as 
$b_2=\bar {c}_2(h_2-h_1)$, in which 
$\bar {c}_2$ is the 
$z$-averaged solute concentration in the dispersed phase, (2.5) can be further simplified to
\begin{equation} \phi\frac{\partial b_2}{\partial t} + \frac{\partial }{\partial x}(u_2 b_2) = w_{e1}c_s. \end{equation}
Because the bulk buoyancy, 
$b_1=c_sh_1$, equals the source buoyancy, solute conservation in the bulk phase is trivial.
Similar to the classical entrainment hypothesis that was proposed for flow in free jets by Ellison & Turner (Reference Ellison and Turner1959), we consider that the entrainment of ambient fluid is proportional to the gravity current characteristic velocity. Of greater relevance to buoyancy-driven flow in porous media, Sahu & Neufeld (Reference Sahu and Neufeld2020) also used a linear relationship between the entrainment and characteristic velocities in their study of dispersive gravity currents. Motivated by this latter work most especially, we define 
$w_{e1}=\varepsilon _1 u_1$ and 
$w_{e2}=\varepsilon _2 u_2$, where 
$\varepsilon _1$ and 
$\varepsilon _2$ are entrainment coefficients that account for the effects of dispersive mixing. Theoretically speaking, there is no reason that 
$\varepsilon _1$ and 
$\varepsilon _2$ have to be different. Therefore, we assume 
$\varepsilon _1=\varepsilon _2=\varepsilon$ so as to reduce the number of free parameters in our problem. (The preliminary accuracy of this assumption can be assessed in the context of the agreement between theory and numerical simulation to be presented later. A more detailed assessment of the relative magnitudes of 
$\varepsilon _1$ and 
$\varepsilon _2$ requires a dedicated study and so is left for future work.) On the other hand, and motivated by analogous studies of turbulent free gravity currents (Ellison & Turner Reference Ellison and Turner1959; Reeuwijk, Holzner & Caulfield Reference Reeuwijk, Holzner and Caulfield2019), we allow the possibility that 
$\varepsilon$ varies with the inclination angle of the bottom boundary, 
$\theta$. Such a dependence will be explored below.
Pressure in the bulk and dispersed phases is defined as
$$\begin{gather} p_1(x,z,t)=[\Delta \bar{\rho}_2 g h_2+ ( \Delta \rho_1-\Delta \bar{\rho}_2 )g h_1 -\rho_s g z]\cos \theta +\rho_0 g x\sin \theta + P_0 \quad 0 \leq z \leq h_1 , \end{gather}$$
$$\begin{gather}p_2(x,z,t)=[\Delta \bar{\rho}_2 g h_2 - \bar{\rho}_2g z]\cos \theta +\rho_0 g x \sin \theta +P_0 \quad h_1 \leq z \leq h_2, \end{gather}$$
in which g is the gravitational acceleration, 
$P_0$ is a reference pressure evaluated at 
$x=z=0$ outside of the gravity current, 
$\rho _s$ is the source (or bulk) fluid density and 
$\bar {\rho }_2$ is the 
$z$-averaged density in the dispersed phase. Moreover, 
$\Delta \rho _1=\rho _0\beta c_s$ is the density difference between the bulk and ambient phases and 
$\Delta \bar {\rho }_2=\rho _0\beta \bar {c}_2$ is the corresponding density difference for the dispersed phase, averaged over 
$z$. In deriving the previous expressions for 
$\Delta \rho _1$ and 
$\Delta \bar {\rho }_2$, reference is made to a linear equation of state of the form 
$\rho = \rho _0(1+\beta c)$ in which 
$\rho _0$ is a reference density and 
$\beta$ is the solute contraction coefficient. There is, in fact, a more subtle assumption associated with the derivation of (2.8) namely that vertical variation of 
$\rho _2$ small enough such that
\begin{equation} \frac{1}{h_2-z}\int_z^{h_2} \rho_2 \, \mathrm{d}z \simeq \frac{1}{h_2-h_1}\int_{h_1}^{h_2} \rho_2 \, \mathrm{d}z \equiv \bar{\rho}_2 .\end{equation}Stated differently, the above assumption suggests a dispersed phase hydrostatic balance of the form
\begin{equation} \frac{\partial p_2}{\partial z}={-} \bar{\rho}_2g,\end{equation}rather than
\begin{equation} \frac{\partial p_2}{\partial z}={-} \rho_2g .\end{equation}Our assumption that (2.10) and (2.11) are approximately equal is, of course, consistent with the principle of ignoring the vertical variation of concentration and of velocity in the dispersed phase.
By combining (2.7) and (2.8) with Darcy's law, i.e.
\begin{equation} \boldsymbol{V} ={-}\frac{k}{\mu}(\boldsymbol{\nabla} p - \rho \boldsymbol{g}), \end{equation}
where 
$\boldsymbol {V}$ is the Darcy flux, 
$\mu$ is the dynamic viscosity and 
$k$ is the (assumed constant) medium permeability, the calculation steps of Appendix A suggest that
$$\begin{gather} u_1(x,t)={-}\frac{kg\beta}{\nu}\left[\frac{\partial b_2}{\partial x}\cos \theta +c_s\left(\frac{\partial h_1}{\partial x}\cos \theta-\sin \theta\right)\right] , \end{gather}$$
$$\begin{gather}u_2(x,t)={-}\frac{kg\beta}{\nu}\left[\frac{\partial (\bar{c}_2 h_2)}{\partial x}\cos \theta-\bar{c}_2\sin \theta\right] \equiv{-}\frac{kg\beta}{\nu}\left[\frac{\partial }{\partial x} \left( \frac{b_2h_2}{h_2-h_1} \right)\cos \theta-\bar{c}_2\sin \theta\right]. \end{gather}$$
Here, 
$\nu = \mu /\rho _0$ is the kinematic viscosity.
If we assume drainage to be hydrostatically driven through a fissure having permeability 
$k_f$, width 
$\xi$ and depth 
$l$, application of Darcy's law (2.12) similar to Neufeld et al. (Reference Neufeld, Vella and Huppert2009) yields the following expression for the drainage velocity:
\begin{equation} w_{d}(x,t)=\frac{k_f g\beta }{\nu} \left(\frac{c_sh_1+b_2}{l}+c_s\right)\cos \theta; \end{equation}
(see Appendix B). Upon substituting ((2.13)–(2.15)) and the expressions for the entrainment velocities 
$w_{e1}$ and 
$w_{e2}$ into (2.2), (2.4) and (2.6), the following modified governing equations result:
\begin{gather} \phi\frac{\partial h_1}{\partial t} +\frac{k g\beta }{\nu}\frac{\partial (h_1U)}{\partial x} ={-}\varepsilon \frac{k g\beta }{\nu}U -\frac{k_f g\beta }{\nu} \left(\frac{c_sh_1+b_2}{l}+c_s\right)\cos \theta\times F(x,x_f,\xi), \end{gather}
\begin{gather} \phi\frac{\partial h_2}{\partial t} -\frac{k g\beta }{\nu} \frac{\partial }{\partial x} \left[(h_2-h_1)\left(\frac{\partial \varPsi }{\partial x}-C \right) -h_1 U \right] \nonumber\\ ={-} \varepsilon \frac{k g\beta }{\nu}\left(\frac{\partial \varPsi }{\partial x}-C\right) -\frac{k_f g\beta }{\nu} \left(\frac{c_sh_1+b_2}{l}+c_s\right)\cos \theta\times F(x,x_f,\xi), \end{gather}
\begin{gather} \phi\frac{\partial b_2}{\partial t} -\frac{k g\beta }{\nu} \frac{\partial }{\partial x} \left[b_2 \left(\frac{\partial \varPsi}{\partial x}-C \right)\right] = \varepsilon \frac{k g\beta c_s }{\nu}U . \end{gather}In the above equations, and for the sake of notational economy, we have introduced the following symbols:
$$\begin{gather} U={-}\left(\frac{\partial b_2}{\partial x} +c_s\frac{\partial h_1}{\partial x}\right)\cos \theta+c_s \sin \theta, \end{gather}$$
$$\begin{gather}\varPsi=\frac {b_2 h_2}{h_2-h_1}\cos \theta, \end{gather}$$
$$\begin{gather}C=\frac {b_2 }{h_2-h_1}\sin \theta. \end{gather}$$
The variables 
$U$, 
$\varPsi$ and 
$C$ are introduced only to simplify the notation; we do not regard these variables as having a noteworthy physical significance. Equations ((2.16)–(2.18)) contain three unknowns, namely 
$h_1$, 
$h_2$ and 
$b_2$. The equations are solved with the boundary conditions listed below.
2.2. Boundary conditions
Boundary conditions for ((2.16)–(2.18)) are
$$\begin{gather} -\frac{kg\beta}{\nu}\left[\left(\frac{\partial b_2}{\partial x} +c_s\frac{\partial h_1}{\partial x}\right)h_1\cos \theta -c_s h_1 \sin \theta \right]_{0}=q_s, \quad h_1\vert_{x_{N_b}} =0, \end{gather}$$
$$\begin{gather}h_2\vert_{0} =h_1\vert_{0}, \quad h_2\vert_{x_{N_d}} =0, \end{gather}$$
$$\begin{gather}b_2\vert_{0} =0 , \quad b_2\vert_{x_{N_d}} =0 . \end{gather}$$
Here, 
$x_{N_b}$ and 
$x_{N_d}$ are the bulk and dispersed nose positions, respectively as indicated in figure 1. Note that (2.22a) represents the influx boundary condition, 
$q_s=( u_1h_1)_{0}$, at the source. It is not required to take 
$b_1\vert _{x_{N_b}} = 0$ because the source concentration is fixed (and finite) so 
$b_1\vert _{x_{N_b}} = 0$ is automatically satisfied by (2.22b). From boundary condition (2.22c), we assume that the thickness, 
$h_2-h_1$, of the dispersed phase is zero at 
$x=0$; we investigate the validity of this assumption below. Finally, the expressions of global volume balance in the bulk phase and the expression of global buoyancy/solute balance for the combined bulk and dispersed phases can be written as
$$\begin{gather} q_st =\phi \int_{0}^{x_{N_b}} h_1\,\mathrm{d} x - \xi\int_{0}^{t} w_d\,\mathrm{d}t -\int_{0}^{t} \int_{0}^{x_{N_b}} w_{e1}\,\mathrm{d} x\,\mathrm{d}t , \end{gather}$$
$$\begin{gather}q_s c_st = \phi c_s\int_{0}^{x_{N_b}} h_1\,\mathrm{d} x - \xi c_s\int_{0}^{t} w_d\,\mathrm{d}t + \phi\int_{0}^{x_{N_d}} (h_2-h_1)\bar{c}_2\,\mathrm{d} x . \end{gather}$$
The first term on the right-hand side of (2.23a) represents the volume of the bulk phase fluid, the second term represents the volume of bulk fluid drained through the fissure and the third term represents the volume of fluid lost by the bulk phase to the dispersed phase. The first term on right-hand side of (2.23b) represents the buoyancy in the bulk phase, the second term represents the buoyancy lost by the bulk phase due to fissure drainage and the third term represents buoyancy in the dispersed phase where we have assumed implicitly that 
$h_1 = 0$ for 
$x_{N_b} < x < x_{N_d}$.
2.3. Non-dimensional governing equations
Consistent with Neufeld et al. (Reference Neufeld, Vella and Huppert2009), we respectively define the characteristic downdip length scale, the characteristic across-dip length scale and the characteristic time scale as
\begin{equation} X = x_f, \quad H = \left(\frac{x_f\,q_s\phi\nu}{kg^\prime_s}\right)^{1/2}, \quad \mbox{and} \quad T = \left(\frac{x_f^3\phi\nu}{kg^\prime_sq_s}\right)^{1/2}, \end{equation}
where 
$g^\prime _s=g^\prime _1= g\beta c_s$ is the source reduced gravity and 
$x_f$ denotes the fissure position – see figure 1. Thus do we define the following non-dimensional (starred) variables:
\begin{align} c_2^* = \frac{\bar{c}_2}{c_s}, \quad x^* = \frac{x}{X}, \quad \xi^* = \frac{\xi}{X}, \quad h_1^* = \frac{h_1}{H}, \quad h_2^* = \frac{h_2}{H}, \quad l^* = \frac{l}{H}, \quad t^* = \frac{t}{T}. \end{align} Neufeld et al. (Reference Neufeld, Vella and Huppert2009) defined a parameter to characterize the drainage through an isolated fissure of width 
$\xi$. With reference to this parameter and their (2.12), we define, for the flow depicted in figure 1, an upstream flow parameter 
$\varGamma$ and a fissure permeability ratio 
$K$ as
$$\begin{gather} \varGamma = \frac{\frac{x_f h_0}{q_s}}{\frac{ x_f}{u_b}} \frac{x_f}{h_0} =\frac{u_b x_f}{q_s} \equiv \frac{kg^\prime_s x_f}{\phi \nu q_s} , \end{gather}$$
$$\begin{gather}K = \frac{k_f}{k} , \end{gather}$$
respectively. In (2.26a), 
$h_0$ is the height of the gravity current at the source and 
$u_b={kg^\prime _s}/{\phi \nu }$ is the buoyancy velocity associated with a source concentration 
$c_s$. There are a variety of different ways to interpret the upstream flow parameter 
$\varGamma$. Firstly, 
$\varGamma$ can be thought of as the analogue of the Richardson number because its definition includes the ratio of the time, 
$t_s = {x_f h_0}/{q_s}$, for fluid to flow from the source to the fissure based on the source volume flux, to the time, 
$t_f ={x_f}/{u_b}$, for fluid to flow from the source to the fissure based on the source reduced gravity. Keeping with a ratio of time scales, 
$\varGamma$ can also be interpreted as a ratio including 
$t_f$ and 
$t_b = {q_s}/{u_b^2}$, which is a characteristic time of the flow. Finally, 
$\varGamma (=x_f/({q_s}/{u_b}))$ can be thought of as the ratio of the fissure distance, 
$x_f$, to the flow thickness, 
${q_s}/{u_b}$. As the above definitions of 
$\varGamma$ make clear, the larger the value of 
$\varGamma$, the more the gravity current flow is influenced by its density contrast with the ambient. The fissure permeability ratio 
$K(<1)$ corresponds to the ratio of the flow resistance through the porous medium to the drainage resistance through the fissure. Therefore, the larger 
$K$, the greater the volume of fluid that can drain through the fissure for a fixed depth of gravity current.
Applying the above definitions, ((2.16)–(2.18)) may be rewritten in non-dimensional form as
\begin{equation} \frac{\partial h_1^*}{\partial t^*} +\frac{\partial (h_1^*U^*)}{\partial x^*} ={-}\varepsilon \varGamma^{1/2} U^* -K\,\varGamma \left(\frac{h_1^*+b_2^*}{l^*}+1\right)\cos \theta\times F(x^*,1,\xi^*),\end{equation}
\begin{align} &\frac{\partial h_2^*}{\partial t^*} -\frac{\partial }{\partial x^*} \left[(h_2^*-h_1^*)\left(\frac{\partial \varPsi^* }{\partial x^*}-C^* \right) -h_1^* U^* \right] \nonumber\\ &\quad ={-} \varepsilon \,\varGamma^{1/2}\left(\frac{\partial \varPsi^* }{\partial x^*}-C^*\right) -K \varGamma \left(\frac{h_1^*+b_2^*}{l^*}+1\right)\cos \theta\times F(x^*,1,\xi^*), \end{align}
\begin{equation} \frac{\partial b_2^*}{\partial t^*} -\frac{\partial }{\partial x^*} \left[b_2^* \left(\frac{\partial \varPsi^*}{\partial x^*}-C^* \right)\right] = \varepsilon \varGamma^{1/2} U^* . \end{equation}Here,
$$\begin{gather} b_2^*=c_2^* (h_2^*-h_1^*), \end{gather}$$
$$\begin{gather}U^*={-}\left(\frac{\partial b_2^*}{\partial x^*} +\frac{\partial h_1^*}{\partial x^*}\right) \cos \theta+\varGamma^{1/2} \sin \theta, \end{gather}$$
$$\begin{gather}\varPsi^*=\frac {b_2^* h_2^*}{h_2^*-h_1^*}\cos \theta, \end{gather}$$
$$\begin{gather}C^*=\varGamma^{1/2}\frac {b_2^* }{h_2^*-h_1^*}\sin \theta . \end{gather}$$Also, the boundary conditions (2.22) now read
$$\begin{gather} \phi\left[\left(\frac{\partial b_2^*}{\partial x^*} +\frac{\partial h_1^*}{\partial x^*}\right)h_1^*\cos \theta -\varGamma^{1/2} h_1^* \sin \theta \right]_{0}={-}1, \quad h_1^*\vert_{x_{N_b}^*} =0, \end{gather}$$
$$\begin{gather}h_2^*\vert_{0} =h_1^*\vert_{0}, \quad h_2^*\vert_{x_{N_d}^*} =0, \end{gather}$$
$$\begin{gather}b_2^*\vert_{0} =0 , \quad b_2^*\vert_{x_{N_d}^*} =0 . \end{gather}$$
Assuming a fixed value for the medium porosity, 
$\phi$, there are five dynamically significant dimensionless groups in ((2.27)–(2.34)), namely 
$\varGamma$, 
$K$, 
$\xi ^*$, 
$l^*$ and 
$\theta$. These dimensionless groups characterize the fluid, medium and fissure properties.
The governing equations are solved using an explicit finite difference algorithm where spatial derivatives are discretized using backward finite differences because the source is situated on the upstream side (see Appendix C for more details). Sample results are shown in figure 2 for 
$\theta =0^{\circ }$ and for 
$\theta =5^{\circ }$. Because bulk fluid drains through the fissure but not so dispersed fluid, significant separation of the bulk and dispersed interfaces occurs only downstream of 
$x^* = 1$. Figure 2 illustrates a sharp change in the leading edge profile of the dispersed phase, especially at late times. The slope of this leading edge is set by a balance between the advection and dispersion. When, as is the case here, advection dominates, the nose of the dispersed phase is expected to change abruptly (though not discontinuously) as 
$x^* \to x_{N_d}^*$.
Figure 2. Theoretical predictions showing, for different times, gravity current profiles for (a) 
$\theta =0^\circ$ and (b) 
$\theta =5^\circ$. The thick line represents the bulk interface and the thin line represents the dispersed interface. The location of the fissure is as indicated. Here, 
$\varGamma =35$, 
$K=0.5$, 
$l^*=0.79$ and 
$\xi ^*=0.04$. As we will justify in § 3.5, we consider 
$\varepsilon = 0.0125$ when 
$\theta = 0^\circ$ and 
$\varepsilon = 0.0086$ when 
$\theta = 5^\circ$.
Accordingly, we focus on the dispersed phase and its fraction, relative to the gravity current as a whole, of buoyancy (per unit box width) and of area (volume per unit box width). In symbols, these quantities are denoted as 
$B$ and 
$A$, respectively. In performing the requisite calculations, we first evaluate the area enclosed by the bulk interface (the thick lines in figure 2) and by the dispersed interface (the thin lines in figure 2). Areas are calculated from
\begin{equation} A^*_{bulk} = \int_0^{x^*_{N_b}} h_1^*\,\mathrm{d} x^* \quad \mbox{and} \quad A^*_{disp} = \int_0^{x^*_{N_d}} (h_2^*-h_1^*)\,\mathrm{d} x^*. \end{equation}
The dispersed area fraction 
$\tilde {A}^*_{disp}$ is found from
\begin{equation} \tilde{A}^*_{disp}=\frac{A^*_{disp}}{A^*_{bulk}+A^*_{disp}} .\end{equation}Buoyancies in the bulk and dispersed phases are respectively determined from
\begin{equation} B^*_{bulk} = \int_0^{x^*_{N_b}} h_1^*\,\mathrm{d} x^*= A_{bulk}^* \quad \mbox{and} \quad B^*_{disp} = \int_0^{x^*_{N_d}} b^*_2\,\mathrm{d} x^*. \end{equation}
The quantities 
$B^*_{bulk}$ and 
$A^*_{bulk}$ are equal because concentrations are scaled by 
$c_s$. Similar to (2.36), the dispersed buoyancy fraction is given by
\begin{equation} \tilde{B}^*_{disp}=\frac{B^*_{disp}}{B^*_{bulk}+B^*_{disp}} .\end{equation} Figures 3(a) and 3(b) show time series of 
$\tilde {A}^*_{disp}$ and 
$\tilde {B}^*_{disp}$ for 
$\theta =0^{\circ }$ and 
$\theta =5^{\circ }$. As these figures make clear, dispersion increases when the bottom boundary is inclined and more solute will then reside in the dispersed phase. This phenomenon can be understood with reference to the significantly larger characteristic value for 
$u_1$ measured in the case of the sloping boundary, i.e. at 
$t^*=60$, the magnitude of 
$u_1$ when 
$\theta =5^{\circ }$ is approximately 30 % larger than that for 
$\theta =0^{\circ }$. A parametric study that more carefully documents the impact of the non-dimensional parameters 
$\theta$, 
$\varGamma$, 
$K$, 
$\xi ^*$ and 
$l^*$ on the evolution of the gravity current is included in § 4 below.
Figure 3. Percentage of the gravity current (a) volume and (b) buoyancy that remains in the dispersed phase for 
$\theta =0^\circ$ with 
$\varepsilon =0.0125$ and 
$\theta =5^\circ$ with 
$\varepsilon =0.0086$. Here, consistent with figure 2, 
$\varGamma =35$, 
$K=0.5$, 
$l^*=0.79$ and 
$\xi ^*=0.04$.
3. Numerical investigation
COMSOL simulations were performed so as to illustrate the effects of dispersive mixing in gravity currents within porous media and to assess our mathematical model in various scenarios. COMSOL utilizes the finite element method to discretize the governing equations (given below).
Our COMSOL model is validated in two complementary ways. First, we model the flow of a porous medium gravity current along an impermeable boundary (in which dispersion is comparatively small) and thereby demonstrate excellent agreement with the sharp interface solution of Huppert & Woods (Reference Huppert and Woods1995). Second, we confirm that our COMSOL model correctly predicts the degree of dispersion in a scenario where fluid density differences are absent, i.e. the scalar is passive rather than active. More specifically, we model the mixing of two miscible fluids in a long capillary tube. As considered by Bear (Reference Bear1972), § 10.6, the up- and downstream fluids have initial solute concentrations of 0 and 
$c_s$, respectively, but mix as a result of dispersion. Here, again, we observe excellent agreement between the theoretical solution and the corresponding COMSOL-based numerical result.
3.1. Numerical set-up
Simulations are conducted in a two-dimensional rectangular box 400 cm long and 25 cm deep filled with a porous medium saturated with water having a density 
$\rho _0=0.998 \ \mbox {g cm}^{-3}$. The medium porosity is 
$\phi =0.38$ based on the assumption of a random close packing of beads (Happel & Brenner Reference Happel and Brenner1991). The permeability is 
$2.18\times 10^{-4}\ \mbox {cm}^2$ and, by inverting Rumpf 
$\&$ Gupte's relation (Rumpf & Gupte Reference Rumpf and Gupte1971),
\begin{equation} k=\frac{1}{5.6}d_p^2 \phi^{5.5},\end{equation}
we determine that the equivalent bead diameter measures 
$d_p=5\ \mbox {mm}$, which is broadly consistent with the related experiments of Sahu & Flynn (Reference Sahu and Flynn2017) and Bharath et al. (Reference Bharath, Sahu and Flynn2020) for which the Reynolds number is of the order of 0.3. The salt water of fixed concentration is discharged at a constant rate from a source located in the bottom-left corner of the numerical domain – see figure 4. The source has a vertical expanse of 1 cm; broadly comparable to Neufeld et al. (Reference Neufeld, Vella and Huppert2009) the fissure is situated at a distance of 
$x_f=7.5$ cm from the source. At this location, deviation from a hydrostatic pressure gradient is small.
Figure 4. Schematic of the numerical set-up.
Two different COMSOL physics interfaces are used, i.e.
(i) The Darcy's law interface is used to model fluid flow within the porous medium specifically by solving the following mass and momentum equations:
(3.2a)
$$\begin{gather} \frac{\partial u}{\partial x}+\frac{\partial w}{\partial z}=0, \end{gather}$$(3.2b)$$\begin{gather}\frac{1}{\rho_0}\frac{\partial p}{\partial x}+\frac{\nu}{k}u ={-}\frac{\rho}{\rho_0}g \sin \theta, \end{gather}$$(3.2c)$$\begin{gather}\frac{1}{\rho_0}\frac{\partial p}{\partial z}+\frac{\nu}{k}w=\frac{\rho}{\rho_0}g \cos \theta . \end{gather}$$(ii) Solute transport is modelled using the transport of diluted species in porous media interface where the underlying equation to be solved reads
(3.3)
\begin{align} \phi\frac{\partial
c}{\partial t}+ u\frac{\partial c}{\partial
x}+w\frac{\partial c}{\partial
z}=\phi\left[\frac{\partial}{\partial x}
\left({\mathsf{D}}_{xx}\frac{\partial c}{\partial
x}+{\mathsf{D}}_{xz}\frac{\partial c}{\partial
z}\right)+\frac{\partial}{\partial z}
\left({\mathsf{D}}_{xz}\frac{\partial c}{\partial
x}+{\mathsf{D}}_{zz}\frac{\partial c}{\partial z}\right)\right] .
\end{align}
Here, 
$c$ is concentration and 
${\mathsf{D}}_{ij}$ is a dispersion coefficient. The dispersion tensor 
${\mathsf{D}}_{ij}$ is defined as
\begin{equation} {\mathsf{D}}_{ij}=D_{mol} +{\mathsf{a}}_{ijlm}\frac{V_lV_m}{|\boldsymbol{V}|} . \end{equation}
Here, 
$D_{mol}$ is the coefficient of molecular diffusion, 
$V_{()}$ is a component of the velocity whose overall magnitude is given by 
$|\boldsymbol {V}|$ and 
${\mathsf{a}}_{ijlm}$ is the geometrical dispersivity of the porous medium. Scheidegger (Reference Scheidegger1961) showed that 
${\mathsf{a}}_{ijlm}$ is a sparse matrix in which only terms 
${\mathsf{a}}_{1111}={\mathsf{a}}_{2222}=a_L$, 
${\mathsf{a}}_{1122}={\mathsf{a}}_{2211}=a_T$ and 
${\mathsf{a}}_{1212}={\mathsf{a}}_{2121}={\mathsf{a}}_{1221}={\mathsf{a}}_{2112}=\frac {1}{2}(a_L-a_T)$ are non-zero for two-dimensional porous medium flow. The dispersion coefficients for such a two-dimensional flow can therefore be defined as follows:
$$\begin{gather} {\mathsf{D}}_{xx}=D_{mol}+a_L\frac{u^2}{|\boldsymbol{V}|}+a_T\frac{w^2}{|\boldsymbol{V}|}, \end{gather}$$
$$\begin{gather}{\mathsf{D}}_{zz}=D_{mol}+a_L\frac{w^2}{|\boldsymbol{V}|}+a_T\frac{u^2}{|\boldsymbol{V}|}, \end{gather}$$
$$\begin{gather}{\mathsf{D}}_{xz}=D_{mol}+(a_L-a_T)\frac{|uw|}{|\boldsymbol{V}|} . \end{gather}$$
The variables 
$a_L$ and 
$a_T$ are called the longitudinal dispersivity and the transverse dispersivity, respectively. The dispersivities do not assume universal values and are, instead, resolved by curve fitting relevant experimental data corresponding to various regions of the parameter space e.g. as defined by the Schmidt (
$Sc$) and Péclet (
$Pe$) numbers. Notwithstanding this complication, we can adapt 
$(3)$ and 
$(4)$ of Delgado (Reference Delgado2007) to derive reasonable predictions for these two parameters for the Péclet numbers relevant to our flow. For mathematical convenience in practical applications and as suggested by Sahimi (Reference Sahimi2011), the power of 
$Pe$ in 
$(3)$ of Delgado (Reference Delgado2007) is considered to be unity for 
$5< Pe <300$. Therefore, the longitudinal and transverse dispersivities in this region are
$$\begin{gather} a_L = 0.5 d_p , \end{gather}$$
$$\begin{gather}a_T = 0.025 d_p , \end{gather}$$
respectively. For 
$300< Pe <10^5$, the relevant equations are
$$\begin{gather} a_L = (1.8 \pm 0.4) d_p , \end{gather}$$
$$\begin{gather}a_T = 0.025 d_p . \end{gather}$$
In this work, we take 
$a_L$ to be 
$1.8\,d_p$ for 
$300< Pe <10^5$.
3.2. Initial conditions
Consistent with the theory of § 2, we assume the porous medium is saturated with quiescent fresh water at 
$t=0$. The initial pressure distribution is therefore hydrostatic and the initial solute concentration is everywhere zero. The source concentration is related to the source density, 
$\rho _s$, via the equation of state, i.e.
\begin{equation} c_s=\frac{\rho_s-\rho_0}{\rho_0\beta} .\end{equation}
Here, 
$\beta$ is considered constant so the effects of temperature and pressure are ignored. According to the study of Millero & Poisson (Reference Millero and Poisson1981), 
$\beta$ is 0.82 
$\mbox {cm}^{3}\ \mbox {g}^{-1}$. After defining the source concentration for desired 
$\rho _s$, the linear equation of state 
$\rho =\rho _0(1+\beta c)$ is used to relate the density to the concentration field 
$c$.
3.3. Meshing and solver
Equations (3.2) and (3.3) are discretized using an unstructured triangular mesh. As shown in figure 4, grid refinement is applied in the vicinity of the source and of the fissure because these are regions of significant velocity shear. Figure 5 shows, for 
$t^*=40$, 
$\tilde {A}^*_{disp}$ for various grid sizes from coarse to fine. The dispersed phase area fraction is sufficiently close to its asymptotic value when the grid is comprised 
$81\, 715 \simeq 10^{4.91}$ triangles; at and beyond this point, we deem the numerical results to be grid independent.
Figure 5. Numerically determined estimates for the dispersed phase area fraction for different grid sizes. Here, consistent with figure 2, 
$\theta =0^\circ$, 
$\varGamma =35$, 
$K=0.5$, 
$l^*=0.79$ and 
$\xi ^*=0.04$.
To discretize the equations in space, cubic shape functions are chosen for (3.2), whereas quadratic shape functions are selected for (3.3). For this latter equation, a third-order implicit backward differentiation formula is applied such that
\begin{align} &\phi
\frac{11c^{n+1}-18c^n+9c^{n-1}-2c^{n-2}}{6\Delta t}
\nonumber\\ &\quad =\left[\phi \frac{\partial}{\partial x}
\left({\mathsf{D}}_{xx}\frac{\partial c}{\partial
x}+{\mathsf{D}}_{xz}\frac{\partial c}{\partial z}\right)+\phi
\frac{\partial}{\partial z} \left({\mathsf{D}}_{xz}\frac{\partial
c}{\partial x}+{\mathsf{D}}_{zz}\frac{\partial c}{\partial
z}\right)-u\frac{\partial c}{\partial x}-w\frac{\partial
c}{\partial z}\right]^{n+1} ,
\end{align}
where 
$n$ is the time increment.
We implement a two-step sequential method to solve (3.2) and (3.3) by using a two-step segregated solver within COMSOL. In the first step, (3.2a–c) are solved by considering the fluid density, 
$\rho$, as known. Then, in the second step, velocities calculated in step one are applied to solve (3.3) for concentration. Thus the Darcy and species equations are solved in sequence at each time step until convergence is achieved. In this work, we consider a relative convergence tolerance of 0.001.
3.4. Qualitative observations (horizontal bottom boundary)
One of the key assumptions applied in the model of § 2.1 is that there persists a bulk gravity current within which the solute concentration is effectively the same as the source concentration, 
$c_s$. The numerical simulations afford us the opportunity to test the validity of this assumption in different regions of the parameter space. To this end, numerical simulations indicate that the bulk phase concentration decreases downstream of the fissure and is less than 
$c_s$ due to the dispersion that arises in conjunction with drainage. This discrepancy between the bulk phase concentration and the source concentration can, for sufficiently vigorous dispersive mixing, affect the accuracy of our theoretical model. It is necessary, therefore, to estimate the parametric regime where such vigorous dispersive mixing is or is not significant. The degree to which the bulk phase concentration deviates from 
$c_s$ depends on the upstream flow parameter, 
$\varGamma$, the permeability ratio, 
$K$, and the non-dimensional fissure width, 
$\xi ^*$ – see figure 6. Although there is an additional dependence on the vertical extent, 
$l^*$, of the fissure, this dependence is weak and so is not considered in the figure. Concentrations in the regime diagram of figure 6 are measured at 
$x^*=2$, 
$z^*=0$ and at a time 
$t^*$ when the dispersed nose position 
$x^*_{N_d}=10$. Based on figure 6, we surmise that the theoretical model predicts the bulk interface with good accuracy for arbitrary 
$\varGamma$ and 
$K$. On the other hand, our model does a comparatively poor job of predicting the location of the dispersed interface when 
$c^*(2,0) \equiv c(2,0)/c_s < 0.8$. Among other challenges when 
$c^*(2,0) < 0.8$ is the fact that the bulk phase terminates at the location of the fissure, i.e. any gravity current fluid appearing downstream of 
$x^* = 1$ has a density non-trivially less than 
$\rho _s$. This limitation notwithstanding, there remains a large parametric domain over which our theoretical model works well. Figure 6 suggests that as the upstream flow parameter 
$\varGamma$ decreases, so too does the front speed. Less dispersive mixing is therefore observed and the solute concentration after the fissure decreases relatively slowly. As a result, there is a broader range of 
$K$ over which our theoretical model generates predictions in reasonable agreement with the output of the numerical model. Of course, the degree of agreement between theory and numerics is related to the numerical value of the entrainment coefficient, 
$\varepsilon$, which appears e.g. in ((2.16)–(2.18)). We discuss the procedure for determining 
$\varepsilon$ in the following subsection.
Figure 6. Bulk phase concentration reduction beyond the fissure for 
$\theta =0^\circ$ and 
$l^*=0.79$. The red area shows 
$c^*>0.9$, the green area shows 
$0.8< c^*<0.9$ and the blue area shows 
$c^*<0.8$. Boundaries are drawn based on an interpolation performed over a total of 14 simulations for each of the lower and upper surfaces. The inset images show a comparison between theory and numerical simulations for different combinations of 
$\varGamma$ and 
$K$. The thick (thin) white line is the bulk (dispersed) interface as predicted by the theoretical model of § 2. Meanwhile coloured contours show the output of the COMSOL numerical model. Red dashed lines indicate the location 
$x^*=2$, where concentrations are evaluated in constructing the regime diagram.
3.5. Determining the entertainment coefficient
Comparisons such as that depicted in figure 6 (and also figures 8 and 9 below) require that a value be specified for the entrainment coefficient, 
$\varepsilon$. This coefficient is found based on our numerical results. More specifically, we specify 
$\varepsilon$ such that the separation distance, 
$x^*_{N_d} - x^*_{N_b}$, between the dispersed and bulk nose positions matches as closely as possible the distance measured numerically. A mean temporal error is therefore defined as
\begin{equation} \bar{E} = \int_{t^*_1}^{t^*_2} \frac{| (x^*_{N_d}-x^*_{N_b})_{{theory}}- (x^*_{N_d}-x^*_{N_b})_{{num}}| }{(x^*_{N_d}-x^*_{N_b})_{{num}}}\,\mathrm{d}t^*, \end{equation}
where 
$(x^*_{N_d}-x^*_{N_b})_{num}$ is evaluated from the numerical model of § 3 and 
$(x^*_{N_d}-x^*_{N_b})_{theory}$ is evaluated from the theoretical model of § 2 for various 
$\varepsilon$. The (unique) 
$\varepsilon$ that minimizes 
$\bar {E}$ is referred to as the optimum entrainment coefficient. For simplicity, we assume that this optimum value does not depend on 
$\varGamma$ and 
$K$; a justification for this assumption is given in the next paragraph. However, and motivated by the work of Ellison & Turner (Reference Ellison and Turner1959), we allow the error-minimizing value of 
$\varepsilon$ to vary with the inclination angle, 
$\theta$, of the bottom boundary. Results associated with (3.10) and the minimization of 
$\bar {E}$ are displayed in figure 7. They suggest that the optimum value of 
$\varepsilon$ experiences a non-trivial decrease as the slope angle is increased and the gravity current propagates more rapidly downdip.
Figure 7. Error-minimizing value of 
$\varepsilon$ vs 
$\theta$. Here, we consider 
$\varGamma =45$, 
$K=0.3$, 
$l^*=0.79$ and 
$\xi ^*=0.04$. Blue circles consider 
$\varGamma =45$, 
$K=0.2$ and red crosses consider 
$\varGamma =30$, 
$K=0.3$ Also, and with reference to (3.10), 
$t_1^* = 20$ and 
$t_2^* = 70.$
Also included in figure 7 are data corresponding to different values of 
$\varGamma$ and 
$K$. The blue circles indicate the same value of 
$\varGamma$ but a different value of 
$K$. The red crosses indicate the same value of 
$K$ but a different value of 
$\varGamma$. In all cases, we see but a minor deviation from the quantitative data indicated by the black line. In principle, larger changes of 
$K$ can be imagined, however, these would be inconsistent with figure 6, i.e. we limit ourselves to changes of 
$K$ or 
$\varGamma$ that keep us strictly within the red or green sections of the regime diagram.
4. Results and discussion
Within the region of model validity defined in subsection 3.4, theoretical results are compared in figure 8 against COMSOL numerical output. Also included in figure 8 is a sharp interface solution that is obtained by setting 
$\varepsilon = 0$ in (2.16) and (2.17). The sharp interface model over-predicts the nose position while under-predicting the height of the gravity current, especially when the bottom boundary is inclined. A comparison of inclined vs. horizontal gravity currents reveals that the height of the dispersed interface has a monotonic variation with 
$x^*$ when 
$\theta =0^\circ$ but a non-monotonic variation with 
$x^*$ when 
$\theta >0^\circ$. The non-monotonic variation in question becomes more pronounced as time increases.
Figure 8. Gravity current profiles as predicted theoretically and numerically. Line types are as follows: thick solid line – bulk interface; thin solid line – dispersed interface; dashed line – sharp interface solution obtained by setting 
$\varepsilon = 0$ in (2.16) and (2.17). Numerical output is indicated by the colour contours. Panels (a–d) show 
$\theta =0^\circ$, 
$\varGamma =35$, 
$K=0.5$, 
$\xi ^*=0.04$ and 
$l^*=0.79$. Panels (e–h) show 
$\theta =5^\circ$, 
$\varGamma =70$, 
$K=0.3$, 
$\xi ^*=0.04$ and 
$l^*=1.11$. The variation of parameter values between the left- and right-hand side panels is deliberate and illustrates model predictions over a broad range of the parameter space. Note that the scale of the horizontal axis in the left- and right-hand side images is different.
To quantify dispersion effects in gravity currents, we examine the time variation of 
$\tilde {A}^*_{disp}$ from (2.36) and 
$\tilde {B}^*_{disp}$ from (2.38) – see figure 9(c–e). Each panel includes both theoretical and numerical data and considers a different value for 
$K$. Whether 
$K$ is comparatively large (0.4) or small (0.2), the same general trend appears: theory underpredicts the area and buoyancy fraction at early times, however, for sufficiently large 
$t^*$, good agreement is achieved. The early time discrepancy likely appears because the theoretical model assumes a long and thin gravity current. Although the gravity current evolves into a shape that is long and thin for large 
$t^*$, such a large aspect ratio does not apply initially. Furthermore, and for mathematical convenience, our theory is predicated on the assumption that 
$h_2(x=0) = h_1(x=0)$, however, a careful inspection of our numerical results (not shown) indicates some non-trivial initial height of the dispersed phase, at least for not small times. Fortunately, the consequence of ignoring dispersed fluid in the neighbourhood of the source becomes smaller as time progresses and the cumulative volume occupied by the (elongating) gravity current grows.
Figure 9. (a) Area fraction (2.36) and (b) buoyancy fraction (2.38) as a function of the upstream flow parameter (2.26a) for three different values of the fissure parameter (2.26b). Here, we consider a horizontal bottom boundary such that 
$\theta = 0^\circ$ and 
$t^* = 40$. Crosses indicate the solutions for 
$\varGamma = 45$, for which corresponding time series data are given in panels (c–e) for 
$K = 0.2$, 0.3 and 0.4, respectively. These same three 
$K$ values are considered in the time series of panels ( f–h), which consider, again for 
$\varGamma = 45$ and 
$\theta = 0^\circ$, the difference of nose position between the bulk and the dispersed gravity currents. This nose position difference is shown as a function of 
$\varGamma$ in panel (i), where we again consider 
$t^* = 40$.
Whereas figure 9(c–e) assumes the same value for 
$\varGamma$, i.e. 
$\varGamma = 45$, the influence of the upstream flow parameter is explored in figure 9(a,b), where 
$\varGamma$ appears as the 
$x$-axis variable. The different curves of figure 9(a,b) correspond to 
$K = 0.2, 0.3$ and 0.4 such that, with 
$\varGamma = 45$, the red and green regions of figure 6 are spanned appropriately. Curves are drawn at 
$t^* = 40$ by which time the gravity current is indeed long and thin, i.e. the comparison is not unduly influenced by effects related to the early time evolution of the flow. The area and buoyancy fractions increase with 
$\varGamma$ and also with 
$K$. When the upstream flow parameter is large, there is more dispersive mixing because of larger velocities in the gravity current so the fraction of buoyancy and area associated with the dispersed phase increases. By increasing 
$K$, draining of bulk fluid is more robust causing the gravity current to elongate more slowly. Although this has a secondary effect on the dispersed phase (which is, after all, fed by the bulk phase), the overall impact of increasing 
$K$ is to likewise increase the area and buoyancy fractions occupied by the dispersed phase in comparison with the bulk phase. Correspondingly, one might expect that, as 
$K$ increases, so too does the separation distance between the bulk and dispersed nose positions. Figure 9( f–h) confirms this hypothesis and indicate that the theory matches well with the analogous numerical result. Whereas figure 9( f–h) assumes the same value for 
$\varGamma$, i.e. 
$\varGamma = 45$, the influence of the upstream flow parameter is explored in figure 9(i), where 
$\varGamma$ again appears along the abscissa. The different curves of figure 9(i) correspond to 
$K = 0.2, 0.3$ and 0.4 and are drawn at 
$t^* = 40$. Figure 9(i) indicates that 
$x^*_{N_d} - x^*_{N_b}$ increases with both 
$\varGamma$ and 
$K$.
Analogous results have been generated for 
$\theta =5^\circ$ to quantify dispersion effects for the case of inclined gravity currents – see figure 10. Comparing this figure against figure 9 shows that, as expected, dispersion is more robust when 
$\theta >0^\circ$. For instance, and when 
$t^*=60$, figure 10(g) suggests the nose separation for 
$\theta =5^\circ$ is almost twice that in figure 9(g) for 
$\theta =0^\circ$. This is because of the larger characteristic velocity in the inclined case, which leads to more entrainment to the dispersed phase either from the ambient or the bulk phase. This point notwithstanding, very similar trends are observed in figures 9 and 10, e.g. in both cases the separation between the bulk and dispersed noses increases with 
$K$ (due to a stronger drainage) and also with 
$\varGamma$ (due to a larger characteristic velocity). Moreover, theory and numerical simulation demonstrate satisfactory agreement with generally better overlap observed for larger 
$t^*$.
Figure 10. As in figure 9 but with an inclined boundary (
$\theta = 5^\circ$).
Further to the comparison between figures 9 and 10, the sensitivity of our model predictions to the slope of the bottom boundary is more thoroughly explored in figure 11. Figure 11(a) shows, for 
$t^* = 45$, the difference of nose positions for the bulk vs. dispersed phases. Meanwhile figure 11(b) shows the corresponding area and buoyancy fractions for dispersed phase fluid, i.e. 
$\tilde {A}^*_{disp}$ and 
$\tilde {B}^*_{disp}$. Both panels of figure 11 include theoretical and numerical data and confirm the hypothesis that dispersion increases with 
$\theta$.
Figure 11. Difference of (a) nose positions and of (b) dispersed phase area and buoyancy fractions vs. 
$\theta$. Here, 
$t^* = 45$, 
$\varGamma =35$, 
$K=0.5$, 
$l^*=0.79$ and 
$\xi ^*=0.04$.
Because fissure dimensions directly impact drainage, 
$\xi ^*$ and 
$l^*$ also influence the degree of dispersion. Figure 12 confirms that increasing the fissure width, 
$\xi ^*$, leads to more dispersion, whether measured in terms of 
$x^*_{N_d}-x^*_{N_b}$ or 
$\tilde {A}^*_{disp}$. On the other hand, bulk fluid drainage decreases by increasing the vertical extent, 
$l^*$, of the fissure. As a result, there is less dispersion as 
$l^*$ is increased – see figure 13. Consistent with our previous results, figures 12 and 13 indicate that inclined gravity currents experience more dispersion than do horizontal gravity currents.
Figure 12. Difference of (a) nose position and (b) area fraction in the dispersed phase for 
$\theta =0^\circ$ with 
$\varepsilon =0.0125$ and 
$\theta =5^\circ$ with 
$\varepsilon =0.0086$ for various 
$\xi ^*$. Here, 
$t^* = 45$, 
$\varGamma =35$, 
$K=0.5$ and 
$l^*=0.79$.
Figure 13. As in figure 12 but considering the influence of 
$l^*$ for 
$\xi ^* = 0.04$.
5. Application to UHS
To illustrate the application of our results, we return to the example of UHS considered in § 1. More specifically, and for the idealized case of an unbounded reservoir, we wish to estimate the fraction of 
$\mathrm {H}_2$ that will be lost to dispersion as a function of, say, the source volume flow rate. To this end, we consider a line, rather than a point, source such that 
$q_s$ is expressed in units of standard m
$^3$ m
$^{-1}$ day
$^{-1}$ or (S m
$^2$) day
$^{-1}$. Motivated by the numerical investigation of Feldmann et al. (Reference Feldmann, Hagemann, Ganzer and Panfilov2016), we further suppose that 
$\mathrm {H}_2$ is injected into a sandstone layer, bounded above and below by clay layers, where the cushion gas consists of 
$80$ mol
$\%$ 
$\mathrm {N}_2$ and 
$20$ mol
$\%$ 
${\rm CH}_4$. The reservoir pressure is 170 bar such that the density contrast between the 
$\mathrm {H}_2$ and the cushion gas is approximated as 118 kg m
$^{-3}$. The sandstone layer has a porosity of 
$\phi = 13.08$ % and a permeability of 
$k = 22.4$ mD. Meanwhile the clay layer through which the injected 
$\mathrm {H}_2$ drains is idealized as being impermeable except for an isolated fissure situated at a horizontal distance 
$x_f = 50$ m from the source. The fissure is assumed to have a width and depth of 2 m and 1 m, respectively, and is characterized by 
$K = 0.3$, 0.5 or 0.7. Finally, we assume 
$\theta = 0^\circ$ and consider the evolution of the flow over a 10 year period.
Given all of the above parameters, figure 14 shows 
$\tilde {B}^*_{disp}$ as a function of 
$q_s$. As expected from the model predictions of the previous section, the proportion of 
$\mathrm {H}_2$ that mixes with the cushion gas through dispersion decreases with the source volume flow rate. Moreover, and as expected, 
$\tilde {B}^*_{disp}$ is larger when more 
$\mathrm {H}_2$ is allowed to drain, i.e. when 
$K$ is comparatively large. Results such as those shown in figure 14 are helpful because they can, for given 
$K$, identify the minimum source volume flow rate necessary to limit losses by dispersion to a particular value. For example, if, as suggested by the dashed line of figure 14, the maximum loss fraction were set to 5 %, the minimum possible 
$q_s$ could be identified for different fissure permeabilities. Note that this minimum value of the source volume flow rate, 
$(q_s)_{min}$ decreases with 
$K$. Obviously, as 
$K$ tends to zero (indicating a fissure of very limited outflow capacity), 
$(q_s)_{min}$ also tends to zero.
Figure 14. Fraction of 
$\mathrm {H}_2$ lost to dispersion vs. source volume flow rate for the example of § 5. The horizontal dashed line assumes a maximum loss fraction of 5 %.
The preceding analysis can be criticized for prioritizing injectate losses due to dispersion over those due to drainage. Indeed, 
$\mathrm {H}_2$ losses by either mechanism have the potential to make otherwise profitable ventures unattractive economically. On the other hand, there are scenarios such as the ‘selective technology’ advocated by Feldmann et al. (Reference Feldmann, Hagemann, Ganzer and Panfilov2016) for which 
$\mathrm {H}_2$ injection or withdrawal occur simultaneously to/from adjacent sandstone layers. In such a scenario, 
$\mathrm {H}_2$ drained from one layer can be extracted from another layer and so is not necessarily lost to the geological formation. Rather different considerations apply to dispersion because any attempt to produce 
$\mathrm {H}_2$ that has mixed with cushion gas requires the ability, at the surface, to separate 
$\mathrm {H}_2$ from, say, 
${\rm CH}_4$ or 
$\mathrm {N}_2$. The expenses associated with such surface separation operations justify our emphasis on dispersion vs. drainage as a key loss mechanism for 
$\mathrm {H}_2$.
Notwithstanding the conclusion of the previous paragraphs, it should be recalled that our analysis neglects a concentration dependence on viscosity. Strictly speaking, this assumption is incorrect for UHS-type applications; the viscosity of pure 
$\mathrm {H}_2$ is less than that of a mixture of 
$\mathrm {H}_2$ and cushion gas. Because we do not account for the greater mobility of the bulk vs. the dispersed phase, our model likely overestimates the volume of the latter relative to the former, although by how much is not straightforward to quantify. In a similar, though more complicated, spirit our model obviously falls well short of permitting the kinds of fingering instabilities that may arise as a result of a Taylor–Saffman-type instability and the injection of a less viscosity fluid into a more viscous fluid. The modification of our momentum equations to include a concentration-dependent viscosity shall be the subject of future investigations.
6. Summary and conclusions
A theoretical model is developed for a porous medium gravity current consisting of a bulk phase and a dispersed phase – see figure 1. Our theoretical model of § 2 considers local drainage along the bottom boundary, which may be either horizontal or inclined. Model equations are robust enough to capture the essential physical processes of draining and dispersion but are simple enough to be solved using a straightforward numerical algorithm. To this end, we solve the non-dimensional governing equations by defining five non-dimensional parameters namely the inclination angle, 
$\theta$, the upstream flow parameter, 
$\varGamma$, of (2.26a), the permeability ratio, 
$K$, of (2.26b), the fissure width, 
$\xi ^*$, and the fissure length, 
$l^*$. We surmise that all five non-dimensional parameters influence the degree of dispersion. However, 
$l^*$ exerts a subordinate influence compared with 
$\varGamma$, 
$K$, 
$\xi ^*$ and 
$\theta$. Increasing one or both of 
$\varGamma$ and 
$\theta$ increases the gravity current front speed and so increases the degree of dispersion. With reference to the definition of the entrainment velocities 
$w_{e1}$ and 
$w_{e2}$, increasing the gravity current speed makes the entrainment more robust. This supports the idea that the volume of the dispersed phase is significantly larger when we increase parameters such as 
$\varGamma$ or 
$\theta$ that increase the driving force for gravity current flow. Dispersion may also be augmented by causing more (bulk phase) fluid to drain through the fissure, which is realized as either of 
$K$ or 
$\xi ^*$ is increased or 
$l^*$ is decreased. Because drainage directly removes mass from the bulk phase, increasing the drainage leads to more separation between the leading edges (or nose positions) of the bulk and dispersed phases.
Complementing our theoretical results, a COMSOL-based numerical model is developed – see § 3. The numerical model is exploited to estimate the approximate optimum value of the entrainment coefficient, 
$\varepsilon$, which appears as a parameter in the theoretical model e.g. ((2.16)–(2.18)). Through this analysis, we find that the error-minimizing value of 
$\varepsilon$ is a function of the inclination angle, 
$\theta$, as depicted graphically in figure 7. With the appropriate value of 
$\varepsilon$ so selected, we find from figures such as 8–10 generally good agreement between theory and numerical simulation. In other words, our model of § 2 does a reasonable job of predicting the fractions of fluid or solute that appear in the dispersed phase. Our theoretical model also provides generally accurate estimates of the separation distance between the noses of the bulk and dispersed phases. Note, however, that comparisons are restricted to the region of the parameter space for which the degree of draining and subsequent dispersion is not too severe. The results of figure 6 suggest that our theoretical model makes inaccurate predictions of the shape of the dispersed phase when 
$K$ is so large and draining is so vigorous that little or no bulk fluid appears downstream of the source. The modelling of this more complicated case is left for future investigations.
Our study is motivated by the need to address uncertainties in 
$\mathrm {H}_2$ storage in depleted natural gas reservoirs. Mixing of 
$\mathrm {H}_2$ with resident cushion gas is an inevitable facet of depleted reservoir-based UHS systems, particularly in the medium to long term (Feldmann et al. Reference Feldmann, Hagemann, Ganzer and Panfilov2016). Granted our theoretical and numerical models are predicated on a number of simplifying assumptions e.g. ignoring viscosity variations, compressibility effects, ambient counterflow or possible bio-geochemical reactions involving 
$\mathrm {H}_2$. Progressively relaxing these (and related) assumptions are topics for future study. It would also be interesting to consider not localized but rather distributed drainage cf. Pritchard et al. (Reference Pritchard, Woods and Hogg2001), Goda & Sato (Reference Goda and Sato2011) and Bharath et al. (Reference Bharath, Sahu and Flynn2020). Work on this latter problem is already underway and will be reported upon in a future publication.
Funding
This research is supported by NSERC (Discovery Grant program).
Declaration of interests
The authors report no conflict of interest.
Appendix A. Derivation of the bulk and dispersed velocity in the theoretical model
By using the definition of 
$\Delta \rho _1$, 
$\Delta \bar {\rho }_2$ and 
$b_2$ in terms of concentration, (2.7) and (2.8) can be rewritten as
$$\begin{gather} p_1(x,z,t)=\rho_0 g \beta ( b_2 + c_s h_1 )\cos \theta -\rho_{s} g z \cos \theta+\rho_0 g x \sin \theta + P_0 \quad 0 \leq z \leq h_1 , \end{gather}$$
$$\begin{gather}p_2(x,z,t)=[\rho_0 g \beta \bar{c}_2 h_2 - \bar{\rho}_2 g z]\cos \theta +\rho_0 g x \sin \theta +P_0 \quad\quad h_1 \leq z \leq h_2. \end{gather}$$Moreover, Darcy's law (2.12) in the x-direction indicates that
$$\begin{gather} u_1(x,t)={-}\frac{k}{\mu}\left(\frac{\partial p_1}{\partial x}- \rho_s g \sin \theta\right), \end{gather}$$
$$\begin{gather}u_2(x,t)={-}\frac{k}{\mu}\left(\frac{\partial p_2}{\partial x}- \bar{\rho}_2 g \sin \theta\right) . \end{gather}$$
If we insert 
$p_1$ and 
$p_2$ from (A1) and (A2) into (A3) and (A4), the bulk and dispersed velocities then read
$$\begin{gather} u_1(x,t)={-}\frac{kg\beta}{\nu}\left[\frac{\partial b_2}{\partial x}\cos \theta +c_s\left(\frac{\partial h_1}{\partial x}\cos \theta-\sin \theta\right)\right] , \end{gather}$$
$$\begin{gather}u_2(x,t)={-}\frac{kg\beta}{\nu}\left[\left(\frac{\partial (\bar{c}_2 h_2)}{\partial x}-z\frac{\partial \bar{c}_2 }{\partial x}\right)\cos \theta-\bar{c}_2\sin \theta\right]. \end{gather}$$
COMSOL results show that, for the conditions relevant to our analysis, the term 
$z({\partial \bar {c}_2 }/{\partial x})$ in (A6) is two orders of magnitude smaller than 
${\partial (\bar {c}_2 h_2)}/{\partial x}$ and can therefore be ignored. Accordingly, velocity in the dispersed phase can be simplified to
\begin{equation} u_2(x,t)={-}\frac{kg\beta}{\nu}\left[\frac{\partial (\bar{c}_2 h_2)}{\partial x}\cos \theta-\bar{c}_2\sin \theta\right].\end{equation}Appendix B. Derivation of the drainage velocity in the theoretical model
Using (2.7), the pressure at the bottom boundary of the gravity current is expressed as
\begin{equation} p(x,0,t)=[\Delta \bar{\rho}_2 g h_2+ ( \Delta \rho_1-\Delta \bar{\rho}_2 )g h_1 ]\cos \theta +\rho_0 g x \sin \theta + P_0 . \end{equation}
Moreover, and assuming a hydrostatic pressure balance, the pressure at 
$z=-l$ corresponding to the base of the fissure is given by
\begin{equation} p(x,-l,t)=\rho_0gl \cos \theta +\rho_0 g x \sin \theta + P_0 . \end{equation}
Analogously to Acton, Huppert & Worster (Reference Acton, Huppert and Worster2001), and by considering pressure continuity at 
$z=0$, the pressure distribution within the fissure is described by the following linear function:
\begin{align} p(x,z,t)&=[\Delta \bar{\rho}_2 g h_2+ ( \Delta \rho_1-\Delta \bar{\rho}_2 )g h_1 ]\left(1+\frac{z}{l}\right)\cos \theta\nonumber\\ &\quad -\rho_0 g z \cos \theta+\rho_0 g x \sin \theta + P_0 \quad\quad -l \leq z \leq 0. \end{align}Applying Darcy's law, the vertical velocity in the fissure reads
\begin{equation} w_{d}(x,t)={-}\frac{k_f}{\mu}\left(\frac{\partial p}{\partial z}+\rho_0\, g\, \cos \theta\right)={-}\frac{k_f g}{\mu} \left[\frac{\Delta \bar{\rho}_2 h_2}{l}+( \Delta \rho_1-\Delta \bar{\rho}_2 )\frac{ h_1}{l}+\Delta \rho_1 \right]\cos \theta.\end{equation}
If we insert 
$\Delta \rho _1=\rho _0\beta c_s$ and 
$\Delta \bar {\rho }_2=\rho _0\beta \bar {c}_2$ into (B4), it can be shown that
\begin{equation} w_{d}(x,t)={-}\frac{k_f g\beta}{\nu} \left(\frac{c_sh_1+b_2}{l}+c_s\right) \cos \theta . \end{equation}Appendix C. Method of solution for the theoretical model
The finite difference method is used to discretize equations (2.27)–(2.29) in space. First-order derivatives in space are discretized using backward finite differences and a central finite difference is used to discretize second-order derivatives. Although implicit methods are more stable, they produce extra diffusion in our solution; we therefore apply an explicit method for time discretization. Thus (2.27)–(2.29) may be rewritten in discrete form as
\begin{align} &(h^*_{1,i})^{n+1}=
(h^*_{1,i})^{n}-\Delta t^* (h^*_{1,i} \, {\rm d}U^*_i)^n
-\Delta t^* (U^*_i)^n \left( \frac{h^*_{1, i} -h^*_{1,
i-1}}{\Delta x^*} + \varepsilon \varGamma^{1/2}\right)^n
\nonumber\\ &\quad -\Delta t^* K \varGamma
\left(\frac{h^*_{1, i}+b^*_{2, i}}{l^*}+1\right)^n \cos
\theta \times F(x^*,1,\xi^*) ,\end{align}
\begin{align}
&(h^*_{2,
i})^{n+1} =(h^*_{2, i})^n+\frac{\Delta t^*}{(\Delta
x^*)^2}[(h^*_{2, i}-h^*_{1, i})(\varPsi^*_{i-1}-2
\varPsi^*_i +\varPsi^*_{i+1}- (C^*_i-C^*_{i-1})\Delta
x^*)\nonumber\\ &\quad +(h^*_{2, i}-h^*_{2, i-1}-h^*_{1,
i}+h^*_{1, i-1}) (\varPsi^*_{i}-\varPsi^*_{i-1}-C^*_i
\Delta x^*)]^n \nonumber\\ &\quad -\Delta t^* \left(h^*_{1, i} \, {\rm
d}U^*_i +U^*_i \frac{h^*_{1, i}-h^*_{1, i-1}}{\Delta
x^*}\right)^n - \varepsilon
\varGamma^{1/2} \Delta t^* \left(
\frac{\varPsi^*_i-\varPsi^*_{i-1}}{\Delta x^*} -C^*_i
\right)^n\nonumber\\ &\quad -\Delta t^* K \varGamma \left(\frac{h^*_{1,
i}+b^*_{2, i}}{l^*}+1\right)^n \cos \theta \times
F(x^*,1,\xi^*) ,\end{align}
\begin{align}
&(b^*_{2,i})^{n+1}=(b^*_{2, i})^n+\frac{\Delta
t^*}{(\Delta x^*)^2}[b^*_{2, i}
(\varPsi^*_{i-1}-2\varPsi^*_i +
\varPsi^*_{i+1}-(C^*_i-C^*_{i-1})\Delta x^*) \nonumber\\
&\quad +
(b^*_{2,i}-b^*_{2,i-1})(\varPsi^*_{i}-\varPsi^*_{i-1}-C^*_i
\Delta x^*) ]^n +\Delta t^* \varepsilon \varGamma^{1/2}(U^*_i)^n ,
\end{align}
respectively. Here, 
$i$ and 
$n$ are non-negative indices that respectively correspond to space and time. In addition, 
$(U^*_i)^n$ and 
$(dU^*_i)^n$ are defined as
\begin{align}
&(U^*_i)^n={-}\left(\frac{b^*_{2, i}-b^*_{2, i-1}}{\Delta
x^*}+ \frac{h^*_{1, i}-h^*_{1, i-1}}{\Delta x^*} \right)^n
\cos \theta +\varGamma^{1/2} \sin \theta ,
\end{align}
\begin{align}
&(dU^*_i)^n={-}\left(\frac{\partial^2 b^*_2}{\partial
x^*\,^2} + \frac{\partial^2~h^*_1}{\partial
x^*\,^2}\right)_i^n \cos \theta\nonumber\\
&\quad
={-}\left(\frac{b^*_{2, i-1}-2b^*_{2, i}+b^*_{2,
i+1}}{(\Delta x^*)^2}+\frac{h^*_{1, i-1}-2 h^*_{1,
i}+h^*_{1, i+1}}{(\Delta x^*)^2} \right)^n \cos \theta .
\end{align}
Equation (C1) applies for 
$i\geq 1$. When 
$i=1$, 
$(h^*_{1, 0})^{n}$ in (C1) is found based on the discrete form of the influx boundary condition in (2.34a), such that
\begin{equation} (h^*_{1, 0})^{n}=\frac{\Delta x^*}{\phi (h^*_{1, 1})^{n} \cos \theta}+(h^*_{1, 1})^n-(b^*_{2, 0})^n- \varGamma^{1/2} \tan \theta \Delta x^* .\end{equation}
To find 
$(h^*_{1, 0})^{n}$ in (C6), we assume some amount for 
$(b^*_{2, 0})^n$ and solve (C2) and (C3) for 
$i=1$ to recover 
$(h^*_{2, 1})^{n+1}$ and 
$(b^*_{2, 1})^{n+1}$. We then iterate using the secant method to satisfy boundary conditions (2.34c) and (2.34e), i.e.
$$\begin{gather} (h^*_{2, 1})^{n+1}=(h^*_{1, 1})^{n+1} , \end{gather}$$
$$\begin{gather}(b^*_{2, 1})^{n+1}=0 . \end{gather}$$
Then the expressions in (C2) and (C3) apply for 
$i\geq 2$. Finally, for 
$i=N_b$ and 
$i=N_d$, the bulk and dispersed nose positions in (2.34c–f) read
$$\begin{gather} (h^*_{1, N_b})^{n+1} = 0 , \end{gather}$$
$$\begin{gather}(h^*_{2, N_d})^{n+1} = (b^*_{2, N_d})^{n+1}=0 . \end{gather}$$
In the above equations, 
$\Delta x^*$ and 
$\Delta t^*$ indicate the grid spacing and the time step, respectively. Our discretized equations are solved with 
$\Delta x^*= 10^{-2}$ and 
$\Delta t^*=10^{-5}$. We estimate 
$\Delta x^*$ by fixing 
$\Delta t^*$ and then performing a grid-independence test. Once the largest value of 
$\Delta x^*$ that preserves grid independency is determined, 
$\Delta t^*$ is increased slightly, but not beyond a value where computed results vary with the magnitude of the time step. With suitable values for 
$\Delta x^*$ and 
$\Delta t^*$ selected, we find that the run time to produce a figure such as figure 2 is approximately 0.65 core hours using an Intel Core i7-9700 CPU (3.00 GHz and 16 GB memory).
Open access
