2.1. Governing equations

We consider an axisymmetric geometry and describe the thin film flow using cylindrical coordinates. In particular, we first examine the inward spreading of a converging gravity current on a permeable boundary of thickness $d$ and permeability $k_{d}$ , see figure 2. The viscosity of the fluid is ${\it\mu}$ , and the density of the fluid is ${\it\rho}_{a}+{\rm\Delta}{\it\rho}$ , where ${\it\rho}_{a}$ is the density of air. The height of the air–liquid interface is denoted as $h(r,t)$ . We neglect the effect of surface tension and fluid mixing at the interface and contact line effects at the moving boundary, as well as viscous stresses at the air–liquid boundary. When the spreading current is long and thin, it is natural to use the lubrication approximation, in which case we neglect the vertical velocity and the fluid has a hydrostatic pressure distribution. Following standard steps, the vertically averaged velocity along the horizontal direction $u(r,t)$ is given by

Next, we assume that the thickness of the permeable boundary $d$ is small compared with the height of the gravity current $h$ , i.e.  $h/d\gg 1$ ; this assumption is only violated in a region very close to the front. Since fluid is allowed to drain from the permeable boundary, the gravitationally driven drainage rate $w(r,t)~(w<0)$ is estimated by

where $k_{d}$ is the effective permeability of the substrate. This model has been employed in several previous studies on the spreading of gravity currents with vertical fluid drainage (see e.g. Acton et al. Reference Acton, Huppert and Worster2001; Pritchard et al. Reference Pritchard, Woods and Hogg2001; Pritchard & Hogg Reference Pritchard and Hogg2002; Pritchard Reference Pritchard2007; Spannuth et al. Reference Spannuth, Neufeld, Wettlaufer and Worster2009; Neufeld et al. Reference Neufeld, Vella, Huppert and Lister2011; Vella et al. Reference Vella, Neufeld, Huppert and Lister2011). The applicability of this model has been discussed in other viscous film spreading problems on porous bases (see e.g. Davis & Hocking Reference Davis and Hocking1999, Reference Davis and Hocking2000). Note that if the permeable boundary is a perforated plate, e.g. as in our laboratory experiments presented below, when estimating $k_{d}$ , the resistance from entering and exiting the holes may also need to be considered if the thickness of the plate is small (see e.g. Sampson Reference Sampson1891; Weissberg Reference Weissberg1962; Dagan, Weinbaum & Pfeffer Reference Dagan, Weinbaum and Pfeffer1982; Jensen, Valente & Stone Reference Jensen, Valente and Stone2014).

Considering fluid drainage, the local continuity equation is written as

where the drainage rate $w$ is given by (2.2). Substituting (2.1) and (2.2) into (2.3), we obtain a PDE that describes the time evolution of the air–liquid interface:

Note that in previous studies on outward-spreading viscous gravity currents over permeable substrates, the governing equations for the fluid–fluid interface are the same as or similar to (2.4) (see e.g. Acton et al. Reference Acton, Huppert and Worster2001; Pritchard et al. Reference Pritchard, Woods and Hogg2001; Spannuth et al. Reference Spannuth, Neufeld, Wettlaufer and Worster2009). However, converging currents are distinguished by distinct boundary conditions, and a different solution procedure also needs to be described, as given below.

To complete the problem, appropriate boundary and initial conditions are needed. We assume that a viscous fluid of height $h_{0}$ initially fills the space between $r=r_{0}$ , i.e. the location of the cylindrical lock gate, and $r=r_{out}$ , i.e. the radius of a cylindrical tank. Thus, the initial condition is given by

The boundary conditions for this problem are given by

where $r_{f}(t)$ denotes the location of the moving front. In order to develop second-kind similarity solutions, where dimensional inspection is not sufficient to determine the functional form of $r_{f}(t)$ , it is necessary to consider the momentum (2.1) and continuity (2.3) equations in their primitive forms. Our main goal is to describe the propagation of the fluid towards $r\rightarrow 0$ , i.e. to determine $r_{f}(t)$ for $k_{d}\neq 0$ . Note that when $k_{d}=0$ , i.e. when the substrate is not permeable, the problem statement reduces to the previous studies of second-kind similarity solutions given by Gratton & Minotti (Reference Gratton and Minotti1990) and Diez et al. (Reference Diez, Gratton and Gratton1992).

2.2. Non-dimensionalization

