2. Previous work

In work by Huppert (Reference Huppert1982), the flow of a viscous gravity current over a flat surface, inclined at angle $\alpha$, was considered (see figure 1). The following equation describing the flow was derived:

(2.1a,b)\begin{equation} \frac{\partial h}{\partial t} ={-} 3\beta \sin \alpha h^2 \frac{\partial h}{\partial s},\quad \beta = \rho g/(3\mu), \end{equation}

where $h(s,t)$ is the fluid thickness at the point $s$, measured down the plane, $\rho$ is the fluid density, $g$ is the acceleration due to gravity and $\mu$ is the dynamic viscosity. The exact, self-similar solution of (2.1a,b) is

(2.2)\begin{equation} \left.\begin{gathered} h(s) = \sqrt{s /\left( 3 t \beta \sin \alpha \right)},\quad 0\leq s \leq s_N,\\ s_N = \left(27 A^2 \beta t \sin \alpha /4\right)^{{1}/{3}},\quad v_N = 9A^2 \beta \sin \alpha/\left(4 s_N^2\right), \end{gathered}\right\} \end{equation}

where $v_N$ is the front speed and $A$ is the total fluid area. We note that (2.1a,b) is approximate. It is based on the assumption that the effect of the hydrostatic pressure gradient is much smaller than the direct action of gravity. This assumption is violated near the front ($s=s_N$) and the tail ($s=0$) for any $\alpha$, while it has reasonable accuracy far from these boundaries. Additionally, the solution (2.2) is clearly inaccurate near the front, with a discontinuity in the fluid thickness $h(s, t)$.

Figure 1.

Schematic representation of the gravity current down a flat inclined surface, being unconstrained (a) and constrained by a wall (b).

Here we examine the influence of these approximations on the accuracy of the solution.

3. Governing equations

In the present section, we derive a general equation, describing the flow of an incompressible fluid over a smooth surface under the action of gravity. The surface is given by the function $\boldsymbol {r}(s)$, where $s$ is a curvilinear coordinate along the surface. A flat inclined surface is a particular case.

The fluid is modelled by a one-dimensional continuum. Each point of the continuum has a velocity $\boldsymbol {v}$, directed parallel to the tangent to the surface at the point:

(3.1a,b)\begin{equation} \boldsymbol{v}(s,t) = v(s, t)\boldsymbol{e}_{\tau}(s),\quad \boldsymbol{e}_{\tau} = \frac{\partial \boldsymbol{r}}{\partial s}, \end{equation}

where $\boldsymbol {e}_{\tau }$ is a unit vector, tangent to the surface at point $s$.

Then the mass balance equation for the fluid has the form

(3.2)\begin{equation} \frac{\partial h}{\partial t} ={-}\frac{\partial}{\partial s}\left(h v\right). \end{equation}

At time $t$, the fluid occupies the domain $s \in [s_T(t); s_N(t)]$, where $s_T(t)$, $s_N(t)$ are the positions of the fluid tail and the fluid front, respectively. We assume that the motion of the fluid is quasi-stationary in the sense of slow, low-Reynolds-number motions. If the fluid is bounded by a wall at $s=s_T$ (as considered by Huppert (Reference Huppert1982)) then

(3.3a,b)\begin{equation} v(s_T)=0,\quad h(s_N) = 0. \end{equation}

Otherwise

(3.4a,b)\begin{equation} h(s_T)=0,\quad h(s_N) = 0. \end{equation}

The total area of the fluid, $A$, is conserved, i.e.

(3.5)\begin{equation} \int_{s_T(t)}^{s_N(t)} h(t)\,\textrm{d} s = A. \end{equation}

Here we implicitly assume that $h$ is much less than the local radius of curvature of the surface. If $h$ is not small then the elementary volume is not equal to $h\,\textrm {d}s$. All further derivations also rely on this assumption.

Neglecting inertial terms in the equation of momentum balance, we obtain

