2.1. Problem description

To theoretically model the flow described in figure 1, we consider a two-layer porous medium with notation as illustrated schematically in figure 2. The permeability and porosity in the upper layer are $k_1$ and $\phi _1$, respectively, while the corresponding values in the lower layer are $k_2$ and $\phi _2$, respectively. The permeability jump makes a constant angle $\theta$ with the horizontal. The natural coordinate system in two-dimensional space is represented by (${X,Z}$), where the vertical Z-axis is aligned anti-parallel to the gravitational acceleration $g$. For reference, note that the coordinate system (${x,z}$) associated with the along- and cross-jump directions is obtained from (${X,Z}$) by a clockwise rotation of $\theta$ about the Y -axis, $Y$ being the direction normal to the page. Coordinate rotation can be expressed mathematically using a transformation matrix, i.e.

(2.1)\begin{equation} \begin{bmatrix}{x} \\ {z}\end{bmatrix} = \begin{bmatrix} \mathrm{cos}\ \theta\ & \sin \theta \\ -\sin \theta\ & \cos \theta \end{bmatrix} \begin{bmatrix} {X} \\ {Z} \end{bmatrix}. \end{equation}

We position the origin of both coordinate systems on the permeability jump, directly below the dense line source, which is itself is positioned at ($X=0, Z=H$), see figure 2.

Figure 2.

Definition schematic showing the propagation of discharged plume fluid in the up- and downdip directions along the inclined permeability jump. Draining into the lower layer is also indicated.

Our theoretical model is predicated on the following simplifying assumptions. (i) Although the upper and lower layers are assumed to support different permeabilities, the porosities are assumed equal, i.e. $\phi _1 = \phi _2 \equiv \phi$. (ii) The upper and lower layers are assumed to be very deep in vertical extent (Goda & Sato Reference Goda and Sato2011); the consequences of assuming otherwise will be briefly highlighted in § 4. (iii) Initially, the entire porous medium is assumed to be uniformly saturated with ambient fluid of density $\rho _{o}$. (iv) The source and ambient fluids have equal dynamic viscosities and are assumed to be fully miscible so that effects due to capillarity can be ignored. (v) When the plume strikes the permeability jump at $x=0$, all of its fluid is then discharged in the form of flows propagating up- and down-dip; this assumption is defensible provided $k_2/k_1 \ll 1$, see for example figure 3 of Sahu & Flynn (Reference Sahu and Flynn2017). (vi) The so produced up- and downdip gravity currents remain long and thin. (vii) Consistent with Bear (Reference Bear1972), Woods & Mason (Reference Woods and Mason2000), Pritchard et al. (Reference Pritchard, Woods and Hogg2001) and De Loubens & Ramakrishnan (Reference De Loubens and Ramakrishnan2011), gravity current fluid is at all times separated from overlying ambient fluid by a sharp interface. (viii) The density difference between the source fluid and the ambient fluid is moderate and the flow remains Boussinesq everywhere in the domain. (ix) Spatial variations in the density within the up- and downdip gravity currents and within the contaminated fluid consisting of discharged plume fluid that has drained into lower layer are modest and can be ignored to leading order.

The plume line source discharges fluid of density $\rho _s>\rho _o$ at a constant volume flux per unit width, $q_{s}$. The source buoyancy flux per unit source width is $F_s=q_{s} {g}^\prime _{s}$, where $g'_s= g (\rho _s-\rho _o)/\rho _o\ll {g}$ is the source reduced gravity. As the dense fluid falls downwards, entrainment occurs as a result of which the volume flux of the resulting plume increases with the vertical distance from the source. Using a boundary layer approximation, Wooding (Reference Wooding1963) derived an expression to predict this variation in the limit of small Péclet number, i.e. ${Pe}=U d_o\tau /D_m \ll O(1)$, where $U$ is a characteristic velocity, $d_o$ is the mean grain diameter, $D_m$ is the molecular diffusivity. Further, $\tau > 1$ is the (hydraulic) tortuosity, which is defined as the square of the hydraulic flow path length to the corresponding straight line length (Carman Reference Carman1937) or as the product of porosity and a formation factor (Clennell Reference Clennell1997; Ghanbarian et al. Reference Ghanbarian, Hunt, Ewing and Sahimi2013). Later, Sahu & Flynn (Reference Sahu and Flynn2015) derived a similar relation for ${Pe}\gg O(1)$ and also proposed a relation for the variation of the plume reduced gravity with $Z$. Herein, we model the flow assuming ${Pe}\gg O(1)$, see § 2.4.