Let us now define dimensionless variables $\tilde{h}=h/h_{0}$ , $\tilde{r}=r/r_{0}$ , $\tilde{t}=t/t_{c}$ and $\tilde{u} =u/u_{c}$ , with

Then, the equation for the horizontal velocity (2.1) and local continuity equation (2.3) can be rewritten in dimensionless form,

where the constant ${\it\beta}$ , which measures the strength of drainage from the permeable boundary, is defined as

Note that ${\it\beta}$ only involves geometric parameters.

In addition, the governing equation (2.4) for the height of the gravity current $h$ takes the dimensionless form

The dimensionless initial condition for (2.10) is

where $\tilde{r}_{out}\equiv r_{out}/r_{0}$ , and the dimensionless boundary conditions are

where $\tilde{r}_{f}(\tilde{t})\equiv r_{f}(t)/r_{0}$ . Equation (2.10) together with the boundary and initial conditions (2.11) and (2.12) can be solved numerically to provide the time evolution of the fluid–fluid interface $\tilde{h}(\tilde{r},\tilde{t})$ .

2.3. Transformation of the problem

The free-boundary problem of (2.10), subject to (2.11) and (2.12), can be reorganized by introducing new dimensionless variables $f$ and ${\it\tau}$ based on the following transform (see e.g. Murray Reference Murray1989; Pritchard et al. Reference Pritchard, Woods and Hogg2001; Seminara et al. Reference Seminara, Angelini, Wilking, Vlamakis, Ebrahim, Kolter, Weitz and Brenner2012):

Note that ${\it\tau}\rightarrow 0^{+}$ as $\tilde{t}\rightarrow 0^{+}$ , and ${\it\tau}\rightarrow 1/(3{\it\beta})$ as $\tilde{t}\rightarrow +\infty$ . Then, (2.10) is transformed to a well-known nonlinear diffusion equation:

The initial condition (2.11) can be rewritten as

and the boundary conditions (2.12) can be rewritten as

Thus, the change of variables (2.13) eliminates ${\it\beta}$ from the problem statement and effectively reduces the problem to that of propagation of gravity currents over a rigid impermeable boundary, as has been noted in many of the papers listed in table 1. Now (2.14) can be solved numerically subject to initial condition (2.15) and boundary conditions (2.16) using previously reported algorithms (see e.g. Gratton & Minotti Reference Gratton and Minotti1990; Diez et al. Reference Diez, Gratton and Gratton1992; Zheng et al. Reference Zheng, Christov and Stone2014).

2.4. Phase-plane analysis

In order to further study the analytical behaviour of this problem in terms of a second-kind similarity solution, let us now eliminate ${\it\beta}$ from the independent momentum and mass equations (2.8) by introducing a new dimensionless variable $v$ ,

and recall (2.13), so (2.8) can be rewritten as

This pair of equations is independent of ${\it\beta}$ and is defined over the time interval $0\leqslant {\it\tau}\leqslant 1/(3{\it\beta})$ . We note again that ${\it\tau}\rightarrow 0^{+}$ as $\tilde{t}\rightarrow 0^{+}$ , and ${\it\tau}\rightarrow 1/(3{\it\beta})$ as $\tilde{t}\rightarrow +\infty$ . Now, (2.18) recovers the same governing equations for converging gravity currents on an impermeable horizontal boundary (Gratton & Minotti Reference Gratton and Minotti1990). For completeness, we give details on conducting a phase-plane analysis based on (2.18), and seek a self-similar solution of the second kind for the time evolution of the fluid–fluid interface.

We assume that there is a finite time $t_{f}$ , and hence ${\it\tau}_{f}$ from (2.13b ), for the front to reach the origin, and we introduce $s\equiv {\it\tau}_{f}-{\it\tau}$ , which represents the time remaining for the front to arrive at the origin. We define

and substitute (2.19) into (2.18). After some algebra, we obtain the equivalent equations in terms of $V(\tilde{r},s)$ and $F(\tilde{r},s)$ :

Next, we define the similarity variable ${\it\xi}=\tilde{r}/s^{{\it\delta}}$ , where the scaling exponent ${\it\delta}$ is to be determined. Then, $F(\tilde{r},s)$ and $V(\tilde{r},s)$ become functions of only the similarity variable ${\it\xi}$ , i.e. $F=F({\it\xi})$ and $V=V({\it\xi})$ . In this case, (2.20a,b ) take the self-similar form: We can further eliminate ${\it\xi}$ in (2.21) and obtain the autonomous forms Equation (2.22a ) is a first-order nonlinear ordinary differential equation with an unknown parameter ${\it\delta}$ . Given the value of the parameter ${\it\delta}$ , the phase portrait of (2.22a ) can be obtained; different solution curves in the phase portrait represent different self-similar flows subject to different boundary conditions, see figure 3.