(3.6)\begin{equation} \frac{\partial \boldsymbol{F}}{\partial s} + T\boldsymbol{e}_{\tau} + N \boldsymbol{e}_{n} + \rho \boldsymbol{g} h = 0, \end{equation}

where $\boldsymbol {F}$ is the force per unit width acting on the cross-section of the fluid at point $s$, $\boldsymbol {e}_{n}$ is the unit vector normal to the surface, $T$, $N$ are shear and normal stresses acting on the fluid from the surface and $\boldsymbol {g}$ is the gravitational acceleration.

The system of equations (3.2), (3.5), (3.6) is not closed. To close the system, constitutive equations, relating $\boldsymbol {F}$ and $T$ with $h$ and $\boldsymbol {v}$, are required, while $N$ does not enter the resulting equation. We use the following constitutive relations:

(3.7)\begin{equation} \left.\begin{gathered} \boldsymbol{F} ={-}\rho g h^2 \cos\alpha\boldsymbol{e}_{\tau}/2,\quad T ={-}3\mu v/h,\\ \sin \alpha = \boldsymbol{e}_\tau \boldsymbol{\cdot} \boldsymbol{g}/g,\quad \cos \alpha ={-}\boldsymbol{e}_n \boldsymbol{\cdot} \boldsymbol{g}/g, \end{gathered}\right\} \end{equation}

where $\mu$ is the fluid viscosity and $g = |\boldsymbol {g}|$. According to (3.7), $\boldsymbol {F}$ is proportional to the average hydrostatic pressure multiplied by the fluid thickness, while the expression for $T$ guarantees that (3.9) for $\alpha =0$ reduces to the result of Huppert (Reference Huppert1982).

Substituting (3.7) into (3.6) and multiplying the resulting equation by $\boldsymbol {e}_{\tau }$, we obtain the expression for the fluid velocity:

(3.8)\begin{equation} v = \beta h\left[h \sin \alpha-\frac{1}{2} \frac{\partial}{\partial s} \left(h^2 \cos \alpha\right)\right]. \end{equation}

Then the mass balance equation (3.2) takes the form

(3.9)\begin{equation} \frac{\partial h}{\partial t} = \beta\frac{\partial}{\partial s} \left[\frac{h^2}{2} \frac{\partial}{\partial s} \left(h^2 \cos \alpha\right) - h^3 \sin \alpha\right]. \end{equation}

Note that in general $\alpha$ may be a function of $s$. A different form of (3.9) is