When the plume strikes the permeability jump, discharged plume fluid is divided into equal ($\theta ={0}^{\circ }$) or unequal ($\theta \neq {0}^{\circ }$) flows to the right and left. The discharged plume fluid propagates as a pair of leaky gravity currents under a balance of buoyancy and viscosity with simultaneous draining into the lower layer. The pressure is continuous along $z=h(x,t)$ where $h$ denotes the gravity current height, measured perpendicular to the permeability jump, see figure 2. Because of the continual addition of fluid from above, $h$ steadily increases with time, $t$. Consequently, and because the density, $\rho _p$, and plume volume flux per unit width, $q_p$, vary with the vertical coordinate due to entrainment, the fluid feeding the gravity current has a density that slowly increases with time and a volume flux that slowly decreases with time. The influx density is denoted by $\rho _p(h_{0})$ with a corresponding volume flux per unit width denoted by $q_p(h_{0})$ where $h_0 \equiv h(x=0,t) \equiv h_0(t)$ is the time-dependent gravity current height measured directly below the plume source. The volume influx $q_p(h_0)$ is obviously the same $q_c$. In a similar spirit, and neglecting spatial variations of the gravity current density, we have that $\overline {\rho _c} = \rho _p(h_{0})$ where $\rho _o < \overline {\rho _c} <\rho _s$. Although the spatially uniform approximation is clearly a theoretical simplification relative to the real flow, we expect internal stratification effects to be relatively minor provided $H$ is not small. With this assumption to hand, the average density difference between the up- and downdip gravity currents and the ambient fluid is ${\rm \Delta} \overline {\rho _c} = \overline {\rho _{c}}-\rho _{o}$ and the corresponding reduced gravity is $\overline {g'_c}=g ({\rm \Delta} \overline {\rho _c}/\rho _o)$.

Below the permeability jump, there is a draining of discharged plume fluid into the lower layer. The draining flow displaces less dense ambient fluid thereby creating an unstable interface leading to Rayleigh–Taylor type fingering (Saffman & Taylor Reference Saffman and Taylor1958; Homsy Reference Homsy1987), see figure 1. Occurrence of similar fingering phenomena has been observed previously for gravity currents propagating along a horizontal permeability jump (Sahu & Flynn Reference Sahu and Flynn2017) or along the underside of a sloping, impermeable boundary (MacMinn & Juanes Reference MacMinn and Juanes2013). A consequence of the fingers is that some mixing of draining discharged plume fluid and ambient fluid must occur. To distinguish the less dense discharged plume fluid in the lower layer from that more dense discharged plume fluid in the upper layer, we shall, in the lower layer, make reference to contaminated fluid, which in figure 2 is characterized by a depth $l$. By contrast, when we refer to discharged plume fluid, it should hereafter be understood that we refer specifically to the upper layer.

As a further consequence of the mixing described above, we cannot, strictly speaking, assume a sharp interface in the lower layer, cf. figure 1. Nonetheless, and for theoretical expediency, we assume a spatially uniform mixing process in the lower layer and thereby define an average interface location for the contaminated fluid relative to the ambient, see figure 2. The average density of the contaminated fluid is $\overline {\rho _d}<\overline {\rho _c}$ and the corresponding density difference with the ambient fluid is ${\rm \Delta} \overline {\rho _d} = \overline {\rho _d }- \rho _o$. Thus the mean reduced gravity of the draining fluid is expressed as $\overline {g'_d}=g ({\rm \Delta} \overline {\rho _d}/\rho _o)<\overline {g'_c}$. Determination of $\overline {\rho _d}$ and $\overline {g'_d}$ will be discussed below.

2.2. Gravity currents

Consider, as in figure 2, a two-dimensional flow of the gravity currents propagating along the permeability jump in the upper layer along the up- and downdip directions. Let the Darcy velocity of the gravity currents be ${{\boldsymbol {u}}}_c\equiv ({u}_{c},{w}_{c})$, with components $u_c$ and $w_c$ in the along- and cross-jump directions, respectively. As noted above, the gravity currents are assumed long and thin (small aspect ratio, $\varepsilon = \text {height/length} \ll 1$). As a consequence, and consistent with the Dupuit approximation (see Bear (Reference Bear1972) and appendix A.3) ${w}_{c}$ can be considered small compared with ${u}_{c}$ and the cross-layer pressure gradient within each gravity current can be considered approximately hydrostatic (Huppert & Woods Reference Huppert and Woods1995). Considering the along-jump flows in a rotated coordinate system for $\theta > {0}^{\circ }$, we define the hydrostatic pressure in the cross-flow direction as

(2.2)\begin{equation} {P}_c(x,z,t)={p}_{c}(x,t)-\overline{\rho_c} g z\cos \theta. \end{equation}

Here, $p_c(x,t)$ is given by

(2.3)\begin{equation} {p}_{c}(x,t)= P_o+ (\overline{\rho_c}-{\rm \Delta}\overline{\rho_c}) {g}x\sin \theta+ {\rm \Delta} \overline{\rho_c} g h \cos \theta, \end{equation}

where $P_o$ is the pressure measured at the origin. Using Darcy's law, it can be shown that the along-jump flow velocities with the up- and downdip gravity currents are

(2.4)\begin{equation} u_c= -\frac{k_{1} {\rm \Delta} \overline{\rho_c} g}{\mu} \begin{cases} \displaystyle \frac{\partial{h}}{\partial{x}}\cos \theta + \sin \theta, & [\text{updip}, -x_{N_u} < x < 0] \\ \displaystyle \frac{\partial{h}}{\partial{x}}\cos \theta - \sin \theta, & [\text{downdip}, 0 < x < x_{N_d}] \end{cases} \end{equation}

where $\mu$ is dynamic viscosity and $x_{N_u}$ and $x_{N_d}$ are the nose locations in the up- and downdip directions, respectively, see figure 2. The velocities prescribed by (2.4) apply for $z < h$. Because the upper layer is assumed semi-infinite, the velocity within the ambient is considered negligible (cf. De Loubens & Ramakrishnan Reference De Loubens and Ramakrishnan2011). To derive an evolution equation for $h$, we take the depth average of the mass conservation equation, similar to (2.6) of Huppert & Woods (Reference Huppert and Woods1995), and consider mass loss due to draining from the gravity current undersides. Upon substituting (2.4) into the resulting depth-averaged equation, we obtain

(2.5)\begin{equation} \phi \frac{\partial{h}}{\partial{t}}= \begin{cases} \displaystyle \frac{k_{1} \overline{g'_c}}{\nu} \frac{\partial{}}{\partial{x}} \left(h\frac{\partial{h}}{\partial{x}}\cos \theta +h \sin \theta \right) + w_{drain}, & [\text{updip}, -x_{N_u} < x < 0] \\ \displaystyle \frac{k_{1} \overline{g'_c}}{\nu} \frac{\partial{}}{\partial{x}}\left( h\frac{\partial{h}}{\partial{x}}\cos \theta -h \sin \theta \right) + w_{drain}, & [\text{downdip}, 0 < x < x_{N_d}] \end{cases} \end{equation}

where $\nu =\mu /\rho _o$ denotes the kinematic viscosity and $w_{drain}$ is the draining velocity, an expression for which is provided in the following subsection.

2.3. Draining flow in the lower layer

Here, we evaluate the time evolution of the average position of the interface between the contaminated and the ambient fluids in the lower layer. The Darcy velocity of the contaminated fluid in the lower layer is ${{\boldsymbol {u}}}_d\equiv ({u}_{d},{w}_{d})$, having components $u_d$ and $w_d$ in the along- and cross-jump directions, respectively. For $\theta ={0}^{\circ }$, the along-jump velocity in the lower layer $u_{d}$ is considered by Acton et al. (Reference Acton, Huppert and Worster2001) to depend only on the horizontal gradient of the hydrostatic pressure exerted by the gravity current, i.e. $\partial {h}/\partial {x}$. Meanwhile, the cross-jump component of velocity, $w_{d}$, depends not on $\partial h/\partial x$ but rather on $h$. Because $\partial {h}/\partial {x}$ can be shown to be both approximately constant and small compared with unity, $|{u}_{d}|\ll \lvert {w}_{d} \rvert$. Consequently, the influence of ${u}_{d}$ on the draining velocity has typically been ignored in previous works. However, for inclined permeability jumps ($\theta \neq {0}^{\circ }$), the along-jump velocities may be significant because they include gravitational acceleration projected into the along-jump direction. Using Darcy's law, the along-jump velocities just below the jump boundary are given by

(2.6)\begin{equation} u_d= -\frac{k_{2} {\rm \Delta}\overline{\rho_d }g}{\mu} \begin{cases} \displaystyle \frac{\partial{h}}{\partial{x}}\cos \theta + \sin \theta, & [\text{updip}, -x_{N_u} < x < 0] \\ \displaystyle \frac{\partial{h}}{\partial{x}}\cos \theta - \sin \theta, & [\text{downdip}, 0 < x < x_{N_d}]. \end{cases} \end{equation}

Meanwhile, the cross-jump component is given by

(2.7)\begin{equation} {w}_{d}=-\frac{k_{2} {\rm \Delta}\overline{\rho_d }g}{\mu} \left( 1 + \frac{h}{l'} \right)\cos \theta, \end{equation}

where $l'$ is the depth of the contaminated fluid as measured below and perpendicular to the permeability jump, see figure 2. In contrast to the case of a horizontal permeability jump, $|{u}_{d}| \approx \lvert {w}_{d} \rvert \tan \theta$, so that $|{u}_{d}|$ can be neglected only for relatively modest $\theta$. For larger $\theta$, ${u}_{d}$ becomes significant and needs to be considered when modelling the draining flow. In such cases, the net draining flow velocity $w_{drain}$ is influenced by both the along-jump velocity and the hydrostatic head in the cross-jump direction. To evaluate $w_{drain}$, we take the resultant of these two velocities, i.e. $\sqrt {{u}^2_{d}+{w}^2_{d}}$. Consistent with figure 2, it can be shown that for draining lengths significantly greater than the gravity current height, the direction of $w_{drain}$ remains vertical and its length is given as $l=l'/\cos \theta$. Taking these factors into consideration, the draining velocity can be written as

(2.8)\begin{equation} w_{drain}=-\frac{k_{2} {\rm \Delta}\overline{\rho_d}g}{\mu} \left(1 + \frac{h}{l} \cos \theta \right). \end{equation}

In the limiting case when $\theta =0^{\circ }$, (2.8) reverts back to the expression used in previous works (Goda & Sato Reference Goda and Sato2011; Sahu & Flynn Reference Sahu and Flynn2017).

At the contaminated–ambient fluid interface, the volume flow rate per unit width, $q_{entr}$, of the ambient fluid that mixes into the draining gravity current fluid can be expressed in terms of $w_{drain}$ as

(2.9)\begin{equation} q_{entr}=\frac{{\rm \Delta}\overline{\rho_c} - {\rm \Delta}\overline{\rho_d}}{{\rm \Delta} \overline{\rho_d}}\int_{-x_{N_u}}^{x_{N_d}} w_{drain} \,\mathrm{d}x, \end{equation}

see appendix A.1. The entrainment prescribed by $q_{entr}$ is important because it, along with $w_{drain}$, dictates the time rate of increase of the draining fluid length, $l$. Following the derivation of appendix A.1 and incorporating (2.8), it can ultimately be shown that $l$ satisfies the following evolution equation:

(2.10)\begin{equation} \displaystyle \phi \frac{\partial{l}}{\partial{t}} = - \frac{\overline{g'_c}}{\overline{g'_d}}w_{drain}= \frac{k_{1} \overline{g'_c}}{\nu}\frac{k_{2} }{k_{1}}\left( 1 + \frac{h}{l} \cos \theta \right). \end{equation}

As with (2.5), (2.10) is valid in the range $-x_{N_u}<x<x_{N_d}$.

2.4. Gravity currents and draining flows fed by a descending plume

Recall that the density of the fluid supplied to the gravity currents slowly increases with time as a result of the gradual increase of $h$. Wishing to account for this fact in the governing equations (2.5) and (2.10) we adopt (2.25) of Sahu & Flynn (Reference Sahu and Flynn2015), and write

(2.11)\begin{equation} \displaystyle {g}^\prime_{c}=\left[\left(\frac{{\rm \pi} F_{s} \nu}{16 k_{1}} \right) ^2 \frac{1}{ \phi \alpha (H+ {Z}_{s} - h_0\cos \theta)}\right]^{{1}/{4}}. \end{equation}

Here, the source correction term ${Z}_{s}= ({1}/{\phi \alpha}) ({{\rm \pi} \nu }/{16 F_{s} k_{1}})^2 q_{s}^{4}$ is evaluated considering the volumetric flow rate $q_s$ at the plume source. Also, $\alpha$ is the transverse dispersivity, whose connection to the porous medium grain size is discussed below in § 5.1. Substituting (2.11) into (2.5) and (2.10) yields, after some simplification,

(2.12)\begin{equation} \frac{\partial{h}}{\partial{t}}= \displaystyle \begin{cases} \displaystyle \beta (1- \chi h_{0^-}\cos \theta) ^{-{1}/{4}} \left[\frac{\partial{}}{\partial{x}}\left( h\frac{\partial{h}}{\partial{x}}\cos \theta +h \sin \theta \right) - K G^{\prime} \left(1 + \frac{h}{l}\cos \theta \right) \right],\\ \quad [\text{updip}, -x_{N_u} < x < 0], \\ \displaystyle \beta (1- \chi h_{0^+} \cos \theta) ^{-{1}/{4}} \left[\frac{\partial{}}{\partial{x}}\left(h\frac{\partial{h}}{\partial{x}}\cos \theta -h \sin \theta \right) - K G^{\prime} \left(1 + \frac{h}{l} \cos \theta \right) \right], \\ \quad [\text{downdip}, 0 < x < x_{N_d}], \end{cases} \end{equation}

and

(2.13)\begin{equation} \frac{\partial{l}}{\partial{t}}= \displaystyle \beta K \begin{cases} \displaystyle (1- \chi h_{0^-} \cos \theta)^{-{1}/{4}} \left(1 + \frac{h}{l} \cos \theta \right) , & [\text{updip}, -x_{N_u}< x < 0], \\ \displaystyle (1- \chi h_{0^+} \cos \theta) ^{-{1}/{4}} \left(1 + \frac{h}{l} \cos \theta \right) , & [\text{downdip}, 0 < x < x_{N_d}]. \end{cases} \end{equation}

Here $\beta$ is a velocity parameter, $\chi$ is a source parameter, $K$ is the ratio of the lower to upper layer permeabilities and $G^{\prime }$ is the ratio of the lower to upper layer reduced gravities. More precisely, and in symbols,

(2.14a–d)\begin{align} \displaystyle \beta= \frac{k_{1} }{\phi \nu} \left[\left(\frac{{\rm \pi} F_{s} \nu}{16 k_{1}} \right) ^2 \frac{ 1}{ \phi \alpha (H+Z_{s})}\right]^{{1}/{4}} ,\quad \displaystyle \chi =\frac {1}{H+Z_{s}},\quad \displaystyle K=\frac{k_{2}}{k_{1}},\quad \displaystyle G^{\prime} = \frac{\overline{g'_d}}{\overline{g'_c}}. \end{align}

Regarding the permeability jump angle $\theta$ as a further independent parameter, there are a total of five variables in the governing equations (2.12) and (2.13).

2.5. Initial and boundary conditions

Recalling the initial condition of the porous medium to be uniformly saturated by ambient fluid, we initialize the gravity current height and draining length to zero, i.e. $h=l=0$ at $t=0$. For $t>0$, (2.12) and (2.13) are solved using an influx boundary condition, which requires the time rate of volume increase of the up- and downdip gravity currents to balance the volume of fluid supplied by the descending plume. The plume volume flux per unit width, $q_c$, supplied at $x=0$, increases with distance from the source. Analogous to (2.11), $q_c$ can be expressed in terms of $h_0$ using

(2.15)\begin{equation} {q}_{c}= \left[\left(\frac{16 F_{s}k_{1}}{{\rm \pi} \nu} \right) ^2 \phi \alpha (H + {Z}_{s} - h_{0}\cos \theta) \right]^{{1}/{4}}, \end{equation}

see (2.24) of Sahu & Flynn (Reference Sahu and Flynn2015). The $q_c$ in (2.15) is divided into unequal up- and downdip components for $\theta \neq 0^{\circ }$. To this end, and borrowing the notation of Rayward-Smith & Woods (Reference Rayward-Smith and Woods2011), we consider that the dimensionless fraction of the flow propagating downdip is $f_a$, while the remaining fraction travelling updip is $1-f_a$. The influx boundary conditions are then represented as

(2.16)\begin{align} \left.\begin{array}{@{}cc@{}} \displaystyle \left.\beta ^2 \left( h\frac{\partial{h}}{\partial{x}}\cos \theta +h \sin \theta \right) \right|_{0^-} = - (1-f_a) \varGamma (1-\chi h_{0^-} \cos \theta)^{{1}/{2}}, & [\text{updip}, -x_{N_u} < x < 0], \\ \displaystyle \left.\beta ^2 \left( h\frac{\partial{h}}{\partial{x}}\cos \theta -h \sin \theta \right) \right|_{0^+} = - f_a \varGamma (1-\chi h_{0^+} \cos \theta)^{{1}/{2}}, & [\text{downdip}, 0 < x < x_{N_d}], \end{array}\right\} \end{align}

where $\varGamma = (k_{1} F_{s})/(\phi ^2 \nu )$ is a buoyancy flux factor. The two components of (2.12) are coupled by insisting that the gravity current height remains continuous at $x=0$, i.e.

(2.17)\begin{equation} \displaystyle h_{0^-}= h_{0^+} . \end{equation}

By enforcing this condition at each time step, $f_a$ can be determined as a function of time $t$. A final boundary condition is applied at the noses of the up- and downdip gravity currents, such that

(2.18a,b)\begin{equation} h_{-x_{N_u}}=l_{-x_{N_u}}=0\quad \text{and} \quad h_{x_{N_d}}=l_{x_{N_d}}=0. \end{equation}

Finally, the global mass balance equation can be written symbolically as

(2.19)\begin{align} \int_{-x_{N_u}}^{x_{N_d}} (h+\left|l\right|) \,\mathrm{d}x= \left[\frac{\varGamma}{\beta} (1-\chi h_{0} \cos \theta)^{{1}/{4}} + \beta K (1-G^{\prime}) \int_{-x_{N_u}}^{x_{N_d}} \left(1 + \frac{h}{l} \cos \theta \right) \mathrm{d}x \right]t. \end{align}

The former term on the right-hand side represents the volume of fluid discharged by the plume while the latter term corresponds to the volume of lower layer ambient fluid entrained into the contaminated fluid.

2.6. Dimensionless governing equations

Following Goda & Sato (Reference Goda and Sato2011), we define characteristic spatial and temporal variables, ${\varPi _x}$ and ${\varPi _t}$, as follows:

(2.20a,b)\begin{equation} \displaystyle \varPi_x=\frac{q_c|_{h=0}}{\phi \beta} \quad \text{and}\quad \displaystyle \varPi_t=\frac{q_c|_{h=0}}{\displaystyle (1- \delta \cos \theta)^{-{1}/{4}} \phi \beta^2}, \end{equation}

where

(2.21)\begin{equation} \delta= \frac{16}{\rm \pi} \left(\frac {\phi \alpha}{H+Z_{s}}\right)^{{1}/{2}}. \end{equation}

Note that $\varPi _x$ characterizes the distance measured along the permeability jump, whereas $\varPi _t$ characterizes the speed of draining into the lower layer. Using the above characteristic variables, we non-dimensionalize other variables as follows:

(2.22a–d)\begin{equation} x^\ast=\frac{x}{\varPi_x}, \quad h^\ast=\frac{h}{\varPi_x},\quad l^\ast=\frac{l}{\varPi_x},\quad t^\ast=\frac{t}{\varPi_t} \end{equation}

Thus (2.12) and (2.13) may be rewritten, respectively, as

(2.23)\begin{equation} \frac{\partial{h^\ast}}{\partial{t^\ast}}= \begin{cases} \displaystyle \left(\frac{1- \delta h^\ast_{0^-} \cos \theta}{1- \delta \cos \theta}\right)^{-{1}/{4}} \left[\frac{\partial{}}{\partial{x^\ast}} \left(h^\ast\frac{\partial{h^\ast}}{\partial{x^\ast}}\cos\theta +h^\ast \sin \theta \right)- K G^{\prime} \left( 1 + \frac{h^\ast}{l^\ast} \cos \theta \right) \right],\\ \quad [\text{updip}, -x^\ast_{N_u} < x^\ast < 0], \\ \displaystyle \left(\frac{1- \delta h^\ast_{0^+} \cos \theta}{1- \delta \cos \theta}\right)^{-{1}/{4}} \left[\frac{\partial{}}{\partial{x^\ast}}\left(h^\ast\frac{\partial{h^\ast}}{\partial{x^\ast}}\cos \theta -h^\ast \sin \theta \right) - K G^{\prime}\left( 1 + \frac{h^\ast}{l^\ast} \cos \theta \right) \right],\\ \quad [\text{downdip}, 0 < x^\ast < x^\ast_{N_d}], \end{cases} \end{equation}

and

(2.24)\begin{equation} \frac{\partial{l^\ast}}{\partial{t^\ast}}= K \begin{cases} \displaystyle \left(\frac{1- \delta h^\ast_{0^-} \cos \theta}{1- \delta \cos \theta}\right) ^{-{1}/{4}} \left( 1 + \frac{h^\ast}{l^\ast} \cos \theta \right), & [\text{updip}, -x^\ast_{N_u} < x^\ast < 0], \\ \displaystyle \left(\frac{1- \delta h^\ast_{0^+} \cos \theta}{1- \delta \cos \theta}\right)^{-{1}/{4}} \left( 1 + \frac{h^\ast}{l^\ast} \cos \theta \right), & [\text{downdip}, 0 < x^\ast < x^\ast_{N_d}]. \end{cases} \end{equation}

The initial condition reads $h^\ast =l^\ast =0$. Meanwhile, the boundary conditions (2.16) become

(2.25)\begin{align} \left.\begin{array}{cc@{}} \displaystyle \left.\left( h^\ast\frac{\partial{h^\ast}}{\partial{x^\ast}}\cos \theta +h^\ast \sin \theta \right)\right|_{0^-} = - (1-f_a) (1-\delta h^\ast_{0^-} \cos \theta)^{{1}/{2}}, & [\text{updip}, -x^\ast_{N_u} < x^\ast < 0],\\ \displaystyle \left.\left( h^\ast\frac{\partial{h^\ast}}{\partial{x^\ast}}\cos \theta -h^\ast \sin \theta \right)\right|_{0^+} = - f_a (1-\delta h^\ast_{0^+} \cos \theta)^{{1}/{2}}, & [\text{downdip}, 0 < x^\ast < x^\ast_{N_d}]. \end{array}\right\}\end{align}

Note that the choice of scalings associated with (2.20a,b) eliminates the factor of $\varGamma$ present in (2.16) but absent directly above. The height continuity and nose conditions, respectively, read as

(2.26)\begin{equation} h^\ast_{0^-}= h^\ast_{0^+} \end{equation}

and

(2.27a,b)\begin{equation} h^\ast_{-x^\ast_{N_u}} =l^\ast_{-x^\ast_{N_u}} =0 \quad \text{and}\quad h^\ast_{x^\ast_{N_d}}=l^\ast_{x^\ast_{N_d}}=0. \end{equation}

Also, the global mass conservation equation (2.19) now reads

(2.28)\begin{align} \int_{-x^\ast_{N_u}}^{x^\ast_{N_d}} (h^\ast+ \left|l^\ast\right|) \,\mathrm{d}x^\ast & =\frac{1}{(1- \delta \cos \theta)^{-{1}/{4}}} \left[ \vphantom{\left(1 + \frac{h^\ast}{l^\ast} \cos \theta \right)}(1-\delta h^\ast_{0} \cos \theta)^{{1}/{4}}\right.\nonumber\\ &\quad \left. +K (1-G^{\prime}) \int_{-x^\ast_{N_u}}^{x^\ast_{N_d}} \left(1 + \frac{h^\ast}{l^\ast} \cos \theta \right) \mathrm{d}x^\ast \right]t^\ast. \end{align}

Equations (2.23)–(2.28) contain four dimensionless variables, namely $\theta$, $\delta$, $K$ and $G'$. Here, $\theta$ defines the left to right asymmetry of the flow; $\delta$ defines the influence of the plume source; $K$ defines the cross-flow resistance at the permeability jump; and $G'$ defines the degree of entrainment experienced by the draining fluid. The first three of these variables are easily estimated for a given porous medium and plume source location. However, $G^{\prime }$, defined as the ratio of the reduced gravity in the lower versus upper layers, remains to be determined. We adopt an empirical approach in estimating $G'$ as described in § 5.1. (In § 5.1 it will be shown that $G'$ depends on the plume source conditions and $\theta$. The ramifications of this observation are deferred to § 5.3 where we draw comparisons between theory and experiment.) With this value (plus $\theta$, $\delta$ and $K$) to hand, (2.23)–(2.28) may be solved numerically, for example using the explicit finite difference scheme described in appendix A.2.

4.1. Experimental set-up

A transparent acrylic box $118\ \textrm {cm}\ \textrm {long} \times 7.6\ \textrm {cm}\ \textrm {wide} \times 60\ \textrm {cm}$ deep filled with spherical glass beads (Potters Industries A Series Premium) and tap water served as the experimental tank. Glass beads were $d_1=3.0\pm 0.2$ mm and $d_2=1.0\pm 0.2$ mm in diameter and were used to construct the two layer porous medium in the manner depicted in figure 7, with the larger beads in the upper layer and the smaller beads in the lower layer. The beads had a density of $1.54\ \textrm {g}\,\textrm {cm}^{-3}$ as compared with $\rho _o=0.998\ \textrm {g}\,\textrm {cm}^{-3}$ for the tap water. The porosity of the tank was measured and found to be $\phi =0.38 \pm 0.05$. Permeabilities were determined based on the empirical relationship proposed by Kozeny and Carman, which is discussed in Dullien (Reference Dullien1979), i.e.

(4.1)\begin{equation} k_i= \frac {d_i^2 \phi^{3}}{180(1-\phi)^2}, \end{equation}

where $i=1,2$. For all experiments conducted here, the value of the permeability ratio, $K=k_2/k_1 \propto d_2^2/d_1^2$ was kept fixed at 0.11.

Figure 7.

Schematic of the set-up for the laboratory experiments.

As depicted in figure 7, source fluid was supplied at the top of the upper layer using a line nozzle that spanned the tank width. The nozzle was designed in such a manner that it produced a uniform flow along its length even at small flow rates (Roes Reference Roes2014). For all of the experiments to be reported upon below, the nozzle was located at the centre of the tank and at a vertical height of $H=18.3$ cm from the permeability jump. Shown in table 1 are the different values of $\theta$ used and the corresponding maximum and minimum layer heights of the upper and lower layers.

Table 1.

Porous media dimensions. (The notation is described in figure 7.)

Dense source fluid supplied through the nozzle was prepared by mixing a precalculated mass of salt into tap water in a 100 l reservoir, whose density was measured to an accuracy of $0.00005\ \textrm {g}\,\textrm {cm}^{-3}$ using an Anton Paar DMA 4500 density meter. Moreover, for the purpose of flow visualization, a small amount of cold-water dye (Procion MX) was added to the saltwater in the reservoir. The dye concentration (determined from the calibration curves of appendix B.1) was small enough that it did not significantly alter the plume source density. The dyed, saltwater was then pumped into an overhead bucket using a hydraulic pump (Little Giant Pump Co.). The bucket contained a cylindrical weir that helped to maintain a constant hydrostatic pressure. Moreover, a manual flow control valve and a flowmeter (Gilmont GV-2119-S-P) were used to ensure a constant flow rate through the nozzle.

4.2. Experimental parameters and flow visualization

Experiments were conducted for four permeability jump angles as listed in table 1 and four source conditions as listed in table 2. For each jump angle, all four source conditions were considered such that we performed 16 experiments in total. It took approximately one hour for each experiment to complete.

Table 2.

Conditions at the plume source.

For flow visualization, experimental images were captured using a Canon Rebel EOS T2i 18.0 PM with an 18–55 mm IS II zoom lens, which collected images every 30 s. Uniform intensity backlighting was achieved using a 3M 1880 overhead projector and by covering the back side of the acrylic box with tracing paper, which served to diffuse the incoming light. The images captured were standard RGB , $720\times 400$ pixels in size and had a resolution of 72 dpi. They were post-processed in MATLAB where all of the images corresponding to a particular experimental set were first cropped to remove unwanted regions outside of the flow domain. The images were then corrected by subtracting away the reference image (collected before the initiation of flow) to remove any systematic spatial variations in the light intensity. Finally, the images were converted to false colour and the pixel intensities were normalized and so ranged from 0 to 1. Images so processed were then compared to categorize various experimental phenomena.

4.3. Experimental observations and interface detection

For qualitatively analysing the experimental results, comparison is made between results obtained with $\theta =15^{\circ }$ and a source flow rate of $Q_s= 1\ \textrm {cm}^{3}\,\textrm {s}^{-1}$, but exhibiting two different source reduced gravities – $g'_s = 20\ \textrm {cm}\,\textrm {s}^{-2}$ and $80\ \textrm {cm}\,\textrm {s}^{-2}$ – corresponding to flow combinations 3 and 4, respectively, from table 2. Representative snapshot images are presented in figure 8(a–c) for $g'_s=20\ \textrm {cm}\,\textrm {s}^{-2}$ and figure 8(d–f) for $g'_s=80\ \textrm {cm}\,\textrm {s}^{-2}$. As expected, the plume, after striking the permeability jump, propagated as an asymmetric pair of gravity currents while simultaneously draining into the lower layer. From our previous discussion, we anticipate that the discharged plume fluid, as it crosses the permeability jump and mixes with lower ambient fluid, will lose its sharp interface. Interestingly, figure 8 suggests that a similar behaviour arises even in the upper layer where dispersion, not accounted for in the model of § 2, results in a blurring of the boundary between discharged plume fluid and ambient fluid. Motivated by this observation, two distinct interfaces were identified in all our experimental images. These are indicated by the red and yellow contours and are defined as the bulk and dispersed interfaces, respectively. (Our method for determining the precise shape of the bulk and dispersed interface is described below.) Within the bulk interface, the pixel intensity is both high and very nearly uniform in space and time. Because pixel intensity is a surrogate for fluid density, we surmise that the density (or reduced gravity) of the fluid within the red contour, which we shall refer to as the bulk fluid, is also approximately uniform. Conversely, the dispersed interface separates the ambient fluid from either discharged plume fluid or from contaminated fluid that has been more substantially diluted through a process of dispersion. The region in-between the bulk and dispersed interfaces shows a non-trivial variation of pixel intensity, suggesting a reduced gravity that varies, not always monotonically, in space. In turn, the fluid located between the red and yellow contours is referred to as the dispersed fluid. Note finally that the bulk interface or red contour was defined by considering pixels having an intensity of 0.85. Meanwhile, the dispersed interface or yellow contour was defined by considering pixels having an intensity of 0.005. These threshold values were chosen such that they gave consistent results while processing all our experimental images.

Figure 8.

False colour experimental images showing discharge along, and draining through, the permeability jump, which here makes an angle $\theta =15^{\circ }$ to the horizontal. Panels (a–c) correspond to flow combination 3 and panels (d–f) correspond to flow combination 4, see table 2. Red, blue and yellow contours are as described in the text. The significance of the dotted lines drawn along the permeability jump in each of panels (a–c) is explained in relation to the inclined time series described below.

In order to analyse the bulk and dispersed interfaces in greater detail, we plot an inclined time series image, which is, in turn, derived from snapshot images. An example appears in figure 9, which is constructed from 390 snapshot images collected from the experiment depicted in figure 8(a–c). Inclined time series images are produced by considering the time evolution of the flow along a sloping line located two to three pixels above the permeability jump. Pixel intensities along this sloping line (shown as dotted lines in figure 8a–c) are extracted as a function of $t$ from the start to the end of the experiment. In figure 9, the abscissa, $x''$, is normalized by the along-jump length where $x''=0$ coincides with the origin, indicated schematically in figure 2. Thus $x''<0$ corresponds to the updip flow while $x''>0$ to the downdip flow. In figure 9, and consistent with figure 8, the red and yellow contours show the nose positions of the bulk and dispersed interfaces, respectively.

Figure 9.

Inclined time series image for the experiment considered in figure 8(a–c). The red contour shows the nose position of the bulk interface whereas the yellow contour shows the nose position of the dispersed interface. The normalized intensity bar indicated on the right shows the transition from pure discharged plume fluid (intensity of 1) to pure ambient fluid (intensity of 0). Time intervals are marked by the horizontal dashed lines to identify different stages of the flow dynamics, which are described in text. The analogue theoretical prediction is indicated by the blue contour.

The flow dynamics of the up- and downdip gravity currents are analysed by tracking in time the position of the noses corresponding to the bulk and dispersed interfaces represented in inclined time series plots such as figure 9. Accordingly, we can describe the evolution of the flow as follows. At $t^* \equiv t/\varPi _t = 0$, plume fluid first reaches the permeability jump. For $t^\ast >0$, gravity currents consisting of discharged plume fluid propagate up- and downdip and simultaneously drain into the lower layer. The nose corresponding to the bulk interface of the downdip gravity current becomes arrested at $t^*=t_1^*$ at which point the downdip run-out length $L_{N_d}$ is reached. By contrast, the updip run-out length $L_{N_u}$ is reached at an earlier point in time because $L_{N_u}$ is so often substantially less than $L_{N_d}$, see figures 3 and 5. For $t^\ast > t^\ast _1$ discharged plume fluid continues to drain into the lower layer during which time the noses of the bulk interface, both up- and downdip, remain fixed. Eventually, given the finite height of our experimental box, contaminated fluid makes contact with the bottom (impermeable) boundary. We designate this point in time as $t_2^*$. For $t^* > t_2^*$, there form a pair of secondary (horizontal) gravity currents at the base of the lower layer. The left and right propagation of these secondary-gravity currents has the effect of remobilizing the previously arrested gravity currents in the upper layer. In figure 9, for instance, such a remobilization occurs at a non-dimensional time just larger than $10^3$. The associated details and dynamics of the flow post-run-out are beyond the scope of the present inquiry and shall be explored in a forthcoming study. In the present analysis, all the results discussed below correspond to dimensionless times strictly below $t^*_2$, this is to ensure the contaminated fluid has not made contact with the bottom boundary.

4.4. Qualitative comparison between theory and experiment for the gravity current shapes

Figure 8 suggests that non-trivially different flow behaviour may be observed depending on $g'_s$. For $g'_s =20\ \textrm {cm}\,\textrm {s}^{-2}$, a larger volume of discharged plume fluid propagates along the permeability jump, whereas for $g'_s = 80\ \textrm {cm}\,\textrm {s}^{-2}$, a larger fraction drains into the lower layer. These observations are consistent with the efficiency curves presented in § 3, see figure 6(b). This leads to longer up- and, more especially, downdip gravity currents in the former case versus the latter. Also shown in each of the images in figure 8 are complementary theoretical results predicted using the (sharp interface) model of § 2. These are plotted as the blue contours where, consistent with the experimental measurements to be summarized in § 5.3, we have considered $G'$ values of 0.43 and 0.66 in panels (a–c) and (d–f), respectively. In general, we find good correspondence with the red contours, particularly in the upper layer. However, and as we will explore in further detail in § 5.3, theory slightly overpredicts the length of the gravity currents, both up- and downdip. In the lower layer, the agreement is typically less robust. At least part of the reason for this discrepancy comes from the general neglect of dispersion by the analytical model which is more prominent in the lower layer because of flow instabilities. With this being the case, we cannot everywhere expect good overlap of the blue and red contours because, in practice, some non-trivial fraction of the dense fluid that drains into the lower layer is mixed into the ambient and so appears between the red and yellow contours. This is especially true in figure 8(d–f).