Figure 3.

The phase plane for the first-order ordinary differential equation (2.22a ) for different values of ${\it\delta}$ : (a) ${\it\delta}=0.712$ , (b) ${\it\delta}=0.762$ , (c) ${\it\delta}=1.06$ . In (b), the integral curve connects points O and A, thus ${\it\delta}=0.762\dots$ is the solution to the nonlinear eigenvalue problem.

The boundary conditions of a specific physical problem are associated with the critical points in (2.22a ) in the phase plane (see e.g. Gratton & Minotti Reference Gratton and Minotti1990). By letting the numerator and denominator be equal to zero, we obtain three finite critical points in the phase plane of (2.22a ):

Point A describes a state of zero height and finite velocity, which corresponds to the moving front of the gravity current, where $h(r_{f}(t),t)=0$ and $\text{d}r_{f}(t)/\text{d}t\neq 0$ . Point O corresponds to the state when $s=0$ , i.e. ${\it\tau}={\it\tau}_{f}$ (see e.g. Gratton & Minotti Reference Gratton and Minotti1990). Physically, point O corresponds to the state when the front arrives at the origin, i.e. the centre of the cylinder. Thus, the integral curve connecting points O and A in the phase plane describes an inwardly converging viscous gravity current when the front is near the origin.

In order to obtain the integral curve connecting points O and A, we first linearize (2.22a ) near point O, and we obtain the eigenspace:

To solve (2.22a ) subject to (2.23a,b ), we start from an initial condition ( $F_{i},V_{i}$ ) slightly perturbed from point O in the neutral eigendirection given by $\boldsymbol{e}_{2}$ . Given a value of ${\it\delta}$ , the integral curve starting from ( $F_{i},V_{i}$ ) can be calculated numerically using the built-in Mathematica subroutine NDSolve, as shown in figure 3. Then, we adjust the values of ${\it\delta}$ , and we observe that the integral curve connects to point A only when ${\it\delta}=0.762\dots$ , see figure 3, which agrees with the calculation of Gratton & Minotti (Reference Gratton and Minotti1990). Thus, we have found a second-kind self-similar solution for the converging gravity currents with vertical fluid leakage through a permeable substrate.

Table 2.

Summary of experiments on viscous converging gravity currents spreading over impermeable and permeable horizontal substrates. In all of the experiments golden syrup was used as the viscous fluid, and the viscosity ${\it\mu}$ and density ${\it\rho}$ were measured: ${\it\mu}=7.6~\text{Pa}~\text{s}$ and ${\it\rho}=1.47\times 10^{3}~\text{kg}~\text{m}^{-3}$ . The outer radius of the cylindrical tank is $r_{out}=0.147~\text{m}$ , and the location of the lock gate is $r_{0}=0.112~\text{m}$ . The ‘touch-down’ time $t_{f}$ is defined as the time for the front to arrive at the centre of the tank. The dimensionless parameter ${\it\beta}$ is calculated based on (2.9).

3.1. Permeability estimates

In estimating the permeability $k_{d}$ (see (2.2)) of the perforated plates in our experiments, we use

where ${\it\phi}$ denotes the porosity of the substrate and $k_{hole}$ is the permeability of each hole (see e.g. Bear Reference Bear1972; Phillips Reference Phillips1991). The well-known result from the Poiseuille model gives

where $r_{hole}$ represents the radius of the parallel circular holes. When the holes are short, the resistance to flow entering and exiting the holes also needs to be considered besides the wall friction along the hole; for example, see (5) in Jensen et al. (Reference Jensen, Valente and Stone2014), and for more discussions see e.g. Sampson (Reference Sampson1891), Weissberg (Reference Weissberg1962) and Dagan et al. (Reference Dagan, Weinbaum and Pfeffer1982). The interactions of flows entering and exiting each hole may also need to be considered when the spaces between the holes are small; a correction coefficient $G$ is provided for both square and hexagonal arrays, see (16) and figure 5 in Jensen et al. (Reference Jensen, Valente and Stone2014). In our data analysis, in estimating $k_{hole}$ for each hole, we combined (5) and (16) in Jensen et al. (Reference Jensen, Valente and Stone2014), and we obtain:

where $d$ is thickness of the plate and $G$ is the correction coefficient. For example, one perforated plate in our experiments contains parallel circular holes of radius $r_{hole}\approx 1.6~\text{mm}$ in hexagonal arrays, and the centre-to-centre distance between each pair of holes is $L\approx 2.6~\text{mm}$ . The thickness of the plate is $d\approx 1.6~\text{mm}$ . With $G\approx 0.13$ from figure 5 of Jensen et al. (Reference Jensen, Valente and Stone2014), then based on (3.3), the permeability of each hole is $k_{hole}\approx 1.6\times 10^{-7}~\text{m}^{2}$ . The porosity of the plate is ${\it\phi}\approx 0.42$ , and hence the permeability is $k_{d}\approx 6.6\times 10^{-8}~\text{m}^{2}$ . As a comparison, based on the Poiseuille model (3.2) alone, $k_{hole}\approx 3.2\times 10^{-7}~\text{m}^{2}$ and $k_{d}\approx 1.3\times 10^{-7}~\text{m}^{2}$ . As can be seen, both $k_{hole}$ and $k_{d}$ are smaller than the estimates from the Poiseuille model, because of the second term in the denominator of (3.3), which indicates the effect of resistance to flow entering and exiting the holes.

Figure 5.

Time evolution of the front locations from our experiments on viscous converging gravity currents with fluid drainage through the permeable substrate. In (a), the location of the lock gate is 0.112 m, which provides the initial value for all experiments. In (b), $r_{0}$ is the location of the lock gate; $t_{f}$ represents the time for the front to arrive at the centre of the tank, which can be read in (a) as the intercept with the horizontal axis (except for experiments 11 and 12, which stop at a stationary non-zero radial position).

3.2. Front propagation law

We track the locations of the moving front as a function of time, as shown in figure 5. Then, we plot the time evolution of the front locations in the rescaled coordinates, i.e.  $r_{f}/r_{0}$ versus $({\it\tau}_{f}-{\it\tau})/{\it\tau}_{f}$ , see figure 6. Note that the front is not perfectly circular during the propagation process, so we have used an effective radius $r_{eff}$ in our data analysis, defined as $r_{eff}\equiv (A/{\rm\pi})^{1/2}$ , where $A$ denotes the area formed by the moving front. Numerical solutions of the PDE (2.14) subject to initial condition (2.15) and boundary conditions (2.16) are also plotted. The dashed line in figure 6 indicates the slope of 0.762, which comes from the prediction of the second-kind self-similar solution (§ 2.4). Good agreement is observed among the experimental, numerical and theoretical results.

Figure 6.

Experimental results of the front locations versus the prediction of numerical simulations (axes are log–log). The theoretical slope of 0.762 is shown, which is the prediction of the second-kind self-similar solution. We only show the numerical results for experiment 3 since the calculations only show slight differences for experiments 1–10 in table 2 in this plot (the differences are due to different initial conditions). ${\it\tau}_{f}$ is defined based on (2.13b ).

3.3. Will the front reach the origin?

When the initial height is very small in the experiments, e.g. experiments 11 and 12 in table 2, we have observed that the front stops and does not reach the centre of the tank. This result is different from the prediction of the gravity current model without drainage. Note that the height of the gravity current decreases during the experiments because of horizontal spreading and vertical drainage through the permeable substrate; and the vertical drainage prevents the converging front from reaching the centre of the tank for a small initial height.

A numerical simulation has been conducted by solving (2.14) subject to initial condition (2.15) and boundary conditions (2.16). The simulation stops at ${\it\tau}=1/(3{\it\beta})$ , which corresponds to $\tilde{t}\rightarrow +\infty$ . The final location of the front, defined as