(3.10)\begin{equation} \frac{\partial h}{\partial t} = \beta\frac{\partial}{\partial s} \left[h^3 \left( \frac{\partial h}{\partial s} \cos \alpha - \left(1 + \frac{h\alpha'}{2}\right) \sin \alpha\right)\right],\quad \alpha' = \frac{\textrm{d} \alpha}{\textrm{d} s}. \end{equation}

Here the first term on the right-hand side describes the effect of the hydrostatic pressure gradient (which decreases with increasing $\alpha$), while the second term describes the ‘direct’ influence of the gravitational force.

For the flow along the flat inclined surface ($\alpha ' = 0$), (3.9) takes the form

(3.11)\begin{equation} \frac{\partial h}{\partial t} = \beta\frac{\partial}{\partial s} \left[h^3 \left(\frac{\partial h}{\partial s} \cos \alpha - \sin \alpha\right) \right]. \end{equation}

We note that (3.11) reduces to (2.1a,b) only under the assumption $|\partial h/\partial s| \ll |\tan \alpha |$. The influence of this assumption is discussed in § 6.

4. Near-front asymptotics

In this section, we consider the asymptotic behaviour of the solution near the fluid front ($s \rightarrow s_N$).

At the front, the function $h$ vanishes, as in (3.4a,b). We then assume that in the vicinity of the front $h$ can be represented as a simple power law of the form

(4.1)\begin{equation} h(s) \approx c_1 \left(s_N - s\right)^{\gamma}. \end{equation}

Substituting this expression into (3.8), we obtain

(4.2)\begin{equation} v \approx \beta \left[c_1^2 \left(s_N - s\right)^{2\gamma} \sin \alpha + c_1^3 \gamma \cos \alpha \left(s_N - s\right)^{3\gamma-1} + \frac{c_1^3}{2} \left(s_N - s\right)^{3\gamma} \alpha' \sin \alpha\right]. \end{equation}

The front speed $v_N$ is defined as

(4.3)\begin{equation} v_N = \dot{s}_N = \lim_{s\rightarrow s_N} v. \end{equation}

We assume that the front velocity is bounded. Then from (4.2) and (4.3), it follows that $\gamma \geq 1/3$. Then two regimes are possible. One possibility is that the front is stationary, $v(s_N)=0$. This case is considered in § 5. The second possibility is motion of the front with positive velocity. Then

(4.4a–d)\begin{equation} v(s_N)>0,\quad \gamma =1/3,\quad c_1 = [3v_N/(\beta \cos\alpha_N)]^{{1}/{3}},\quad \alpha_N = \alpha(s_N). \end{equation}

Note that $c_1$ depends on the local inclination of the surface. For $\alpha = 0$, formula (4.4a–d) for $c_1$ coincides with a similar expression obtained for a flat surface by  Huppert (Reference Huppert1982). We also note that a similar asymptotic approach has been employed by Linkov for simulation of hydraulic fracturing (Linkov Reference Linkov2015).

Note that according to (4.2) the front speed is non-negative. This is a somewhat counterintuitive result. For example, consider an unconstrainted current over a flat inclined surface (see figure 1). In this case the fluid has a front and a tail. It is expected that the fluid will move downslope. From the asymptotics, presented above, it follows that the front indeed moves downslope. In contrast, the tail is either stationary (see e.g. the solution of (2.1a,b), obtained in Huppert (Reference Huppert1982)) or it moves up the slope. In the next sections we discuss both front and tail movement.

5. When does the tail not move?

In this section, we address the question of when the tail of fluid on a flat inclined surface is stationary.

Evidently, the tail is stationary if the fluid is bounded by a wall. Then according to (3.3a,b) and (3.8), the boundary condition at the wall is

(5.1)\begin{equation} \left.\frac{\partial h}{\partial s}\right|_{s=0} = \tan \alpha. \end{equation}

The tail may also be stationary in the absence of a constraining wall. If it is so, then the fluid velocity $v(s,t)$, given by (3.8), is positive at least in some vicinity of the tail, i.e.

(5.2)\begin{equation} v = \beta h^2 \cos \alpha \left(\tan \alpha - \frac{\partial h}{\partial s}\right)> 0. \end{equation}

This condition is satisfied only if

(5.3)\begin{equation} \frac{\partial h}{\partial s} < \tan \alpha. \end{equation}

In particular, this condition must be satisfied at $s=s_T$. Then $\partial h/\partial s$ is finite at the tail. Additionally, $h(s_T)=0$ and therefore near the stationary tail

(5.4a–c)\begin{equation} h \approx C(s-s_T)^n,\quad n \geq 1,\quad nC < \tan \alpha. \end{equation}

Thus the tail is stationary if (5.3) is satisfied. Note that the solution (2.2) violates both conditions (5.1) and (5.3). However, we show below that it has reasonable accuracy at large times, determined in § 8. Additionally, we note that $\partial h/\partial s$ depends on time. Therefore even if the initial conditions satisfy (5.3), it can still be violated later on.

6. An approximate solution

In this section, we derive an approximate solution, describing the flow along a flat inclined surface. As in the real experiments (Huppert Reference Huppert1982), the motion of the tail is excluded by putting a wall at $s=0$.

We combine the solution (2.2) with the asymptotics (4.1) and the proper boundary condition at the wall (5.1) as

(6.1)\begin{equation} h=\begin{cases} 2\tan \alpha \sqrt{s_0(s+s_0)}, & 0 \leq s \leq s_*,\\ [3v_N (s_N - s)/(\beta \cos \alpha)]^{{1}/{3}}, & s_{*} \leq s \leq s_N, \end{cases} \end{equation}

where $s_0(t)$, $s_*(t)$ are two auxiliary functions. The function $s_0$ is related to the fluid thickness at the wall as $h(0, t) = 2s_0(t)\tan \alpha$, while $s_*$ is a boundary of the region near the front, where the asymptotics (4.1) is ‘important’. We assume that in this region the effect of pressure gradient (proportional to $|\partial h/\partial s|$) dominates over the effect of direct action of gravity (proportional to $\tan \alpha$). Since $|\partial h/\partial s|$ is infinite at the front, such region always exists. As $s$ decreases, $|\partial h/\partial s|$ also decreases, while $\tan \alpha$ is a constant. Therefore at some point the two terms ($|\partial h/\partial s|$ and $\tan \alpha$) become equal. We assume that the asymptotics is applicable up to this point, i.e.

(6.2)\begin{equation} \left.\frac{\partial h}{\partial s}\right|_{s=s_{*}} ={-}\tan\alpha. \end{equation}

Additionally we require the continuity of the solution at $s_*$, i.e. $h(s_*-0) = h(s_*+0)$. Substituting (6.1) into (6.2) and the continuity condition, we obtain

(6.3)\begin{equation} \left.\begin{gathered} s_{*} = s_N - \sqrt{v_N/\left(9\beta \sin \alpha \tan^2 \alpha \right)}, \\ s_0 = \left[\sqrt{s_*^2+v_N/(\beta \sin \alpha \tan^2\alpha)} - s_*\right]/2. \end{gathered}\right\} \end{equation}

Note that since $v_N$ decreases in time, $s_{*}$ tends to $s_N$, while $s_0$ tends to zero. Therefore intervals at the tail and at the front, where the boundary conditions are important, decrease in time.

Substitution of (6.1) into the global mass balance equation (3.5) yields

(6.4)\begin{equation} v_N/(4\beta \sin \alpha \tan\alpha) + 4\tan\alpha s_0^{{1}/{2}} [(s_* + s_0)^{{3}/{2}}-s_0^{{3}/{2}}]/3 = A. \end{equation}

Equations (6.3) and  (6.4) define the dependence of $v_N(s_N)$ in implicit form. This dependence can be treated as the differential equation with respect to the front position. Therefore the approximate solution $h(s,t)$ of (3.11) for a flat inclined surface is given by (6.1), (6.3) and (6.4).

We estimate the range of applicability of the solution (6.1) as follows. Since the solution is only valid for $s_* > 0$, then from (6.3) it follows that

(6.5)\begin{equation} s_N^2 > v_N/(9\beta \sin\alpha \tan^2\alpha). \end{equation}

We assume that for small $\alpha$, motion of the front may be approximated by the flat surface solution (Huppert Reference Huppert1982), e.g.

(6.6)\begin{equation} s_N = \eta_N (\beta A^3 t)^{{1}/{5}},\quad \eta_N \approx 1.411. \end{equation}

Substituting (6.6) into (6.5), we obtain

(6.7)\begin{equation} \beta A^{{1}/{2}} t> (45\eta_N \sin \alpha \tan^2 \alpha)^{-({5}/{6})}. \end{equation}

Then the solution (6.1) is only valid at large times satisfying the condition (6.7). According to (6.7), this time increases with decreasing inclination as $\alpha ^{-5/2}$.

To obtain the explicit dependence $v_N(s_N)$, we employ the fact that $s_0 \ll s_*$, because $s_* \approx s_N$, $s_0 = h(0)/(2\tan \alpha )$, $h(0) \ll s_N$. Then (6.3) can be approximately replaced by

(6.8)\begin{equation} s_0 \approx v_N/(4 s_* \beta \sin \alpha \tan^2\alpha). \end{equation}

Substituting (6.8) and (6.3) into (6.4), we obtain

(6.9)\begin{equation} \left.\begin{gathered} w_N^2/(24 \tan^2 \alpha) + s_N w_N = 3A/2,\quad w_N = \left(\frac{v_N}{\beta \sin \alpha}\right)^{{1}/{2}}, \\ w_N = 12 \tan \alpha[(s_N^2 + A/(4\tan \alpha))^{{1}/{2}} - s_N]. \end{gathered}\right\} \end{equation}

Since at large times

(6.10)\begin{equation} s_N^2 \gg A/(4\tan \alpha), \end{equation}

the series expansion of (6.9) yields $w_N \approx 3 A/(2 s_N)$ and

(6.11)\begin{equation} v_N \approx 9A^2 \beta \sin \alpha/\left(4 s_N^2\right). \end{equation}

We note that (6.11) coincides with the last formula from (2.2). Using (6.11), we also obtain the approximate expressions for $s_0$, $s_*$ and $h(0)$:

(6.12)\begin{equation} \left.\begin{gathered} s_0/s_N \approx 9A^2/\left(16s_N^4 \tan^2 \alpha\right),\quad s_*/s_N \approx 1 - A/(2s_N^2 \tan \alpha),\\ h(0)/s_N \approx 9A^2/\left(8 s_N^4 \tan \alpha\right). \end{gathered}\right\} \end{equation}

These formulae show how $s_0$ and $h(0)$ tend to zero and $s_*$ tends to $s_N$ with increasing front position $s_N$. Substitution of (6.12) into (6.1) yields the explicit expression for the fluid profile $h(s/s_N)$.

Thus, at large times (at least satisfying (6.7)), the approximate solution of the ‘exact’ equation (3.9) coincides with the exact solution of the approximate equation (Huppert Reference Huppert1982). This fact is demonstrated in the following section. The areas where asymptotics (4.1), (4.4a–d) and the boundary condition (5.1) are important decrease in time ($s_* \rightarrow s_N$, $s_0 \rightarrow 0$). According to (6.12), $s_*$ tends to $s_N$ slower than $s_0$ tends to zero.

7. Evolution of the fluid profile

In this section, we compare predictions of the approximate analytical solution obtained above with numerical results. The solution predicts that for $\alpha > 0$ at large times the fluid profile converges to the square-root shape (2.2). Using numerical simulations, we show that this is indeed the case.

In numerical simulations, we use a uniform spatial mesh with step $\Delta s$. Spatial derivatives in (3.11) are approximated by finite differences. The resulting system of difference equations is solved numerically using the implicit numerical scheme described in Appendix A. Initially, the fluid occupies a square of width $A^{{1}/{2}}$. The following values of dimensionless parameters were used: $\alpha = {\rm \pi}/12$, $\Delta s /A^{{1}/{2}} = 0.003$, $\beta \Delta t A^{{1}/{2}} = 0.033$. Evolution of the fluid profile $h(s,t)$ is shown in figure 2(a). As predicted by the analytical solution (2.2), the area where asymptotics is important shrinks in time. The profile rapidly converges to the square-root-like shape, predicted by both solutions (2.2) and (6.1). The time required for the convergence as a function of the inclination angle $\alpha$ is estimated in the following section. At large times, the solution near the front becomes very steep (see figure 2b). This fact justifies approximation of the solution by the discontinuous function (2.2).

Figure 2.

(a) Normalized fluid profile $h s_N/A$ at short times $\beta t A^{{1}/{2}} = 0$ (blue line), 1.58 (green line), 15.8 (yellow line) and 316.2 (red line). (b) Normalized fluid profile $h s_N/A$ at large times ($\beta t A^{{1}/{2}} = 1.9 \times 10^7$). Numerical solution (blue dots) and analytical solution (2.2) (red line) are shown.

8. What is a ‘small’ inclination angle?

Solutions for flat (Huppert Reference Huppert1982) and inclined (Huppert Reference Huppert1982) surfaces are significantly different. For example, in the flat case, the fluid front propagates as $t^{{1}/{5}}$, while in the inclined case it propagates as $t^{{1}/{3}}$ (see formula (2.2)). Thus two questions arise. Which solution should be used if the inclination angle is small? And what angle can be regarded as small? In this section, we answer these questions by analytical and numerical analysis of (3.11).

Equation (3.11) contains two ‘forces’ driving the flow, namely the pressure gradient proportional to $\partial h/\partial s$ and the direct action of gravity proportional to $\sin \alpha \sim \alpha$. If $|\partial h/\partial s| \gg \alpha$ then the flat surface solution will be valid. However, at large times $\partial h/\partial s$ tends to zero. Therefore the flat surface solution will be valid for a certain time scale only. At larger times, the term proportional to $\alpha$ will always dominate. Therefore at large times the inclined surface solution (2.2) is more relevant.

We estimate the time scale at which the inclination starts to dominate over the pressure gradient as follows. We find a time such that $\partial h/\partial s$ and $\alpha$ are of the same order. Assuming that at short times the flat surface solution is valid, we obtain that the fluid front propagates as $t^{{1}/{5}}$ (see Huppert Reference Huppert1982). Since the total fluid volume is conserved, the fluid thickness decreases as $t^{-({1}/{5})}$. Then the following estimate is valid:

(8.1)\begin{equation} \frac{\partial h}{\partial s} \sim (\beta A^{{1}/{2}}t)^{-({2}/{5})}. \end{equation}

Assuming $\partial h/\partial s \sim \alpha$, we obtain

(8.2)\begin{equation} \beta A^{{1}/{2}} t \sim \alpha^{-({5}/{2})}. \end{equation}

Note that a similar power law follows from (6.7). Formula (8.2) shows that the characteristic time scale at which the inclination becomes dominant rapidly increases with decreasing $\alpha$ and tends to infinity as $\alpha$ tends to zero.

Our numerical simulations, carried out at different inclination angles $\alpha$ from $2^{\circ }$ to $30^{\circ }$, support this conclusion. The simulations show that for all angles the fluid profile converges to the square-root shape. The time, $t_c$, required for convergence strongly depends on $\alpha$. To estimate this dependence we perform simulations with different values of $\alpha$ and compute the difference $\varepsilon$ between the numerical solution $h_{n}$ and the analytical solution $h_{a}$, defined by (2.2). The difference is defined by

(8.3)\begin{equation} \varepsilon(t) = \frac{1}{N(t)}\sum_{i=1}^{N(t)} |h_n^{i} - h_a|/h_n^i, \end{equation}

where $h_n^{i}$ is the fluid thickness at node $i$ and $N(t)$ is the total number of nodes occupied by the fluid at time $t$.

Numerical simulations show that $\varepsilon$ decreases in time (see figure 3a). We define the time required for convergence, $t_c$, as the minimal time satisfying $\varepsilon (t) \leq 0.03$ (see dashed line in figure 3a). Using (8.2), we fit the dependence of $t_c$ on $\alpha$ by the power law

(8.4)\begin{equation} \beta A^{{1}/{2}} t_c = C\alpha^{-({5}/{2})}, \end{equation}

where $C \approx 214$. Figure 3(b) shows that formula (8.4) approximates the numerical results with very good accuracy for all considered inclination angles.

Figure 3.

(a) Relative deviation, $\varepsilon$, of the numerical solution from the square-root solution (2.2) for $\alpha = {\rm \pi}/18$. Dashed line corresponds to $\varepsilon = 0.03$. (b) Time required for convergence of the fluid profile to the square-root shape as a function of the inclination angle $\alpha$. Numerical results (dots) and fitting line (8.4) are shown.

Thus at large times the effect of inclination (direct action of gravity) dominates over the pressure gradient and the fluid profile converges to the square-root shape, predicted by (2.2). The time required for convergence tends to infinity with decreasing inclination angle as $\alpha ^{-5/2}$.