can be obtained through the numerical simulations. To demonstrate the idea, we have chosen $r_{out}/r_{0}\approx 1.31$ (based on our experiments), and we show, in figure 7, the relationship between the final front location $r_{\infty }/r_{0}$ and the drainage constant ${\it\beta}$ . When ${\it\beta}\geqslant {\it\beta}_{c}\approx 0.058$ , the drainage effect becomes important, and the front does not reach the origin. As ${\it\beta}$ continues to increase, the effect of drainage becomes more important, and the final front location increases. Experimental measurements of the final front locations have also been plotted in figure 7, and we have observed good agreement. The error bar indicates a 10 % error for the value of the dimensionless parameters $r_{\infty }/r_{0}$ and ${\it\beta}$ . Note that the theoretical predictions are smaller than the experimental observations. This is likely due to the effect of surface tension and contact angle in our experiments when the height of the gravity current is small at late times. We note that the capillary length scale in a typical experiment is $l_{c}\equiv \sqrt{{\it\gamma}/({\rm\Delta}{\it\rho}g)}\approx 2.4~\text{mm}$ , where we use ${\it\gamma}\approx 80~\text{mN}~\text{m}^{-1}$ for golden syrup. Thus, when the film thickness approaches this length scale at late times, capillary effects become more important.

Figure 7.

The final location of the propagating front $r_{\infty }/r_{0}$ depends on the value of drainage constant ${\it\beta}$ , defined in (2.9), which only involves geometric parameters. We have set $r_{out}/r_{0}\approx 1.31$ in the numerical calculation, based on our experiments. When ${\it\beta}>{\it\beta}_{c}\approx 0.058$ , the front does not reach the origin. The observations from experiments 3–6 and 11–12 are also shown in this figure. The error bar represents a 10 % error in the estimates for the final front location $r_{\infty }/r_{0}$ and the drainage constant  ${\it\beta}$ .

In addition, depending on the values of the two dimensionless parameters $r_{0}/r_{out}$ and ${\it\beta}$ , a regime plot can be obtained that describes whether the front reaches the origin or not (figure 8): a ‘touching’ regime in which the front reaches the origin, and a ‘non-touching’ regime in which the front does not arrive at the origin. For a given $r_{out}/r_{0}$ , the front does not reach the origin when the drainage constant ${\it\beta}$ becomes greater than a critical value ${\it\beta}_{c}$ , which is calculated numerically in this figure (the black curve). Information from our experiments has also been plotted in the regime diagram.

Figure 8.

Regime diagram that describes whether the front reaches the origin (i.e. ‘touching’ regime) or not (i.e. ‘non-touching’ regime). For a given $r_{out}/r_{0}$ , the front does not reach the origin when the drainage constant ${\it\beta}$ is greater than a critical value ${\it\beta}_{c}$ , calculated numerically in this plot (the black curve). The observations from experiments 3–6 and 11–12 are also shown in this diagram.

There are other effects that can contribute to the uncertainty of the experiments besides the effect of surface tension and contact angle. For example, we have neglected the effect of leakage through the gap between the substrate and the bottom of the lock gate at the beginning of the experiments; we also neglected the amount of fluid that attaches to the lock gate and hence leaves the system when we remove the lock gate. The error in the measurements of the initial height, the density and viscosity of golden syrup can also contribute to the uncertainty of our experiments. Overall the errors are relatively small and we obtain good agreement with the numerical and theoretical predictions.

4.1. Hele-Shaw channels with varying gap thickness

We consider viscous gravity currents spreading in Hele-Shaw cells with varying gap thickness $b(x)=b_{1}x^{n}$ in the horizontal direction (figure 9). In our previous work, we found a second-kind self-similar solution for $0<n<1$ , when the gravity current propagates towards the origin, where the gap thickness is zero (Zheng et al. Reference Zheng, Christov and Stone2014). Now, let us consider the same problem but with vertical fluid drainage from a thin permeable bottom boundary of permeability $k_{d}$ and thickness $d$ . The origin of the coordinate system is defined as the location of zero gap thickness, i.e.  $b(0)=0$ . The height of the gravity current and the horizontal velocity are denoted by $h_{a}(x,t)$ and $u_{a}(x,t)$ , and initially fluid of height $h_{a0}$ is confined by a lock gate at location $x_{a0}$ . Assuming $d\ll h_{a}$ for the majority part of the current away from the moving front, using Darcy’s law and the local continuity equation, we obtain

Let us now define dimensionless variables $\tilde{h}_{a}\equiv h_{a}/h_{a0}$ , $\tilde{u} _{a}\equiv 12{\it\mu}u_{a}/({\rm\Delta}{\it\rho}gb_{1}^{2}x_{a0}^{2n-1}h_{a0})$ , $\tilde{x}\equiv x/x_{a0}$ , and $\tilde{t}\equiv {\rm\Delta}{\it\rho}gb_{1}^{2}x_{a0}^{2n-2}h_{a0}t/(12{\it\mu})$ . Then, the dimensionless forms of the governing equations (4.1) become where the dimensionless parameter ${\it\beta}_{a}$ measures the significance of vertical fluid drainage from the permeable horizontal boundary, and is defined as

Now we introduce the following transform:

Immediately, we note that ${\it\tau}\rightarrow 0^{+}$ as $\tilde{t}\rightarrow 0^{+}$ , and ${\it\tau}\rightarrow 1/{\it\beta}_{a}$ as $\tilde{t}\rightarrow +\infty$ . In addition, (4.2) are transformed to

Note that (4.5) has the same form as (2.10) in Zheng et al. (Reference Zheng, Christov and Stone2014), which describes the propagation of viscous gravity currents in a Hele-Shaw cell with an impermeable substrate. Then, a second-kind self-similar solution can be derived following § 2.1.2 of Zheng et al. (Reference Zheng, Christov and Stone2014), when the current propagates towards the origin:

where $s\equiv {\it\tau}_{f}-{\it\tau}$ , indicating the time remaining for the front to reach the origin, and the similarity variable is defined as ${\it\xi}\equiv \tilde{x}/s^{{\it\delta}}$ . The scaling exponent ${\it\delta}(n)$ depends on $n$ , and can be obtained (for $0<n<1$ ) by solving a nonlinear eigenvalue problem using a phase-plane analysis, for example, ${\it\delta}(0.2)\approx 1.14$ , ${\it\delta}(0.5)\approx 1.54$ and ${\it\delta}(0.8)\approx 2.95$ from figure 6 of Zheng et al. (Reference Zheng, Christov and Stone2014).

It should also be noted that when the drainage effect becomes important, the front will not be able to reach the origin. We show the final location of the propagating front $x_{\infty }$ , defined as $x_{\infty }\equiv \lim _{\tilde{t}\rightarrow +\infty }x(\tilde{t})$ , for $x_{out}/x_{a0}\approx 1.31$ in figure 10, based on a numerical calculation. Note that $x_{out}$ represents the horizontal span of the channel, and $x_{a0}$ is the location of the lock gate. When the drainage constant ${\it\beta}_{a}$ , defined in (4.3), becomes greater than a critical value ${\it\beta}_{ac}(n)$ (which depends on $n$ , the shape parameter of the channel), the front does not reach the origin. A regime diagram (figure 11) can also be obtained for different values of the drainage constant ${\it\beta}_{a}$ , the dimensionless parameter $x_{out}/x_{a0}$ and the shape parameter $n$ for the channel. Analogously to the discussion in § 3.3, a ‘touching’ regime and a ‘non-touching’ regime can be identified.

Figure 10.

The final location of the propagating front $x_{\infty }/x_{a0}$ depends on the value of drainage constant ${\it\beta}_{a}$ , defined in (4.3), and the geometric parameter $n$ for the shape of the channel. We have set $x_{out}/x_{a0}\approx 1.31$ in the numerical calculation, to be consistent with § 3. When ${\it\beta}>{\it\beta}_{ac}(n)$ , the front does not reach the origin; the value of ${\it\beta}_{ac}(n)$ depends on $n$ , and we have shown the numerical calculations for $n=0.2,0.5,0.8$ in this figure.

Figure 11.

Regime diagram that describes whether the front reaches the origin (i.e. ‘touching’ regime) or not (i.e. ‘non-touching’ regime). For a given $x_{out}/x_{a0}$ , the front does not reach the origin when the drainage constant ${\it\beta}_{a}$ is greater than a critical value ${\it\beta}_{ac}(n)$ (shown as the curves for $n=0.2,0.5,0.8$ ), which also depends on the geometric parameter $n$ for the shape of the channel.

4.2. Porous media with horizontal permeability and porosity gradients

Let us now study the propagation of a viscous gravity current in porous media with permeability and porosity gradients of power-law forms in the horizontal direction: $k(x)=k_{1}x^{n}$ and ${\it\phi}(x)={\it\phi}_{1}x^{m}$ (figure 12). We note that, in practice, $2<n/m<3$ , where $n/m=2$ represents a porous medium of tubular holes, and $n/m=3$ represents a network of intersecting fissures, as mentioned in Phillips (Reference Phillips1991). We allow slow fluid drainage from the thin permeable substrate of permeability $k_{d}$ and thickness $d$ . The origin of the coordinate system is placed at the location of zero permeability and porosity, i.e.  $k(0)=0$ and ${\it\phi}(0)=0$ . We denote the height of the current and the horizontal velocity as $h_{p}(x,t)$ and $u_{p}(x,t)$ , and we assume that initially the fluid of height $h_{p0}$ is confined by a lock gate at $x_{p0}$ . Then, analogously to the development in § 4.1, the governing equations for the fluid flow are written as

Again, we non-dimensionalize (4.7) by defining dimensionless variables $\tilde{h}_{p}\equiv h_{p}/h_{p0}$ , $\tilde{u} _{p}\equiv {\it\mu}u_{p}/({\rm\Delta}{\it\rho}gk_{1}x_{p0}^{n-1}h_{p0})$ , $\tilde{x}\equiv x/x_{p0}$ and $\tilde{t}\equiv {\rm\Delta}{\it\rho}gk_{1}x_{p0}^{n-m}h_{p0}t/({\it\phi}_{1}{\it\mu})$ . Then, the dimensionless form of governing equation (4.7) becomes where a dimensionless parameter ${\it\beta}_{p}$ is defined as

which measures the significance of the vertical fluid drainage. Then, we can introduce the following transform:

and we note that ${\it\tau}\rightarrow 0^{+}$ as $\tilde{t}\rightarrow 0^{+}$ , and ${\it\tau}\rightarrow 1/{\it\beta}_{p}$ as $\tilde{t}\rightarrow +\infty$ . We can also transform equation (4.8) into

Figure 12.

Propagation of a viscous gravity current towards the origin in a porous medium with horizontal permeability and porosity gradients. The substrate is permeable such that fluid drainage occurs in the vertical direction.

Again, we observe that (4.11) has the same form as (2.27) in Zheng et al. (Reference Zheng, Christov and Stone2014), which describes the propagation of a gravity current on an impermeable substrate. Therefore, a second-kind self-similar solution can be obtained following § 2.2.2 of Zheng et al. (Reference Zheng, Christov and Stone2014). Let us denote $s={\it\tau}_{f}-{\it\tau}$ , representing the remaining time for the front to arrive at the origin; then, the self-similar solution is given by

where the similarity variable ${\it\xi}$ is defined as ${\it\xi}\equiv \tilde{x}/s^{{\it\delta}}$ , and the scaling exponent ${\it\delta}(m,n)$ , depending on both $m$ and $n$ , can be obtained by solving a nonlinear eigenvalue problem, for example, ${\it\delta}(1,2.5)\approx 2.42$ and ${\it\delta}(1.5,3.1)\approx 2.78$ .

We also note that as the drainage constant ${\it\beta}_{p}$ increases, the effect of fluid loss through drainage becomes more significant, and eventually the front may not reach the origin. Similar to the analysis in §§ 3.3 and 4.1, we are able to numerically calculate the final location of the propagating front $x_{\infty }$ , defined as $x_{\infty }\equiv \lim _{\tilde{t}\rightarrow +\infty }x(\tilde{t})$ , and we show the results in figure 13(a). We have chosen ( $m=1$ , $n=2.5$ ) and ( $m=1.5$ , $n=3.1$ ) in the numerical calculations as examples; we have also chosen $x_{out}/x_{p0}\approx 1.31$ to be consistent with the analysis in §§ 3.3 and 4.1. In addition, we obtain a regime diagram (figure 13 b) that characterizes whether or not the front reaches the origin, depending on the drainage constant ${\it\beta}_{p}$ , the dimensionless parameter $x_{out}/x_{p0}$ and the values of $m$ and $n$ that describe the horizontal heterogeneity of the porous medium. Note that $x_{out}$ represents the horizontal span of the porous medium and $x_{p0}$ denotes the location of the lock gate.

Figure 13.

In (a), the final location $x_{\infty }$ of the front depends on the drainage constant ${\it\beta}_{p}$ , and the values of $m$ and $n$ that describe the horizontal heterogeneity of the porous medium. We have chosen ( $m=1$ , $n=2.5$ ) and ( $m=1.5$ , $n=3.1$ ) as examples; we have also chosen $x_{out}/x_{p0}\approx 1.31$ to be consistent with §§ 4.1 and 3.3. In (b), the regime diagram is shown and describes whether or not the front arrives at the origin in a porous medium with horizontal permeability and porosity gradients.