1 Introduction
Sustained inflows of dense fluid into less dense environments generate gravitationally driven motions, termed gravity currents, which propagate along the boundary at the bottom of the domain. These types of flows are common in geophysical and industrial settings, with examples ranging from the continuous discharge of dense gas into the atmosphere, the release of relatively dense industrial waste into waterways and the displacement of one fluid by another of a different density within confined pipes and channels. Gravity currents have been studied for several decades in the laboratory and natural settings (Simpson Reference Simpson1997), through integral and shallow-layer models (Ungarish Reference Ungarish2009) and by direct numerical simulation of the Navier–Stokes equations for fluid flow, coupled to an advection–diffusion equation for a dissolved species that determines the local fluid density (Härtel, Meiburg & Necker Reference Härtel, Meiburg and Necker2000). Each approach has provided valuable understanding into the complex ways in which the relatively dense fluid is transported and drives the motion due to gravitational effects. Much of our insight comes from the investigation of flows generated from ‘lock release’ initial conditions, whereby relatively dense fluid initially at rest behind a lock gate is instantaneously released into a horizontal channel containing less dense fluid. Experiments, numerical simulations and shallow-layer models, which are built upon the assumption of hydrostatic balance, reveal that, provided the effects of viscosity are negligible, the fluid slumps at an initially constant rate of propagation before progressively slowing as the finite volume of released material is spread out along the boundary underlying the fluid domain (Huppert & Simpson Reference Huppert and Simpson1980; Rottman & Simpson Reference Rottman and Simpson1983; Hogg Reference Hogg2006). Subsequently, the motion is retarded by the action of drag (Huppert Reference Huppert1982; Hogg & Woods Reference Hogg and Woods2001) and mixing processes between the two fluids (Johnson & Hogg Reference Johnson and Hogg2013). Of particular importance for applications is that it is possible to identify time scales and length scales at which these transitions in behaviour occur in terms of the properties of the fluid and its initial conditions (see, for example, Huppert & Simpson Reference Huppert and Simpson1980).
Motion due to sustained sources has received less attention and this is the subject of our current investigation. In particular, we investigate the time-dependent motion that occurs as the flow is initiated from a source with the aim of calculating the speed of the front of the flow (i.e. the leading edge of the relatively dense fluid within the less dense ambient fluid). Laboratory experiments have been performed by several investigators including Simpson & Britter (Reference Simpson and Britter1979), Bühler, Wright & Kim (Reference Bühler, Wright and Kim1991) and Hogg, Hallworth & Huppert (Reference Hogg, Hallworth and Huppert2005), with the latter two studies including situations with an externally driven motion in the ambient fluid. These experimental studies, among others, have measured how the speed of the front is related to the density difference between the fluids, and expressed in terms of the reduced gravity and the volume flux per unit width of fluid entering the two-dimensional channel from the source. Modelling studies have sought to capture these motions theoretically; for instance Kranenburg (Reference Kranenburg1993) formulated a two-layer hydraulic model to represent the motion both of the dense fluid and the overlying less dense fluid, while Hogg et al. (Reference Hogg, Hallworth and Huppert2005) introduced a three layer model to account for the motion of the influxed dense fluid and the less dense ambient, as well as a layer of intermediate density formed by mixing processes that occur as the flow evolves. The creation and dynamics of this third, intermediate density layer play a crucial role for modelling aspects of density-driven motion in the presence of adverse ambient flow; the dense currents in these scenarios can be arrested and the layer of intermediate density provides a model for the recirculation of relatively dense fluid. Recently, Shringapure et al. (Reference Shringapure, Lee, Ungarish and Balachandar2013) have performed computational simulations of the motion and interpreted the results using integral models of the steady dynamics. The latter study was among the first to extend direct numerical simulations of the Navier–Stokes equations of gravity currents to ‘open’ flows in which there is a sustained in- and outflow. As opposed to simulations within closed or periodic geometries, such as the lock-exchange configuration, open flows are somewhat more challenging because attention has to be paid to minimising any numerical contamination of the flow by upstream propagation of the outflow boundary conditions. For a broader discussion of the computational challenges associated with modelling gravity currents, we refer the reader to Meiburg, Radhakrishnan & Nasr-Azadani (Reference Meiburg, Radhakrishnan and Nasr-Azadani2015).
In this paper we investigate gravity currents due to a sustained constant influx within a horizontal channel when the motion is such that the height of the layers of dense and less dense fluids are comparable. This means that the motion of the upper less dense fluid cannot be neglected, at least provided the density difference between the two fluids is not large. Importantly then, the motion of the overlying, less dense layer affects the streamwise pressure gradient and because the flow is generated by a sustained source, it does not become increasingly shallow in time, which would allow the motion of the upper layer to be neglected after a sufficient time. Thus this flow configuration is not amenable to study using a single-layer hydrostatic model, in which Hoult (Reference Hoult1972), Gratton & Vigo (Reference Gratton and Vigo1994) and Hogg et al. (Reference Hogg, Hallworth and Huppert2005), amongst others, have calculated the flow speed as a function of the governing parameters. We also note that the flows driven by a sustained influx differ from those arising from lock-release initial conditions, because they do exhibit the same progressive thinning, which ultimately enhances the effects of hydraulic resistance. Our approach extends the two-layer shallow-water models of Kranenburg (Reference Kranenburg1993) and identifies solutions with bores in addition to those that are spatially continuous. Our model neglects mixing between the two fluid layers and drag from either the underlying boundary or the interfacial processes. Such processes may become important, or even dominant, after sufficiently long times and distances of propagation (see, for example, Hogg & Woods Reference Hogg and Woods2001 and Johnson & Hogg Reference Johnson and Hogg2013).
Our methods for calculating the solutions are quasi-analytical: since the motion is due to sustained boundary conditions and the governing equations are hyperbolic, the method of characteristics may be applied to compute the evolution of the dependent variables. In the Boussinesq regime, the governing equations admit analytical expressions for the Riemann invariants that are conserved along characteristics and it will be shown that the velocity and height fields of these currents can be found by the solution of purely algebraic equations. For the non-Boussinesq regime, we could not find analytical expressions for the Riemann invariants and the construction of these solutions entails the numerical integration of a first-order differential equation along characteristics (Rotunno et al. Reference Rotunno, Klemp, Bryan and Muraki2011; Ungarish Reference Ungarish2011b ). However, it will be shown that the same range of solution types exist for the non-Boussinesq regime. Equivalently, these flows may be thought of as similarity solutions to the governing equations because there is a simple gearing between spatial and temporal scales (Kranenburg Reference Kranenburg1993; Gratton & Vigo Reference Gratton and Vigo1994). The type of solution depends upon the values of three dimensionless parameters that measure the magnitude of the volume flux per unit width delivered in the channel at the source, the source Froude number and the relative excess density. Importantly, we find that not all source conditions are consistent with the establishment of the two-layer motion. For some situations, typically when the flux is large, the dense fluid undergoes a near-source transition and fills the channel, with a subsequent further transition downstream as the two-layer motion is re-established.
The paper is structured as follows. First we formulate the governing equations and boundary conditions for the two-layer shallow-water model and identify its characteristic structure (§ 2). Here we treat the front of the motion as a moving discontinuity in which its velocity and height are related according to the model proposed by Benjamin (Reference Benjamin1968). Other theoretical models for the non-hydrostatic motion that occurs at the front have been proposed (Borden & Meiburg Reference Borden and Meiburg2013a ), as well as empirical relationships between the speed of the front and the relative depth (Huppert & Simpson Reference Huppert and Simpson1980). We show that if alternative models are adopted then the same types of solutions emerge but that the boundaries between the parameter values at which solutions change are different (see appendix A). We then analyse the types of motion for Boussinesq currents (§ 3). In this regime the Riemann invariants on the characteristics of the two-layer motion are available analytically and so the determination of the solutions is straightforward. We report results of numerical simulations of the Navier–Stokes equations and compare the ‘bulk’ properties with the predictions from the idealised shallow-layer models (§ 4). We also extend the shallow-layer theory to non-Boussinesq currents in § 5 and show that, although the same types of solutions may be found, the parameter values that demark the transitions between regimes change.
Dambreak flow, in which fluid is set in motion from the collapse of a lock of infinite extent, shares many similarities with the motion due to a sustained inflow, because the unsteady flow that develops exhibits the identical gearing between spatial and temporal scales. Dambreak flow, modelled using the single-layer shallow-water equations, has become an important paradigm in which to investigate unsteady and spatially evolving flow. Here, in appendix B, we apply our methodology to dambreak flow of Boussinesq fluids (see Ungarish Reference Ungarish2009). We extend the earlier works to show that many of the results are available analytically and thus this study provides the analogous results for Boussinesq two-layer dambreak flow to those derived for the single-layer shallow-water equations (Ritter Reference Ritter1892; Whitham Reference Whitham1974).
2 Shallow-layer formulation
We analyse the motion of a sustained current of relatively dense fluid through a less dense environment within a horizontal channel of height
$H$
. The motion is driven from a source that delivers fluid of density
${\it\rho}_{1}$
at a volume flux per unit width
$Q$
, while the environmental fluid is of density
${\it\rho}_{2}$
. The thickness of the flowing dense layer is denoted by
$h_{1}$
and its depth-averaged velocity is given by
$u_{1}$
. The thickness of the overlying layer of less dense fluid, flowing above the dense layer is
$h_{2}$
and its depth-averaged velocity is denoted by
$u_{2}$
(see figure 1). Conservation of mass in each layer is then given by
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial h_{1}}{\partial t}+\frac{\partial }{\partial x}(u_{1}h_{1})=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial h_{2}}{\partial t}+\frac{\partial }{\partial x}(u_{2}h_{2})=0, & \displaystyle\end{eqnarray}$$
where the coordinate axes are aligned so that
$x$
is horizontal and
$z$
is vertical. In this formulation, we have neglected mixing processes between the two layers which would generate a zone of intermediate density and which could become dynamically important sufficiently far downstream from the source when the mixed layer has grown sufficiently thick. The motion is assumed to evolve so that horizontal length scales far exceed vertical length scales and thus the pressure is hydrostatic to leading order. Denoting the gravitational acceleration by
$g$
, we find that
$$\begin{eqnarray}p=p_{T}(x,t)+\left\{\begin{array}{@{}ll@{}}{\it\rho}_{2}g(H-z),\quad & h_{1}<z<H,\\ {\it\rho}_{2}gh_{2}+{\it\rho}_{1}g(h_{1}-z),\quad & 0<z<h_{1},\end{array}\right.\end{eqnarray}$$
where
$p_{T}(x,t)$
is the as yet unknown pressure at
$z=H$
. Then following Long (Reference Long1955), Rottman & Simpson (Reference Rottman and Simpson1983) and others, we may write the inviscid momentum balance in each layer as
$$\begin{eqnarray}\displaystyle & \displaystyle {\it\rho}_{1}\left(\frac{\partial u_{1}}{\partial t}+u_{1}\frac{\partial u_{1}}{\partial x}\right)+\frac{\partial p_{T}}{\partial x}+{\it\rho}_{1}g\frac{\partial h_{1}}{\partial x}+{\it\rho}_{2}g\frac{\partial h_{2}}{\partial x}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\it\rho}_{2}\left(\frac{\partial u_{2}}{\partial t}+u_{2}\frac{\partial u_{2}}{\partial x}\right)+\frac{\partial p_{T}}{\partial x}=0. & \displaystyle\end{eqnarray}$$
Here we have neglected drag at the channel boundaries, noting that it could become a dynamically significant process sufficiently far from the source (Hogg & Woods Reference Hogg and Woods2001). The sustained source of dense fluid leads to the boundary conditions
Figure 1. Configuration of the two-layer gravity current driven by a sustained source of dense fluid in a horizontal channel. In this sketch, the layer of dense fluid includes a bore connecting two regions within one of which the interface is of constant height and within the other it is thinning towards the front.
$$\begin{eqnarray}u_{1}h_{1}=Q\quad \text{and}\quad u_{2}h_{2}=0\quad \text{at }x=0.\end{eqnarray}$$
Since
$h_{1}+h_{2}=H$
, we deduce from (2.1) and (2.2) that
$$\begin{eqnarray}u_{1}h_{1}+u_{2}h_{2}=Q.\end{eqnarray}$$
We may then eliminate the unknown pressure
$p_{T}(x)$
by subtracting (2.4) and (2.5) and further use (2.7) to eliminate
$u_{2}$
. Writing
${\it\rho}_{1}={\it\rho}_{2}(1+S)$
with
$S>0$
, we find that momentum balance in the lower layer is expressed by (cf. Rotunno et al.
Reference Rotunno, Klemp, Bryan and Muraki2011, Ungarish Reference Ungarish2011b
)
$$\begin{eqnarray}\displaystyle & & \displaystyle \frac{\partial u_{1}}{\partial t}+\left(u_{1}-\frac{2h_{1}(Hu_{1}-Q)}{(H+S(H-h_{1}))(H-h_{1})}\right)\frac{\partial u_{1}}{\partial x}\nonumber\\ \displaystyle & & \displaystyle \quad +\left(\frac{H-h_{1}}{H+S(H-h_{1})}\right)\left(gS-\frac{(u_{1}H-Q)^{2}}{(H-h_{1})^{3}}\right)\frac{\partial h_{1}}{\partial x}=0.\end{eqnarray}$$
At this point it is convenient to introduce dimensionless variables. We non-dimensionalise using the length scale
$H$
and the velocity scale
$(gSH)^{1/2}$
and replace
$x$
and
$t$
by their dimensionless counterparts, while we define the dimensionless velocity and height of the lower flowing layer by
$$\begin{eqnarray}u=u_{1}(gSH)^{-1/2}\quad \text{and}\quad h=h_{1}/H.\end{eqnarray}$$
This choice of scales introduces a dimensionless parameter,
$\hat{Q}$
, which measures the magnitude of the flux per unit width in the channel, given by
$$\begin{eqnarray}\hat{Q}=\frac{Q}{(gSH^{3})^{1/2}}.\end{eqnarray}$$
Henceforth, unless stated otherwise, we deal with dimensionless variables.
We now treat dimensionless versions of (2.1) and (2.8) as the coupled governing equations for the system, from which the other dependent variables may be deduced. Potentially they form a hyperbolic system with characteristics
$\text{d}x/\text{d}t={\it\lambda}_{\pm }$
, where
$$\begin{eqnarray}{\it\lambda}_{\pm }=u-{\it\alpha}\pm \sqrt{{\it\alpha}^{2}+h(1-{\it\beta})},\end{eqnarray}$$
and where it is convenient to define
$$\begin{eqnarray}{\it\alpha}=\frac{h(u-\hat{Q})}{(1-h)(1+S(1-h))}\quad \text{and}\quad {\it\beta}=1-\frac{1-h}{1+S(1-h)}\left(1-\frac{(u-\hat{Q})^{2}}{(1-h)^{3}}\right).\end{eqnarray}$$
Real-valued characteristics exist and the system is strictly hyperbolic provided
${\it\alpha}^{2}+h(1-{\it\beta})>0$
, which corresponds to
$$\begin{eqnarray}\frac{h(1-h)}{(1+S(1-h))^{2}}\left(1+S(1-h)-(1+S)\left(\frac{u-\hat{Q}}{1-h}\right)^{2}\right)>0.\end{eqnarray}$$
As observed by Boonkasame & Milewski (Reference Boonkasame and Milewski2011), this condition of hyperbolicity may be compactly written in terms of dimensional variables by
$$\begin{eqnarray}\frac{S{\it\rho}_{2}g}{(u_{2}-u_{1})^{2}}\left(\frac{h_{2}}{{\it\rho}_{2}}+\frac{h_{1}}{{\it\rho}_{1}}\right)\geqslant 1.\end{eqnarray}$$
The loss of hyperbolicity in this two-layer system will be shown to introduce an important constraint on the possible types of gravity current motion. Loss of hyperbolicity may be interpreted as an instability of the two-layer motion in the form of Kelvin–Helmholtz instability at a sufficiently long wavelength so that it is not filtered out by the assumption of hydrostatic balance (Boonkasame & Milewski Reference Boonkasame and Milewski2011).
Riemann invariants may be evaluated along each family of characteristics (Long Reference Long1955). On the characteristic
$\text{d}x/\text{d}t={\it\lambda}_{\pm }$
, they satisfy
$$\begin{eqnarray}\frac{\text{d}u}{\text{d}h}=-\frac{1}{h}\left({\it\alpha}\pm \sqrt{{\it\alpha}^{2}+(1-{\it\beta})h}\right).\end{eqnarray}$$
In the Boussinesq regime
$(S\ll 1)$
, this differential equation may be integrated analytically to yield
$$\begin{eqnarray}R_{\pm }=\sin ^{-1}\left(\frac{u-\hat{Q}}{1-h}\right)\mp \sin ^{-1}(1-2h),\end{eqnarray}$$
and the system evolves such that
$R_{\pm }$
is constant on the characteristic
$\text{d}x/\text{d}t={\it\lambda}_{\pm }$
(Long Reference Long1955).
A second source condition may be applied to the system provided both characteristics propagate into the domain
$x>0$
, a condition given by
${\it\lambda}_{\pm }>0$
. We note that for two-layer flow, the condition of criticality
$({\it\lambda}_{-}=0)$
corresponds in general to
$$\begin{eqnarray}\frac{(\hat{Q}-hu)^{2}}{(1-h)^{2}}+(1+S)\frac{u^{2}}{h}=1.\end{eqnarray}$$
Thus for the source condition analysed in this study, we write
$$\begin{eqnarray}u=F_{0}\sqrt{h}\quad \text{at }x=0,\end{eqnarray}$$
and this corresponds to a supercritical source if
$F_{0}>(1+S)^{-1/2}$
. It is noteworthy that for the boundary condition studied here, this condition is identical to the criterion of supercriticality in a single-layer shallow-water model. If the source is subcritical
$(F_{0}<(1+S)^{-1/2})$
then only one boundary condition (2.6) may be imposed at the source. Finally we need a dynamic condition at the front of the gravity current, which encompasses the non-hydrostatic behaviour there. Applying the model of Benjamin (Reference Benjamin1968) as a jump condition in a frame of reference moving with the front, we find
$$\begin{eqnarray}(u_{N}-\hat{Q})^{2}=h_{N}\frac{(2-h_{N})(1-h_{N})}{1+h_{N}}\quad \text{at }x=x_{N}(t),\end{eqnarray}$$
where
$h_{N}=h(x_{N},t)$
and
$u_{N}=u(x_{N},t)$
. This dynamic condition at the front (2.19) may only be applied provided the speed of propagation that it demands is less than the largest, local characteristic speed
$(u_{N}<{\it\lambda}_{+}(u_{N},h_{N}))$
. If this criterion is not satisfied, then the flow is ‘choked’ and the motion cannot deliver fluid rapidly enough to match the required front speed (Ungarish Reference Ungarish2009, Reference Ungarish2011a
). The onset of ‘choking’ corresponds to the time at which the fastest speed of a characteristic first equals the dynamic front speed given by (2.19),
$u_{N}={\it\lambda}_{\pm }(u_{N},h_{N})$
, which implies that
$$\begin{eqnarray}(u_{N}-\hat{Q})^{2}=(1-h_{N})^{3}.\end{eqnarray}$$
The onset of ‘choked’ flow therefore occurs at
$$\begin{eqnarray}(1-h_{N})^{3}=\frac{h_{N}(2-h_{N})(1-h_{N})}{1+h_{N}},\end{eqnarray}$$
which has solution
$h_{N}=h_{c}\equiv 2\sin ({\rm\pi}/18)=0.3473$
and importantly from (2.20),
$u_{N}=\hat{Q}+0.5273$
. The consequence of choking is that the fluid motion satisfies (2.19) when the height of the front is sufficiently shallow
$(h_{N}<h_{c})$
and that the height of the front cannot exceed
$h_{c}$
. We may therefore combine these conditions to write the condition at the front as
$$\begin{eqnarray}h_{N}\leqslant h_{c}\quad \text{and}\quad (u_{N}-\hat{Q})^{2}=\mathscr{F}(h_{N})\equiv \frac{h_{N}(2-h_{N})(1-h_{N})}{1+h_{N}}.\end{eqnarray}$$
It is important to calculate the energy fluxes transported across the front to ensure that there is net dissipation. This is most conveniently evaluated in a frame moving with the front. Then the difference between the dimensionless energy fluxes transported away from and towards the front is given by
$$\begin{eqnarray}\dot{\mathscr{D}}=\left[\frac{(\hat{Q}-u_{N})^{2}}{2(1-h_{N})^{2}h_{N}}-1\right](\hat{Q}-u_{N})h_{N}.\end{eqnarray}$$
Substituting for the combined dynamic condition at the front, (2.19), we see that there is always dissipation
$(\dot{\mathscr{D}}>0)$
. We emphasise that the results (2.20)–(2.23) are valid for both Boussinesq and non-Boussinesq regimes.
In what follows we examine the unsteady solutions that develop for different values of the three dimensionless parameters that characterise the system, namely the dimensionless volume flux per unit width,
$\hat{Q}$
, the Froude number at the source,
$F_{0}$
and the dimensionless excess density,
$S$
. The height and velocity of the lower layer at the source is then given by
$u=u_{0}\equiv (\hat{Q}F_{0}^{2})^{1/3}$
and
$h=h_{0}\equiv (\hat{Q}/F_{0})^{2/3}$
(provided the source is supercritical). Furthermore the condition that the source must correspond to a strictly hyperbolic system of equations demands that (2.13) is satisfied. This in turn implies that for hyperbolicity,
$$\begin{eqnarray}\hat{Q}\leqslant \frac{F_{0}}{\left(\displaystyle F_{0}^{2}+\frac{S}{1+S}\right)^{3/2}}.\end{eqnarray}$$
In what follows we first analyse Boussinesq flows
$(S\ll 1)$
in detail (§ 3), drawing out the different flow types that may occur for different values of the dimensionless parameters
$(\hat{Q},F_{0})$
. This analysis is straightforward due to the existence of the Riemann invariants (2.16) and so many of the results are available as simple analytical expressions. We then analyse non-Boussinesq flows
$(S=O(1))$
and show that the possible types of solutions are identical to the Boussinesq case, but that the boundaries between the different regimes change (§ 5). Our analysis employs Benjamin’s formula for the dynamic condition at the front of the gravity current (2.19); we note that for Boussinesq currents Borden & Meiburg (Reference Borden and Meiburg2013a
) have developed an alternative condition, while Huppert & Simpson (Reference Huppert and Simpson1980) propose an empirical front condition. The consequences of adopting their condition for the sustained flows studied in the paper are examined in § A.1.
3 Boussinesq currents
$(S\ll 1)$
3.1 Uniform solutions and choked-uniform solutions
Uniform solutions have no spatial variations in their velocity and height fields between the source and the front. Thus they exist when the source and frontal Froude numbers are identical. In terms of the dimensionless variables, by substituting
$h_{N}=h_{0}$
and
$u_{N}=u_{0}$
, this condition demands
$$\begin{eqnarray}\hat{Q}^{2}(1-h_{0})^{2}=h_{0}^{2}\mathscr{F}.\end{eqnarray}$$
This expression may be viewed as a relationship between
$F_{0}$
and
$\hat{Q}$
for which uniform conditions exists; it is plotted in figure 2 and is denoted by
$F_{0}=F_{U}(\hat{Q})$
. Notably, a uniform current in a deep ambient
$(\hat{Q}\ll 1)$
occurs when
$F_{0}=\sqrt{2}$
and this recovers the result from a single-layer model of gravity current motion. Example solutions for the height and velocity fields are plotted in figure 3. Choked flow occurs when
$h_{N}=h_{c}$
and
$u_{N}=u_{c}\equiv Q+(1-h_{c})^{3/2}$
and this condition corresponds to
$\hat{Q}=\hat{Q}_{c}\equiv h_{c}(1-h_{c})^{1/2}=0.2806$
,
$F_{0}=F_{0c}\equiv h_{c}^{-1/2}(1-h_{c})^{1/2}=1.3709$
. This point is plotted in figure 2.
Figure 2. Regimes of Boussinesq solutions, classified in terms of the dimensionless flux per unit width,
$\hat{Q}$
, and the source Froude number,
$F_{0}$
. Uniform solutions are only found along
$F_{0}=F_{U}$
and choked-uniform solutions exist for
$F_{0}=F_{Uf}$
. When
$F_{0}<F_{U}\cup F_{Uf}$
the solutions feature rarefactions, whereas for
$F_{0}>F_{U}\cup F_{Uf}$
the solutions feature shocks. Both rarefaction and shock solutions can form choked flows; the boundary between choked and unchoked flows is given by
$\mathscr{C}_{R}$
and
$\mathscr{C}_{S}$
, respectively. Additionally the curve
$\mathscr{C}_{h}$
is plotted; this curve gives the parameter values for which hyperbolicity is lost. Also plotted are the points
$(\hat{Q}_{c},F_{0c})$
and
$(\hat{Q}_{h},F_{0h})$
, which are the parameter values for which uniform flows first become choked and for which two-layer choked-uniform flows lose their hyperbolic nature and transition to a channel-filling motion.
Beyond these parameter values at which choked flow first occurs, we may no longer find flow solutions that are spatially uniform between source and front. Instead, the flows must thin and accelerate to meet the requirement that the height at the front is given by the critical height,
$h_{c}$
. In this subsection we also construct flow solutions that are uniform within a region attached to the source and then adjust to the front conditions throughout a region attached to the front. These two parts of the solution correspond to a uniform region
$0<x/t<y_{c}$
within which the velocity and height fields are given by
$u_{0}$
and
$h_{0}$
, respectively, and a rarefaction for
$y_{c}<x/t<y_{N}$
, which encompasses an expansion fan of
${\it\lambda}_{+}$
characteristics. Thus the trailing characteristic of this fan is given by
$y_{c}=x/t$
, where
$$\begin{eqnarray}y_{c}=u_{0}-{\it\alpha}_{0}+\sqrt{{\it\alpha}_{0}^{2}+h_{0}(1-{\it\beta}_{0})},\end{eqnarray}$$
and where
${\it\alpha}_{0}$
and
${\it\beta}_{0}$
correspond to
${\it\alpha}$
and
${\it\beta}$
in (2.12) evaluated at
$u=u_{0}$
and
$h=h_{0}$
. At the front of the current
$h_{N}=h_{c}$
and
$u_{N}=u_{c}$
(and so
$y_{N}=u_{c}$
) and the Riemann invariant corresponding to the
${\it\lambda}_{-}$
characteristics is given by
$R_{-}(u_{c},h_{c})$
. Throughout the expansion fan, this Riemann invariant takes the same constant value and so
$$\begin{eqnarray}R_{-}(u,h)\equiv \sin ^{-1}\left(\frac{u-\hat{Q}}{1-h}\right)+\sin ^{-1}(1-2h)=R_{-}(u_{c},h_{c})\end{eqnarray}$$
for
$y_{c}<x/t<y_{N}$
. Then since the
${\it\lambda}_{+}$
characteristics are straight lines in the
$(x,t)$
plane, we deduce that
$$\begin{eqnarray}\frac{x}{t}=u-{\it\alpha}+\sqrt{{\it\alpha}^{2}+h(1-{\it\beta})},\end{eqnarray}$$
and together (3.3) and (3.4) provide an implicit solution for the height and velocity fields within the expansion fan. This form of solution may only be constructed along a curve of parameter values in the
$(\hat{Q},F_{0})$
plane given by
$$\begin{eqnarray}R_{-}(u_{0},h_{0})=R_{-}(u_{c},h_{c}).\end{eqnarray}$$
This curve is denoted
$F_{0}=F_{Uf}(\hat{Q})$
and is plotted in figure 2. It intersects with the curve corresponding to the loss of hyperbolicity at
$$\begin{eqnarray}\hat{Q}=\hat{Q}_{h}\equiv {\textstyle \frac{1}{2}}(1+(4h_{c}-1)(1-h_{c})^{1/2}),\end{eqnarray}$$
and
$F_{0}=F_{0h}\equiv \hat{Q}_{h}^{-1/2}$
; these evaluate to
$(\hat{Q}_{h},F_{0h})=(0.6572,1.2335)$
. We plot this point in figure 2 and solutions at various points along the curve
$F_{U}\cup F_{Uf}$
in figure 3.
Figure 3. The height,
$h$
, and velocity,
$u$
, as functions of
$x/t$
for flows that are uniform or are uniform and choked with a rarefaction at the front. The profiles corresponds to parameter values
$(\hat{Q},F_{0})$
: (i) (0.1, 1.372); (ii) (0.2, 1.366); (iii)
$(\hat{Q}_{c},F_{0c})\equiv (0.2806,1.3709)$
; (iv) (0.4, 1.363); (v) (0.5, 1.333); and (vi) (0.6, 1.279).
3.2 Rarefactions
When the source Froude number is less than the Froude number for which a uniform state would occur
$(F_{0}<F_{U}\cup F_{Uf})$
, the flow must accelerate to match the frontal boundary condition. The form of the solution depends upon whether the flow is choked.
First we examine unchoked flows
$(h_{N}<h_{c})$
and in this case, the solution is continuous and comprises three portions. In the ranges
$0<x/t<y_{a}$
and
$y_{b}<x/t<y_{N}$
, where
$y_{a}$
,
$y_{b}$
and
$y_{N}$
are to be determined, the flow is uniform and exhibits no spatial variations. However between these regions there is a rarefaction, corresponding to an expansion fan of
${\it\lambda}_{-}$
characteristics. The uniform region attached to the source is given by
$$\begin{eqnarray}u=u_{0},\quad h=h_{0},\quad \text{for }0<x/t<y_{a},\end{eqnarray}$$
where the extent of the region,
$y_{a}$
, is given by
$$\begin{eqnarray}y_{a}={\it\lambda}_{-}(u_{0},h_{0}).\end{eqnarray}$$
Within the rarefaction
$(y_{a}<x/t<y_{b})$
, the Riemann invariant,
$R_{+}$
, is constant and so
$$\begin{eqnarray}R_{+}(u,h)=R_{+}(u_{0},h_{0}).\end{eqnarray}$$
For
$y_{b}<x/t<y_{N}$
, there is a uniform region attached to the front of the flow and here
$$\begin{eqnarray}u=u_{N},\quad h=h_{N}\quad \text{and}\quad (u_{N}-\hat{Q})^{2}=\mathscr{F}(h_{N}).\end{eqnarray}$$
Thus using (3.9), evaluated at
$u=u_{N}$
and
$h=h_{N}$
, and (3.10), we have coupled algebraic equations for the height and velocity at the front of the current. These solutions then determine the edge of the expansion
$x/t=y_{b}$
, where
$$\begin{eqnarray}y_{b}={\it\lambda}_{-}(u_{N},h_{N}).\end{eqnarray}$$
To complete the solution we calculate the variation of the velocity and height fields within the expansion fan and this corresponds to coupling (3.9) with the characteristic equation
$$\begin{eqnarray}y\equiv \frac{x}{t}={\it\lambda}_{-}(u,h).\end{eqnarray}$$
Typical profiles are plotted in figure 4. Importantly we may evaluate when this form of solution becomes affected by the ‘choked’ flow constraint, namely that the fastest characteristic speed is less than the required front speed. The onset of choked flow occurs for parameters
$\hat{Q}$
and
$F_{0}$
such that
$$\begin{eqnarray}R_{+}(u_{c},h_{c})=R_{+}(u_{0},h_{0}).\end{eqnarray}$$
This curve, denoted by
$\mathscr{C}_{R}$
, is plotted in figure 2.
Figure 4. The height,
$h$
, and velocity,
$u$
, as functions of
$x/t$
for dimensionless fluxes
$\hat{Q}=0.1,0.2,0.3,0.4,0.5,0.6$
and (a,b)
$F_{0}=1.2$
, (c,d)
$F_{0}=1$
. Increasing values of
$\hat{Q}$
produce successively deeper and faster flowing currents in both (a,b) and (c,d).
When the flow is choked, the construction of rarefaction solutions becomes more complicated because there are now two expansion fans within which the current adjusts to match the critical height,
$h_{c}$
, and velocity,
$u_{c}$
, at the front. The solution thus comprises a uniform region attached to the source, (
$u=u_{0}$
and
$h=h_{0}$
) for
$0<x/t<y_{a}$
; an expansion fan of
${\it\lambda}_{-}$
characteristics for
$y_{a}<x/t<y_{b}$
; a uniform region (
$u=u_{m}$
and
$h=h_{m}$
, where
$u_{m}$
and
$h_{m}$
are yet to be determined) for
$y_{b}<x/t<y_{c}$
; and an expansion fan of
${\it\lambda}_{+}$
characteristics attached to the front of the flow
$y_{c}<x/t<y_{N}$
. The construction of the solutions proceeds as follows:
$y_{a}$
is given by (3.8) and within the
${\it\lambda}_{-}$
fan, the Riemann invariant is given by
$$\begin{eqnarray}R_{+}(u,h)=R_{+}(u_{0},h_{0}).\end{eqnarray}$$
Within the
${\it\lambda}_{+}$
fan of characteristics at the front of the current
$(y_{c}<x/t<y_{N})$
, the Riemann invariant
$R_{-}$
is constant and given by
$$\begin{eqnarray}R_{-}(u,h)=R_{-}(u_{c},h_{c}).\end{eqnarray}$$
Then in the uniform region between the two expansion fans, setting
$u=u_{m}$
and
$h=h_{m}$
, we deduce from (3.14) and (3.15) that
$$\begin{eqnarray}\displaystyle 1-2h_{m}=\sin \left({\textstyle \frac{1}{2}}(R_{-}(u_{c},h_{c})-R_{+}(u_{0},h_{0}))\right) & & \displaystyle\end{eqnarray}$$
and
$$\begin{eqnarray}\displaystyle u_{m}-\hat{Q}=(1-h_{m})\sin \left({\textstyle \frac{1}{2}}(R_{-}(u_{c},h_{c})+R_{+}(u_{0},h_{0}))\right). & & \displaystyle\end{eqnarray}$$
These expressions allow the transition locations,
$y_{b}$
and
$y_{c}$
to be determined by
$$\begin{eqnarray}y_{b}={\it\lambda}_{-}(u_{m},h_{m})\quad \text{and}\quad y_{c}={\it\lambda}_{+}(u_{m},h_{m}).\end{eqnarray}$$
Finally, profiles are determined by simultaneously solving (3.14) with
$$\begin{eqnarray}y={\it\lambda}_{-}(u,h)\quad \text{within }y_{a}<y<y_{b},\end{eqnarray}$$
and (3.15) with
$$\begin{eqnarray}y={\it\lambda}_{+}(u,h)\quad \text{within }y_{c}<y<y_{N},\end{eqnarray}$$
where
$y=x/t$
. Example solutions are plotted in figure 4 for
$F_{0}=1.2$
and
$F_{0}=1$
. Note that when the source is critical
$(F_{0}=1)$
, there is no uniform region before the expansion of
${\it\lambda}_{-}$
characteristics (i.e.
$y_{a}=0$
).
The flow fields in this regime are quite similar to those generated by the motion that follows the instantaneous removal of a lock gate separating fluids of differing densities. For such dambreak flows, provided the height of the reservoir is less than half of the height of the channel, the dense flowing layer produces a rarefaction wave that propagates into the initially undisturbed reservoir, in addition to regions within which the depth of the layer is constant and regions close to the front where the flow is choked (see appendix B).
3.3 Shocks
When the source Froude number is greater than the Froude number for which a uniform state would occur
$(F_{0}>F_{U}\cup F_{Uf})$
, the flow consists of a uniform state connected by a bore, which moves downstream at dimensionless speed
$V$
, to a state that may be uniform if the flow is not choked or may feature a rarefaction attached to the front if the flow is choked. The specification of the change of flow conditions across the bore requires a model of the mass and momentum fluxes beyond that embodied in the governing partial differential equations. Here we conserve the mass and momentum fluxes and assume that energy is conserved in the upper layer (Ungarish Reference Ungarish2009).
Upstream of the bore, the height of the dense fluid is
$h_{0}$
and the flow speed in the lower layer is
$u_{0}$
; the flow speed in the upper layer vanishes. Downstream of the bore we denote the dimensionless height of the layer by
$h_{1}$
and the dimensionless velocities in the lower and upper layers are
$u_{1}$
and
$u_{2}$
, respectively. Conservation of fluid mass in the upper layer yields
$$\begin{eqnarray}u_{2}=V\frac{(h_{0}-h_{1})}{(1-h_{1})},\end{eqnarray}$$
while conservation of mass in the lower layer yields
$$\begin{eqnarray}u_{1}=\frac{\hat{Q}}{h_{1}}+V-V\frac{h_{0}}{h_{1}}.\end{eqnarray}$$
The excess pressure fields, rendered dimensionless by
${\it\rho}_{2}gSH$
, are assumed to be hydrostatic on either side of the bore. Thus upstream of the bore
$$\begin{eqnarray}p_{L}=\left\{\begin{array}{@{}ll@{}}P_{L},\quad & h_{0}<z<1,\\ P_{L}+(h_{0}-z),\quad & 0<z<h_{0},\end{array}\right.\end{eqnarray}$$
and downstream of the bore
$$\begin{eqnarray}p_{R}=\left\{\begin{array}{@{}ll@{}}P_{R},\quad & h_{1}<z<1,\\ P_{R}+(h_{1}-z),\quad & 0<z<h_{1},\end{array}\right.\end{eqnarray}$$
where
$P_{R}$
and
$P_{L}$
are the reduced pressures at the surface
$z=1$
, downstream and upstream of the bore, respectively. Finally energy conservation in the upper layer is enforced by applying the Bernoulli equation along a streamline at the upper boundary and this gives
$$\begin{eqnarray}P_{L}-P_{R}={\textstyle \frac{1}{2}}((u_{2}-V)^{2}-V^{2}).\end{eqnarray}$$
Balancing the momentum flux across a control volume enclosing the bore and simplifying using (3.21) and (3.22) then yields the following condition for the bore speed,
$$\begin{eqnarray}V^{2}\frac{(h_{0}+h_{1}-2h_{0}h_{1})}{(1-h_{1})^{2}}+2(1+S)\left(\frac{\hat{Q}}{h_{0}}-V\right)^{2}\frac{h_{0}}{h_{1}}=h_{1}+h_{0},\end{eqnarray}$$
provided
$h_{1}\neq h_{0}$
. In the Boussinesq regime, (3.26) may be simplified by neglecting terms of order
$S$
. This approximation is made in this subsection, but in § 5, we analyse flows with relatively large density differences between the two fluids and in such cases, the full expression (3.26) must be used.
We may now couple this condition that connects the flow states on either side of a bore with conditions at the front to establish the height and velocity of the flowing layer. If the flow is not choked then the front condition is given by (2.22), which yields
$$\begin{eqnarray}(\hat{Q}(1-h_{1})+V(h_{1}-h_{0}))^{2}=h_{1}^{2}\mathscr{F}(h_{1}).\end{eqnarray}$$
Solving (3.26) and (3.27), we find
$V$
and
$h_{1}$
, given the dimensionless volume flux,
$\hat{Q}$
, and the source Froude number,
$F_{0}$
. Example solutions are plotted in figure 5. This construction remains valid provided the characteristic speed upstream of the shock is greater than the shock speed
$({\it\lambda}_{+}(u_{0},h_{0})>V)$
and that the flow height downstream of the bore is less than the critical height at which the flow becomes choked (i.e.
$h_{1}<h_{c}$
). In the
$(\hat{Q},F_{0})$
plane, the regime for which these solutions exist is bounded by the uniform flow condition,
$F_{0}>F_{U}$
and by a constraint deduced by substituting
$h_{1}=h_{c}$
and simultaneously solving (3.26) and (3.27); this latter condition represents a curve that we denote
$\mathscr{C}_{S}$
(see figure 2).
Figure 5. The height,
$h$
, and velocity,
$u$
, as functions of
$x/t$
for
$F_{0}=1.6$
and dimensionless flux
$\hat{Q}=0.1,0.2,0.2562,0.3,0.35,0.3906$
. Note that
$\hat{Q}\equiv \hat{Q}_{1}=0.2562$
corresponds to the maximum flux with source Froude number
$F_{0}=1.6$
for which the flow is not choked, while
$\hat{Q}=1.6^{-2}=0.3906$
is the maximum dimensionless flux with source Froude number
$F_{0}=1.6$
for which a two-layer flow attached to the source may be found. Increasing values of
$\hat{Q}$
produce successively deeper and faster currents.
It is also possible for the gravity current to feature a choked region at its front (see figure 2). In this case there is a rarefaction attached to the front through which the current thins to match the choked condition. This rarefaction is a fan of
${\it\lambda}_{+}$
characteristics and the Riemann invariant,
$R_{-}$
is constant throughout the region. Thus we have
$$\begin{eqnarray}R_{-}(u,h)=R_{-}(u_{c},h_{c}),\end{eqnarray}$$
and in particular, at the upstream edge of the rarefaction, where it meets with the solution downstream of the bore, we find
$$\begin{eqnarray}u_{1}-\hat{Q}=(1-h_{1})\sin (R_{-}(u_{c},h_{c})-\sin ^{-1}(1-2h_{1})).\end{eqnarray}$$
Thus matching with the bore conditions, we find
$$\begin{eqnarray}(\hat{Q}(1-h_{1})+V(h_{1}-h_{0}))^{2}=h_{1}^{2}(1-h_{1})^{2}\sin ^{2}(R_{-}(u_{c},h_{c})-\sin ^{-1}(1-2h_{1})).\end{eqnarray}$$
The shock speed,
$V$
, and interface height downstream of the bore are then found by solving (3.26) and (3.30). Finally the variation of height and velocity within the rarefaction are given by (3.20) and (3.28). Example solutions are plotted in figure 5.
3.4 Transition to a channel-filling or critical flow
When
$F_{0}^{2}\hat{Q}>1$
, the two-layer model of motion in this problem is no longer hyperbolic. Instead, as demonstrated by Long (Reference Long1955) and Boonkasame & Milewski (Reference Boonkasame and Milewski2011) amongst others, the flows are unstable. In this case, they expand rapidly and fill the channel in regions close to the source. Here we model this effect by treating the transition as instantaneous, localised to the source and occurring without mixing; these assumptions will be examined in the numerical solutions of § 4. There are then two possibilities for the configuration of the flow. Either there will be an expanding region attached to the source
$(0<x/t<y_{r})$
within which the dense fluid fills the channel, accompanied by a region further downstream from the source within which the two-layer current is re-established, or there is an immediate transition at the source to a two-layer flow state in which critical conditions are attained at the source.
First we construct the solutions for the former case in which there is a channel-filling transition attached to the source. At the rearward transition back to a two-layer flow, there is a jump in the height of the interface between the fluid from
$h=1$
to
$h=h_{r}$
. The motion close to this location is governed by a ‘front’ condition akin to that applied at the leading downstream edge of the gravity current (2.19), but here it is applied to the velocity field of the less dense fluid. This ‘front’ is choked since if it were not the fluid could adjust to find a state with a larger backward velocity relative to the mean flow in the channel. Thus
$h_{r}=1-h_{c}$
and the velocity of the less dense fluid is
$u_{t}=\hat{Q}-h_{r}^{3/2}$
. Immediately, therefore, we may deduce the criterion for when this type of flow exists. It requires the velocity of the less dense fluid to be positive and thus
$\hat{Q}>Q_{m}\equiv (1-h_{c})^{3/2}=0.5273$
. Using mass conservation
$(u_{r}h_{r}+u_{t}(1-h_{r})=\hat{Q})$
, we find that the velocity in the dense fluid at the jump is
$$\begin{eqnarray}u_{r}=\hat{Q}+h_{r}^{1/2}(1-h_{r}),\end{eqnarray}$$
and this jump occurs at
$x/t=y_{r}={\it\lambda}_{-}(u_{r},h_{r})=\hat{Q}-h_{r}^{3/2}$
.
The flow solution is then straightforward. For
$y_{r}<x/t<y_{b}$
, there is a rarefaction fan of
${\it\lambda}_{-}$
characteristics within which the Riemann invariant
$R_{+}$
is constant. Thus within this domain
$$\begin{eqnarray}R_{+}(u,h)=R_{+}(u_{r},h_{r})\equiv \sin ^{-1}(h_{c}^{1/2}(4h_{c}-3))\end{eqnarray}$$
and
$y=x/t={\it\lambda}_{-}(u,h)$
. Downstream of this rarefaction is a uniform region within which
$u=u_{m}$
and
$h=h_{m}$
for
$y_{b}<x/t<y_{c}$
and further downstream there is a rarefaction of
${\it\lambda}_{+}$
characteristics attached to the choked front of the flow,
$(y_{c}<x/t<y_{N})$
, such that
$$\begin{eqnarray}R_{-}(u,h)=R_{-}(u_{c},h_{c})\equiv \sin ^{-1}(h_{c}^{1/2}(4h_{c}-3)),\end{eqnarray}$$
and
$y_{N}=u_{c}$
. Thus we deduce that
$R_{+}(u_{m},h_{m})=R_{-}(u_{m},h_{m})$
and so
$h_{m}=1/2$
and
$u_{m}=\hat{Q}+h_{c}^{1/2}(3-4h_{c})/2$
. Also we find that
$$\begin{eqnarray}y_{b}=\hat{Q}-{\textstyle \frac{1}{2}}(1-h_{c})^{1/2}(4h_{c}-1)\quad \text{and}\quad y_{c}=\hat{Q}+{\textstyle \frac{1}{2}}(1-h_{c})^{1/2}(4h_{c}-1).\end{eqnarray}$$
Flow profiles are plotted in figure 6; we note that they are identical, merely being offset by the mean flow
$\hat{Q}$
.
When
$F_{0}^{-2}<\hat{Q}<Q_{m}$
, the flow can no longer transition to a state with the dense fluid filling the channel. Instead it transitions to a critical state at the origin
$({\it\lambda}_{-}=0)$
, with vanishing fluid velocity in the less dense layer. Thus at the origin we have
$$\begin{eqnarray}u_{0}=\hat{Q}^{1/3}\quad \text{and}\quad h_{0}=\hat{Q}^{2/3}.\end{eqnarray}$$
The flow now features a rarefaction fan of
${\it\lambda}_{-}$
characteristics, attached to the origin, a uniform region and potentially a fan of
${\it\lambda}_{+}$
characteristics if the flow is choked at the front. It may be constructed as follows. First for
$0<x/t<y_{b}$
,
$$\begin{eqnarray}R_{+}(u,h)=R_{+}(\hat{Q}^{1/3},\hat{Q}^{2/3}),\end{eqnarray}$$
and
$x/t={\it\lambda}_{-}(u,h)$
. For
$y_{b}<x/t<y_{c}$
, there is a uniform region
$u=u_{m}$
,
$h=h_{m}$
and
$$\begin{eqnarray}y_{b}={\it\lambda}_{-}(u_{m},h_{m})\quad \text{and}\quad y_{c}={\it\lambda}_{+}(u_{m},h_{m}).\end{eqnarray}$$
Finally if the flow is choked at the front
$(h_{m}>h_{c})$
then within
$y_{c}<x/t<y_{N}=\hat{Q}+(1-h_{c})^{3/2}$
,
$$\begin{eqnarray}R_{-}(u,h)=R_{-}(u_{c},h_{c}),\end{eqnarray}$$
and
$x/t={\it\lambda}_{+}(u,h)$
. The limiting condition for choking at the front occurs when the flow adjusts to
$h=h_{c}$
and
$u=u_{c}$
within the rarefaction attached to the source. It is given by
$R_{+}(u_{c},h_{c})=R_{+}(\hat{Q}^{1/2},\hat{Q}^{2/3})$
and this occurs as
$\hat{Q}=Q_{mc}\equiv 0.3001$
. Then for
$\hat{Q}_{m}>\hat{Q}>Q_{mc}$
the flow is choked and there is a rarefaction fan of characteristics attached to the front, whereas for
$\hat{Q}<Q_{mc}$
, the flow is not choked and there is no frontal fan. Some typical profiles are plotted in figure 6.
Figure 6. The height,
$h$
, and velocity,
$u$
as functions of
$x/t$
for two-layer currents that have lost their hyperbolic structure and transition to either channel-filling flows or to flows with critical control at the origin. In this plot,
$\hat{Q}=0.2,Q_{mc}\equiv 0.3001,0.4,Q_{m}=0.5273,1.1$
.
4 Simulations
We now proceed to analyse the flow by means of two-dimensional direct numerical simulation (DNS) in the Boussinesq regime, based on our code TURBINS (Nasr-Azadani & Meiburg Reference Nasr-Azadani and Meiburg2011), which employs a finite-difference discretization combined with a fractional projection method and TVD-RK3 time integration. Numerical details and validation results of the code are provided in Nasr-Azadani & Meiburg (Reference Nasr-Azadani and Meiburg2011) and Nasr-Azadani, Hall & Meiburg (Reference Nasr-Azadani, Hall and Meiburg2013), so that we provide only a brief review here.
The fluid motion is modelled using the dimensionless incompressible Navier–Stokes equations in the Boussinesq regime, given by
$$\begin{eqnarray}\displaystyle & \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial \boldsymbol{u}}{\partial t}+\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{u}=-\boldsymbol{{\rm\nabla}}p+\frac{1}{Re}{\rm\nabla}^{2}\boldsymbol{u}-\hat{{\it\rho}}\hat{\boldsymbol{z}}, & \displaystyle\end{eqnarray}$$
where
$\boldsymbol{u}$
,
$p$
and
$\hat{{\it\rho}}$
are the dimensionless, two-dimensional velocity field, pressure and excess density
$(0\leqslant \hat{{\it\rho}}\leqslant 1)$
, respectively, and where
$\hat{\boldsymbol{z}}$
and
$Re$
denote a unit vector aligned with the direction of gravity and the Reynolds number, respectively. For the purposes of the numerical simulations, it is convenient to non-dimensionalise the variables differently from the shallow-layer analysis of § 2: in the governing equations (4.1) and (4.2), we scale all lengths with respect to the height of the inflow,
$h_{0}^{\ast }\equiv (Q^{2}/(gSF_{0}^{2}))^{1/3}$
, and time with respect to
$t_{0}^{\ast }\equiv (h_{0}^{\ast }/(gS))^{1/2}$
. The Reynolds number thus is obtained as
$$\begin{eqnarray}\displaystyle Re=\frac{h_{0}^{\ast 2}}{t_{0}^{\ast }{\it\nu}}, & & \displaystyle\end{eqnarray}$$
and the dimensionless height of the channel by
${\hat{H}}=H/h_{0}^{\ast }$
. Under these scalings, and keeping the height of the source constant, the rescaled inflow velocity is
$F_{0}$
and variations in the dimensionless flux parameter,
$\hat{Q}$
, are obtained by altering the height of the computational channel, because
$\hat{Q}=F_{0}{\hat{H}}^{-3/2}$
. Although the computations are performed with these new dimensionless variables, we analyse the motion quantitatively in terms of the variables defined in § 2, to enable immediate comparison between the shallow-layer theory and the computations.
To model the evolution of the dimensionless excess density field
$\hat{{\it\rho}}(\boldsymbol{x},t)$
, we employ a continuum description and evolve the density field in an Eulerian manner via
$$\begin{eqnarray}\displaystyle \frac{\partial \hat{{\it\rho}}}{\partial t}+\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\hat{{\it\rho}}=\frac{1}{ReSc}{\rm\nabla}^{2}\hat{{\it\rho}}. & & \displaystyle\end{eqnarray}$$
Here,
$Sc$
denotes the Schmidt number associated with the diffusion of the density field
$\hat{{\it\rho}}$
. In most applications
$Sc\gg 1$
, although test simulations indicate that the precise value of
$Sc$
influences the results only weakly as long as
$Sc\geqslant O(1)$
(Härtel et al.
Reference Härtel, Meiburg and Necker2000). With this in mind, we employ
$Sc=14$
in all of our simulations. We remark that in most simulations
$Re$
is kept at 350 to ensure a stable interface between the dense current and the ambient fluid.
The computational set-up employs a channel of dimensionless size
$L_{x}\times L_{z}=50\times {\hat{H}}$
, with in/outflow boundaries in the horizontal direction, cf. figure 1. The Cartesian grid in the
$x$
- and
$z$
-directions is uniformly spaced with
${\rm\Delta}x=0.1$
and
${\rm\Delta}z=0.016$
. We apply free-slip conditions at the top and bottom walls. At the outflow boundary, all flow variables,
$q$
, are convected out of the domain via the outflow boundary condition
$$\begin{eqnarray}\displaystyle \frac{\partial q}{\partial t}+\bar{U}\frac{\partial q}{\partial x}=0, & & \displaystyle\end{eqnarray}$$
where
$\bar{U}$
represents the maximum
$u$
velocity value in the domain. For the density field, we impose no-flux conditions at the top and bottom walls. At the inlet and for
$z<1$
, the density and horizontal velocity are set to unity and
$F_{0}$
, respectively. These values decay smoothly over a few grid points to a zero value for
$z>1$
. The computational domain is initially filled with ambient fluid at rest.
In figure 7 we plot contours of the concentration field at dimensionless times
$t=20{\hat{H}}^{1/2}$
for varying values of the flux parameters
$\hat{Q}$
, while keeping the source Froude number constant
$(F_{0}=1.2)$
. To reiterate, here, and in the subsequent discussion, the variables are non-dimensionalised as described in § 2. It is noteworthy that the currents maintain a relatively sharp interface between the two fluids in the channel; turbulence and other diffusive processes are insufficient to cause significant mixing on these length and time scales. The other noticeable feature is that for the higher values of the dimensionless flux,
$\hat{Q}$
, the dense current has completely filled the channel depth in an expanding region close to the source, while further downstream the current re-establishes the two-layer structure.
Figure 7. Contours of the dimensionless excess density field from the DNS at
$t=20{\hat{H}}^{1/2}$
for source Froude number,
$F_{0}=1.2$
, and various values of the dimensionless flux,
$\hat{Q}=0.15,0.28,0.42,0.70,0.78,0.94$
. The ten contours in each panel correspond to values of
$\hat{{\it\rho}}$
from 0.1 to 1 with increments of 0.1.
We calculate the rate at which the current progresses downstream by defining a scaled current height
$$\begin{eqnarray}{\it\Phi}=\frac{1}{{\hat{H}}}\int _{0}^{{\hat{H}}}\hat{{\it\rho}}\,\text{d}z.\end{eqnarray}$$
based on the vertical integral of the density field. The front of the current,
$x_{N}(t)$
is then determined by finding the maximum value of
$x$
for which
${\it\Phi}>{\it\epsilon}$
, where
${\it\epsilon}$
is a prescribed threshold, typically
${\it\epsilon}=0.01$
. We plot
$x_{N}(t)$
as a function of time for each computational run in figure 8(a). We observe that the currents move at constant speeds dependent on
$\hat{Q}$
. By fitting a linear trend line to each set of data we measure this speed and compare with the prediction of the shallow-layer model in figure 8(b). The agreement between the model and the simulations is quite good. The Navier–Stokes currents are seen to move slightly more slowly than predicted by the shallow-layer model, which is most likely due to their relatively low Reynolds number of 350, implying the weak, but non-vanishing, action of viscous processes that retard the motion.
We examined the effects of varying the Reynolds number in two ways. First for
$F_{0}=1.2$
and
$\hat{Q}=0.42$
, we performed simulations with Reynolds numbers ranging from 100 to 1500 and we measured the speed of the flow. Computing at higher Reynolds numbers required greater spatial resolution to achieve converged results and the computational grids were adjusted appropriately. The results are given in table 1, where we observe that the measured speed systematically increases with increasing Reynolds number, approaching the shallow-layer value. This supports the view that the slight mismatch between the simulations and the inviscid shallow-layer prediction is due to the weak, but non-vanishing, effects of viscosity.
Figure 8. (a) The position of the front of the current,
$x_{N}(t)$
, as a function of time for each of the simulations with
$F_{0}=1.2$
; and (b) the measured velocity of the current,
$u_{N}$
, as a function of the dimensionless flux,
$\hat{Q}$
at
$Re=350$
(○) and
$Re=1500$
(
$\times$
). Also plotted are the predictions from the shallow-layer analysis (solid line) and the asymptotic result when
$\hat{Q}\ll 1$
(dotted line) and
$\hat{Q}=0.2938$
at which value the flow becomes choked.
Table 1. The dimensionless speed of the front of the current,
$u_{N}$
, measured from simulations conducted at different Reynolds numbers, when the source Froude number,
$F_{0}=1.2$
and the dimensionless flux per unit width,
$\hat{Q}=0.42$
. The shallow-layer prediction for these conditions is
$u_{N}=0.952$
.
Contours of the concentration field from each run at different Reynolds number are plotted in figure 9 at the same instant of time. The increased size and intensity of vortices on the interface between the two fluids becomes more evident as the Reynolds number is increased and these vortices lead to zones within which the density takes intermediate values between the intruding and ambient fluids. However over the time and length scales of these simulations, the emergence of this mixed zone does not significantly modify the speed of the current, presumably because the vertical integral of the excess density field over the depth of the current,
$(\int _{0}^{h}\hat{{\it\rho}}\,\text{d}z)$
, remains approximately constant under mixing and this determines the speed of the density-driven spreading. Johnson & Hogg (Reference Johnson and Hogg2013) have demonstrated that entrainment alone modifies the predictions of shallow-layer models, but only at time and length scales far beyond those computed in these simulations.
We also examine the effects of Reynolds number by performing simulations at
$Re=1500$
for source Froude number
$F_{0}=1.2$
and a range of values of the dimensionless flux,
$\hat{Q}$
. In each case, as for the simulations run with
$Re=350$
, we find the front of the current moves with constant speed and that this speed compares favourably with the shallow-layer predictions (see figure 8
b). For each case, the simulations at higher Reynolds number led to slightly higher front speeds, evidencing the reduced effects of viscosity (table 1).
The shallow-layer predictions always comprise currents with rarefactions, with the transition to choked flows occurring at
$\hat{Q}=0.2938$
and the loss of hyperbolicity, with the associated transition to channel-filling flows, at
$\hat{Q}=0.6944$
. This latter value closely agrees with the value observed in the Navier–Stokes simulations. Also plotted is the asymptotic result for a relatively deep ambient
$(\hat{Q}\ll 1)$
. In this regime, to leading order, the Riemann invariants are given by
$R_{\pm }=u\pm 2\sqrt{h}$
on characteristics
$\text{d}x/\text{d}t=u\pm \sqrt{h}$
. For this source condition there is an expansion fan within which
$R_{+}=\hat{Q}^{1/3}(F_{0}^{2/3}+2F_{0}^{-1/3})$
. The leading edge of the expansion fan corresponds to
$u_{N}=\sqrt{2h_{N}}$
and thus
$$\begin{eqnarray}u_{N}=\frac{\hat{Q}^{1/3}\left(F_{0}^{2/3}+2F_{0}^{-1/3}\right)}{1+\sqrt{2}},\end{eqnarray}$$
a result which is quite accurate for
$\hat{Q}\leqslant 0.3$
when
$F_{0}=1.2$
(see figure 8). This implies that the gravity current motion is relatively unaffected by the motion of the overlying less dense fluid until the flow becomes choked at
$\hat{Q}=0.2938$
. Notably when
$\hat{Q}\ll 1$
, the dimensionless rate of propagation is proportional to
$\hat{Q}^{1/3}$
, whereas for larger values, the speed varies linearly with
$\hat{Q}$
and from (2.20) is given by
$u_{N}=0.5473+\hat{Q}$
.
Figure 9. Contours of the dimensionless excess density field from the DNS at
$t=20{\hat{H}}^{1/2}$
for source Froude number,
$F_{0}=1.2$
, and dimensionless flux,
$\hat{Q}=0.42$
, at various Reynolds numbers. The ten contours in each panel correspond to values of
$\hat{{\it\rho}}$
from 0.1 to 1 with increments of 0.1.
We may also investigate the motion of the interface between the dense and less dense fluids at the upper boundary when the flow is channel filling. In this case, we define the position
$x_{r}(t)$
as the maximum value of
$x$
for which
${\it\Phi}>1-{\it\epsilon}_{1}$
, where
${\it\epsilon}_{1}$
is another prescribed threshold with chosen value
${\it\epsilon}_{1}=0.05$
. We note from figure 7 that for those flows which undergo the transition to a channel-filling motion, although the interface at the top boundary is relatively sharp (i.e. the region of intermediate density is quite narrow), the
$\hat{{\it\rho}}=1$
contour is somewhat displaced from the upper boundary. This means that the threshold
${\it\epsilon}_{1}$
cannot be chosen to be arbitrarily small or else the criterion for determining the position
$x_{r}(t)$
would fail. However, we find that
$x_{r}(t)$
depends only very weakly on
${\it\epsilon}_{1}$
for
$0.04<{\it\epsilon}_{1}<0.1$
. In figure 10 we plot
$x_{r}(t)$
as a function of time for those runs for which the flow becomes channel filling and we overlay the theoretical prediction. It is evident that there is some temporal offset before the channel-filling layer becomes established, presumably due to the development of the interfacial instability. Thereafter, the theoretical model captures the measured speed quite well until the basal front of the two-layer current leaves the domain. Beyond that time, the existence of the outflow boundary prevents the flow in the domain from continuing to be influenced by the motion of the current front. Longer computational domains are expected to produce agreement between the two sets of data for longer times. We examine the motion of this interface at the top of the channel for simulations at higher Reynolds number (
$Re=1500$
, see figure 10). We observe that at higher Reynolds number, the channel-filling part of the flow takes longer to be fully established and that for the lowest flux considered,
$\hat{Q}=0.70$
, it is not completely initiated before the front of the flow leaves the computational domain. For the other cases in which this channel-filling transition occurs, however, the measured positions are broadly in accord with the theoretical predictions after some initial period during which the flow structure is developed.
Figure 10. The position of the interface between the dense and less fluids at the upper boundary of the channel,
$x_{r}(t)$
as a function of time for cases in which the gravity current becomes channel filling close to the source from the DNS results at
$Re=350$
(solid lines), at
$Re=1500$
(dotted lines) and the shallow-layer modelling (dashed lines). The DNS results are curtailed when the basal front leaves the fluid domain.
Gravity currents due to the sustained emission of relatively dense fluid within a channel have recently been studied by Shringapure et al. (Reference Shringapure, Lee, Ungarish and Balachandar2013). These authors conducted numerical simulations of the Navier–Stokes equations and compared their results with a simple analytical model of the motion in which the currents are spatially uniform. Here we quantitatively investigate their results in terms of the shallow-layer formulation developed above (§ 3). We note that the source conditions they treated were subcritical
$({\it\lambda}_{-}(u_{0},h_{0})<0)$
. Thus the flows are subject to downstream influences and adjust close to source to critical conditions
$({\it\lambda}_{-}(u_{0},h_{0})=0)$
, which implies
$F_{0}=1$
. We may then compare the numerically determined speed of propagation and average height over a small frontal region with the shallow-layer model developed here. We plot in figure 11 the dimensionless front position and height at the front as a function of the dimensionless flux, noting that there is very good agreement between the two. This provides further confidence in the shallow-layer model since it is capable of accurately reproducing features in an independently developed numerical simulation.
Figure 11. The dimensionless speed,
$u_{N}$
, and height,
$h_{N}$
, of the front as a function of the dimensionless flux per unit width,
$\hat{Q}$
calculated from the shallow layer with
$F_{0}=1$
(solid lines) and from the numerical simulations of Shringapure et al. (Reference Shringapure, Lee, Ungarish and Balachandar2013) (
$\times$
). Also marked is the value of the flux,
$\hat{Q}_{c}$
at which the motion first becomes choked (dashed line).
5 Non-Boussinesq flows
$(S=O(1))$
For non-Boussinesq currents, the same types of flow structures arise, namely, uniform and choked flows with shocks and rarefactions. However the boundaries between the different states in the
$(\hat{Q},F_{0})$
plane vary with
$S$
and our method of constructing them differs from the Boussinesq regime. Instead of analytical expressions for the Riemann invariants, we must now largely rely on numerically computed solutions to the differential equation along each family of characteristics (2.11).
Non-Boussinesq currents with no spatial variations between the source and the front may be determined through the solutions of coupled algebraic equations; they satisfy (3.1) for
$0\leqslant \hat{Q}\leqslant Q_{c}$
. This produces a current identical to the Boussinesq case (as plotted in figure 3). For larger values of the dimensionless flux
$(\hat{Q}>\hat{Q}_{c})$
, the flow becomes choked at the front and there is a rarefaction of
${\it\lambda}_{+}$
characteristics as the dense fluid accelerates to match the front condition. The curve in the
$(\hat{Q},F_{0})$
plane is mapped out by integrating
$$\begin{eqnarray}\frac{\text{d}u}{\text{d}h}=-\frac{{\it\alpha}^{2}+\sqrt{{\it\alpha}^{2}+(1-{\it\beta})h}}{h},\end{eqnarray}$$
subject to
$u=u_{c}$
when
$h=h_{c}$
and
$u=u_{0}$
when
$h=h_{0}$
, thus providing an implicit relationship between the parameters
$\hat{Q}$
and
$F_{0}$
. The form of this curve is dependent upon the dimensionless excess density,
$S$
, as illustrated in figure 12 for
$S=1$
and
$S=5$
. However it is evident that there always exist some parameter values for which there is ‘choked’ flow because
$\hat{Q}_{c}<F_{0c}(F_{0c}^{2}+S/(1+S))^{-3/2}$
, the condition for hyperbolicity (2.24) evaluated at
$\hat{Q}=\hat{Q}_{c}$
and
$F_{0}=F_{0c}$
. As in § 3.1, we denote the curves on which there is uniform and choked-uniform flow by
$F_{U}$
and
$F_{Uf}$
, respectively.
Figure 12. Regimes of solutions, classified as function of the dimensionless flux,
$\hat{Q}$
and the source Froude number,
$F_{0}$
for non-Boussinesq flows (a)
$S=1$
and (b)
$S=5$
. As in figure 2, uniform flows correspond to
$F_{0}=F_{U}\cup F_{Uf}$
. Smaller, but still supercritical, values of
$F_{0}$
yield rarefactions, while larger values lead to shocks. The curve corresponding to the loss of hyperbolicity is denoted
$\mathscr{C}_{h}$
, while the transition between choked and unchoked flows of rarefactions and shocks are denoted
$\mathscr{C}_{R}$
and
$\mathscr{C}_{S}$
, respectively. Critical conditions are given by
$F_{0}=1/\sqrt{2}$
for
$S=1$
and
$F_{0}=1/\sqrt{6}$
for
$S=5$
; these values are also plotted.
For
$F_{0}<F_{U}\cup F_{Uf}$
, there are rarefactions close to, or attached to, the source and potentially close to the front if the motion is ‘choked’. We may compute the boundary between ‘choked’ and unchoked rarefaction by integration across the rarefaction of
${\it\lambda}_{-}$
characteristics close to source. Thus we integrate
$$\begin{eqnarray}\frac{\text{d}u}{\text{d}h}=-\frac{{\it\alpha}^{2}-\sqrt{{\it\alpha}^{2}+(1-{\it\beta})h}}{h},\end{eqnarray}$$
from
$u=u_{0}$
and
$h=h_{0}$
to
$u=u_{c}$
and
$h=h_{c}$
and this provides an implicit relationship between the dimensionless source flux and source Froude number. The form of this curve,
$\mathscr{C}_{R}$
, is dependent on
$S$
as shown in figure 12.
When
$F_{0}>F_{U}\cup F_{Uf}$
, the flow involves shocks across which the height and velocity vary discontinuously. As with the Boussinesq currents (§ 3), the jump conditions feature the conservation of mass and momentum fluxes, together with the assumption that there is no dissipation in the upper, less dense fluid. Then the currents may either comprise two uniform segments, joined by a shock, or may be choked at the front, in which case the current thins and accelerates downstream of the shock. The boundary between these two possibilities is similar to that deduced above (§ 3.3), although the momentum equations on either side of the shock are modified to take account of the potentially non-negligible difference in densities between the two fluids (see § 2). We may then compute which type of current exists for different values of the dimensionless flux,
$\hat{Q}$
, and the source Froude number,
$F_{0}$
(see figure 12). We note that for the extent of values for which choked flow exists diminishes as the density ratio increases.
In figure 13(a,b), we plot the profiles of height and velocity for non-Boussinesq currents for different ratios of the densities of the fluids. In particular, we evaluate the form of the currents when the source is critical,
$F_{0}=(1+S)^{-1/2}$
and the flux
$\hat{Q}={\it\mu}(1+S)^{-1/2}$
with
$0\leqslant {\it\mu}\leqslant 1$
. (In this range of values the source is always hyperbolic.) We observe that for all of the flows the current immediately transitions from their source conditions through a rarefaction to a uniform state. Furthermore, for the cases plotted, the flow is only choked for relatively large dimensionless fluxes with
$S=1$
. As the density ratio becomes large
$(S\gg 1)$
, the motion of the less dense fluid in the upper layer plays only a negligible role and the motion is modelled by the single-layer shallow-water equation to leading order. In this regime, we must adopt different variables, namely
$\overline{t}=(1+S)^{-1/2}t$
,
$\overline{u}=(1=S)^{1/2}u$
and
$\overline{Q}=(1+S)^{1/2}\hat{Q}$
. The front condition is then given by
$h=0$
to leading order and the velocity and height fields satisfy
$$\begin{eqnarray}\overline{u}=\frac{1}{3}\left(\overline{Q}+2+\frac{2x}{\overline{t}}\right)\quad \text{and}\quad h=\frac{1}{9}\left(\overline{Q}+2-\frac{x}{\overline{t}}\right)^{2}.\end{eqnarray}$$
In figure 13(c) we plot the height and velocity fields for
$\hat{Q}=F_{0}=(1+S)^{-1/2}$
for a range of values of
$S$
in terms of these rescaled variables and we note that as
$S$
increases, the single-layer results are progressively attained.
Figure 13. The height,
$h$
, and velocity,
$u$
as function of
$x/t$
for non-Boussinesq, supercritical source conditions
$F_{0}=(1+S)^{-1/2}$
; (a)
$S=1$
and (b)
$S=5$
with
$\hat{Q}=n(1+S)^{-1/2}/5$
$(n=1,2,3,4,5)$
. Also plotted in (c) are the height and rescaled velocity fields as a functions of
$(S+1)^{1/2}x/t$
for
$F_{0}=\hat{Q}=(1+S)^{-1/2}$
.
6 Conclusions
We have analysed the unsteady motion that results from the sustained release of dense fluid into a horizontal channel which initially contains less dense fluid. We have demonstrated how the ensuing inviscid motion may be analysed using a two-layer shallow-water model and that there are three independent dimensionless groups that characterise the motion, which here we have formulated as the dimensionless flux per unit width,
$\hat{Q}$
, the source Froude number,
$F_{0}$
and the relative excess density,
$S$
. Our results, derived from the solutions of simple algebraic equations in the Boussinesq regime and the integration of a single first-order differential equation in the non-Boussinesq regime, have demonstrated that various types of flows may arise; these include uniform flows, those with internal discontinuities (shocks) and those that vary continuously (rarefactions). Additionally the motion may become choked when the speed of advance is limited by the rate at which dense fluid can be supplied to the front. The identification and demarcation of these variations, including their dependence upon the relative density,
$S$
, adds significantly to what has been reported before and reveals a rich interplay of effects. We also find solutions in which the two-layer structure with less dense fluid overlying dense fluid cannot be sustained. In these situations, the dense fluid undergoes a transition at source to fill the entire depth of the channel before subsequently re-establishing the two-layer form further downstream.
We have reported results from the direct numerical simulation of the full Navier–Stokes equations and have demonstrated very good agreement in the prediction of bulk properties between the idealised shallow models and these more complete descriptions. The simulations may be probed to reveal the spatial and temporal variations of the velocity, pressure and density fields and thus provide much greater information about the unsteady flow associated with any particular source conditions. However they require significant computational resource and so it is costly to determine the parametric dependence that is brought out clearly by the shallow-layer models. Therefore it is the interplay between the two approaches that has provided insight into the underlying dynamics and confidence in the modelling approaches.
One important component of shallow-layer modelling is the dynamic condition that captures the non-hydrostatic motion at the front (Benjamin Reference Benjamin1968; Borden & Meiburg Reference Borden and Meiburg2013a ). Although our shallow-water analysis requires the specification of this condition, in the Boussinesq regime, there is relatively little difference between the form of the solutions or the domains in which they exist. Both yield quantitative results close to the numerical simulations. It would be interesting to compare the non-Boussinesq predictions produced with the Benjamin front condition with other formulations and with numerical simulations when they become available.
Our analytical techniques have also been applied to the related problem of Boussinesq dambreak flows (appendix B). Flows driven by sustained sources and by the collapse of a dam share many features as drawn out for single-layer flows by Gratton & Vigo (Reference Gratton and Vigo1994). Here, we have employed the method of characteristics to construct the height and velocity fields for both partial and full-depth dams; these solutions are of self-similar form and like those due to a sustained source, yield a simple (linear) gearing between the spatial and temporal variables. Our analytical solutions therefore add to the growing literature on dambreak solutions, which as well as being of interest in their own right, provide a test-bed solution for examining non-trivial unsteady flows.
Although we have not pursued this here, we comment that our results have direct application to industrial and environmental problems. In addition to the prediction of flow speeds, our results have shown that, for certain source conditions, the dense fluid completely displaces the less dense material, and that in other situations, the speed of propagation may become choked, increasing linearly with further increases of the source flux. Our shallow-layer formulation is straightforward to extend to non-rectangular cross-sections or release conditions that do not demand outflow at the end of the channel. It would also be interesting to examine how mixing between the fluids modifies the motion. We finally remark that there is a pressing need for further laboratory experiments to investigate the motion in both Boussinesq and non-Boussinesq regimes.
Acknowledgements
A.J.H. acknowledges the financial support of EPSRC (EP/G066353/1). M.M.N.A. and E.M. were funded by the National Science Foundation under grant no. CBET-1335148. This research was also supported in part by the National Science Foundation under grant no. PHY11-25915.
Appendix A
A.1 Circulation-based models of bores
In this section of the appendix, we compute the solutions for Boussinesq gravity current motion due to a sustained influx using the dynamic front and bore conditions for Boussinesq flows proposed by Borden & Meiburg (Reference Borden and Meiburg2013a
) and Borden & Meiburg (Reference Borden and Meiburg2013b
), respectively. These are based upon the conservation of mass and circulation across the front, the latter balancing the vorticity generated at the interface between the fluids with the baroclinic torque due to the excess density (see Borden & Meiburg Reference Borden and Meiburg2013a
). Here we identify the regimes of solutions types in terms of the governing parameters,
$\hat{Q}$
and
$F_{0}$
. The solutions developed here are therefore counterparts to those of § 2 using the Benjamin front condition (2.19).
Adapting the approach of Borden & Meiburg (Reference Borden and Meiburg2013a ) to account for a constant flux along the channel, we find that the front condition is given by
$$\begin{eqnarray}(u_{N}-\hat{Q})^{2}=2h_{N}(1-h_{N})^{2}\end{eqnarray}$$
this formulation differs from the one proposed by Benjamin (Reference Benjamin1968) and there are some immediate consequences. First, choked flow arises when
${\it\lambda}_{+}(u_{N},h_{N})=u_{N}$
and when combined with (A 1), this corresponds to
$h_{N}=h_{c}=1/3$
and
$u=u_{N}=\hat{Q}+0.5443$
. This corresponds to flows that are marginally shallower and faster than the choked flow conditions predicted by the Benjamin front condition ((2.19) and (2.20)).
Uniform flow states occur when
$h_{N}=h_{0}=(\hat{Q}/F_{0})^{2/3}$
and
$u_{N}=u_{0}=(\hat{Q}F_{0}^{2})^{1/3}$
and thus we find that they are given by
$F_{0}=\sqrt{2}$
, with the proviso that the motion is not choked, a condition given by
$\hat{Q}<\hat{Q}_{c}\equiv (2/27)^{1/2}$
. Flows that are choked at the front and uniform up to that point feature a fan of
${\it\lambda}_{+}$
characteristics, within which the Riemann invariant,
$R_{-}$
is constant. The curve in the
$(\hat{Q},F_{0})$
plane corresponding to this type of solution is given by
$$\begin{eqnarray}R_{-}(u_{0},h_{0})=R_{-}(u_{c},h_{c}).\end{eqnarray}$$
As in § 3, we denote the curves along which uniform and choked-uniform solutions exist by
$F_{U}$
and
$F_{Uf}$
, respectively. It is plotted in figure 14.
For
$F_{0}<F_{U}\cup F_{Uf}$
, solutions with rarefaction exist. We compute the onset of choked solution by evaluating the parameter values for which the flow evolves from source through a rarefaction of
${\it\lambda}_{-}$
characteristics to the choked flow state. This corresponds to a curve given implicitly by
$$\begin{eqnarray}R_{+}(u_{0},h_{0})=R_{+}(u_{c},h_{c}).\end{eqnarray}$$
This curve is denoted
$\mathscr{C}_{R}$
and is also plotted in figure 14.
Finally we analyse those solutions with internal bores. To this end, adapting the formulation developed by Borden & Meiburg (Reference Borden and Meiburg2013b ), we enforce mass conservation in both layers (3.21) and (3.22) and then balance the net flux of vorticity with the baroclinic torque generated across the bore. This yields
$$\begin{eqnarray}((u_{2}-V)^{2}-(u_{1}-V)^{2})-(V^{2}-(u_{0}-V)^{2})=2(h_{1}-h_{0}).\end{eqnarray}$$
The onset of choked flows for the parameter values with bores can then be deduced by setting
$h_{1}=h_{c}$
and simultaneously solving (3.21), (3.22), (A 4) and (3.27). This leads to the curve
$\mathscr{C}_{S}$
in the
$(\hat{Q},F_{0})$
plane, which is plotted in figure 14.
Figure 14. Regimes of Boussinesq solutions, classified in terms of the dimensionless flux per unit width,
$\hat{Q}$
, and the source Froude number,
$F_{0}$
, when circulation-based models are used for the conditions at the front and across bores.
A.2 Empirical front conditions
We now examine gravity current motion in the Boussinesq regime as a function of the source Froude number,
$F_{0}$
, and the dimensionless flux,
$\hat{Q}$
, when an empirical model is employed to specify the dynamical condition at the front of the current. Here we examine the consequences of using the front condition proposed by Huppert & Simpson (Reference Huppert and Simpson1980), modified to account for the sustained flux in the channel; in terms of the dimensionless variables defined in this study, this condition is given by
$$\begin{eqnarray}(u_{N}-\hat{Q})^{2}=\left\{\begin{array}{@{}ll@{}}{\textstyle \frac{1}{4}}{\hat{h}}^{-2/3}h_{N},\quad & 0<h_{N}<{\hat{h}},\\ {\textstyle \frac{1}{4}}h_{N}^{1/3},\quad & {\hat{h}}<h_{N},\end{array}\right.\end{eqnarray}$$
where
${\hat{h}}=0.075$
.
Figure 15. Regimes of Boussinesq solutions, classified in terms of the dimensionless flux per unit width,
$\hat{Q}$
, and the source Froude number,
$F_{0}$
, when an empirical model (Huppert & Simpson Reference Huppert and Simpson1980) is used for the front condition.
We may now construct the forms of solution, replacing the Benjamin front condition (2.19) with (A 5). First, we note that choked flow arises when
${\it\lambda}_{=}(u_{N},h_{N})=u_{N}$
, which occurs when
$(1-h_{N})^{3}=h_{N}^{1/3}/4$
. Thus we find that this occurs when
$h_{N}=h_{c}=0.4269$
and
$u_{N}=\hat{Q}+0.4339$
. This implies that choked flows are both deeper and slower than the counterparts predict using the frontal conditions due to Benjamin (Reference Benjamin1968) and Borden & Meiburg (Reference Borden and Meiburg2013a
).
Uniform flow states may be found when the velocity and depth of the current do not vary along its length. They are given by
$h_{N}=h_{0}=(\hat{Q}/F_{0})^{2/3}$
and
$u_{N}=u_{0}=(\hat{Q}F_{0}^{2})^{1/3}$
. Thus from (A 5), we deduce
$$\begin{eqnarray}h_{N}=1-\frac{{\hat{h}}^{-1/3}}{2F_{0}},\quad \text{for }0<h_{N}<{\hat{h}},\end{eqnarray}$$
while
$$\begin{eqnarray}\hat{Q}^{2}=\frac{h_{N}^{7/3}}{4(1-h_{N})^{2}},\quad \text{for }{\hat{h}}<h_{N}<h_{c}.\end{eqnarray}$$
Together, these defined a curve in the
$(\hat{Q},F_{0})$
plane along which a uniform state may be found (see figure 15). We denote this curve
$\mathscr{F}_{U}$
and it is defined for
$0<\hat{Q}<\hat{Q}_{c}\equiv 0.3232$
; note that
$\mathscr{F}_{U}(\hat{Q}_{c})=F_{0c}\equiv 1.1587$
and this is point is plotted in the figure. For larger values of the dimensionless flux, the currents are choked and feature a uniform portion attached to the source, joined to the front by a rarefaction of
${\it\lambda}_{+}$
characteristics within which the flow thins and accelerates. The curve along which currents of this type may be found is constructed by using the constancy of the
$R_{-}$
Riemann invariant. Thus
$R_{-}(u_{0},h_{0})=R_{-}(u_{c},h_{c})$
and this defines the curve
$\mathscr{F}_{Uf}$
in the
$(\hat{Q},F_{0})$
plane (see figure 15). This type of the current ceases to exist when the source conditions are no longer compatible with a two-layer flow. This loss of hyperbolicity occurs when
$F_{0}^{2}\hat{Q}=1$
and occurs when
$\hat{Q}=\hat{Q}_{h}\equiv 0.7678$
(see (3.6)).
When the parameters do not lie on the curve
$\mathscr{F}_{U}\cup \mathscr{F}_{Uf}$
, there are rarefactions and shocks. Rarefactions occur when
$F_{0}<\mathscr{F}_{U}\cup \mathscr{F}_{Uf}$
; through the rarefaction close to the source, the currents accelerate and thin. These flows may become choked and we may compute the transition to choked rarefaction by using the constancy of the
$R_{+}$
Riemann invariant through the fan of
${\it\lambda}_{-}$
characteristics. Thus
$R_{+}(u_{0},h_{0})=R_{+}(u_{c},h_{c})$
and this defines the curve
$\mathscr{C}_{R}$
in the
$(\hat{Q},F_{0})$
plane (figure 15).
Finally, the currents may feature shocks if
$F_{0}>\mathscr{F}_{U}\cup \mathscr{F}_{Uf}$
, over which the current deepens and slows. It is possible for these flows to become choked also if the height at their front becomes too large. The transition to ‘choked shocks’ occurs by simultaneously solving (A 5), together with shock conditions (3.22) and (3.26). This leads to the curve
$\mathscr{C}_{S}$
in the
$(\hat{Q},F_{0})$
plane (figure 15).
The calculations presented in this Appendix demonstrate that our model may be used with any dynamic front condition, including those that provide smaller Froude numbers than predicted by Benjamin’s formula. We note, however, that since the governing equations and balances across internal jumps are based upon inviscid dynamics and do not include effects such as velocity shear and mixing, that it may be inconsistent to apply empirical or semi-empirical formulae, which implicitly incorporate them at the front of the motion. We emphasise that the results presented in the main body of the paper and in § A.1 are self-contained and are closed rigorously without any adjustable constants and that the front condition proposed by Benjamin (2.19) or by circulation conservation (A 1) are self-contained, consistent outgrowths of the same balances that are used for the equations governing the entire current.
Appendix B. Two-layer, Boussinesq, dambreak flow
In this appendix we consider the related problem of two-layer dambreak flow in a channel, a problem that was analysed by Rottman & Simpson (Reference Rottman and Simpson1983) and Ungarish (Reference Ungarish2009). In this study, using the methods presented in § 3, the solutions emerge analytically through the solution of algebraic equations. These dambreak solutions thus form the two-layer counterparts of the dambreak solution derived by Ritter (Reference Ritter1892), which has proven an invaluable prototype for investigating spatially and temporally varying flows.
In terms of the dimensionless governing equations developed in this study, this situation corresponds to
$\hat{Q}=0$
and the initial conditions within the now infinite channel are given by
$u=0$
and
$h=h_{0}$
for
$x<0$
and
$h=0$
for
$x>0$
(see figure 16). For
$t>0$
, the flow evolves so that the dense fluid slumps downstream and a disturbance propagates upstream into the lock. When the lock depth is much less than the ambient
$(h_{0}\ll 1)$
, the motion is governed to leading order by the single-layer shallow-water equations and the dambreak solution is simply derived (Rottman & Simpson Reference Rottman and Simpson1983; Hogg Reference Hogg2006; Ungarish Reference Ungarish2009). When
$h_{0}$
is not negligibly small and the fluid is Boussinesq
$(S\ll 1)$
, complications arise, although the complete solution may still be constructed using the method of characteristics. This is what we undertake here to establish analytical expressions for the velocity and height fields.
Figure 16. The configuration of dambreak flow. The fluid is initially at rest and of depth
$h_{0}$
behind a lock gate at
$x=0$
in an infinitely long channel of dimensionless depth 1. The flow is initiated by the instantaneous removal of the lock gate at
$t=0$
.
First, the motion downstream
$(x>0)$
is governed by a front condition and could potentially be choked as described in § 2. Thus at
$x=x_{N}(t)$
,
$h_{N}\leqslant h_{c}$
and
$u_{N}=\mathscr{F}(h_{N})$
. The lead characteristic propagating upstream
$(x<0)$
into the region
$u=0$
and
$h=h_{0}$
satisfies
$$\begin{eqnarray}\frac{\text{d}x}{\text{d}t}=-\sqrt{h_{0}(1-h_{0})}.\end{eqnarray}$$
This expression is monotonically decreasing with
$h_{0}$
up to
$h_{0}=1/2$
, but if
$h_{0}>1/2$
there is no longer a smooth way in which to connect the solution with the lock conditions. This problem was identified by Rottman & Simpson (Reference Rottman and Simpson1983) and resolved by Klemp, Rotunno & Skamarock (Reference Klemp, Rotunno and Skamarock1994) and Ungarish (Reference Ungarish2009), who showed how to insert a jump that moves upstream as the leading disturbance with velocity,
$V_{s}$
$(V_{s}<0)$
, across which mass and momentum flux are conserved in each layer and energy is conserved in the upper (expanding) layer. Ungarish (Reference Ungarish2009) showed that
$$\begin{eqnarray}V_{s}^{2}=\frac{h_{0}^{2}-h_{1}^{2}}{\displaystyle \frac{2h_{0}^{2}-h_{1}}{h_{1}}+(1-2h_{1})\left(\displaystyle \frac{1-h_{0}}{1-h_{1}}\right)^{2}},\end{eqnarray}$$
where
$h_{1}$
is the depth of the dense fluid immediately downstream of the jump. The solution is completed by identifying the bore with the leading upstream characteristic and this provides a second condition for
$h_{1}$
and
$V_{s}$
. Simultaneously (B 2) and
$V_{s}={\it\lambda}_{-}(u_{1},h_{1})$
where
$u_{1}$
is the velocity of the lower layer just downstream of the bore yields the following equation for
$h_{1}$
in terms of the lock depth,
$h_{0}$
,
$$\begin{eqnarray}h_{1}^{4}-4h_{0}h_{1}^{3}+(-3h_{0}^{2}+6h_{0})h_{1}^{2}+(3h_{0}^{2}-2h_{0})h_{1}-h_{0}^{2}=0.\end{eqnarray}$$
Note that when
$h_{0}=1/2$
,
$h_{1}=1/2$
; thus the jump vanishes precisely when the continuous solution re-emerges. Furthermore, when
$h_{0}=1$
,
$h_{1}=1-h_{c}$
; for a full-depth lock, the upstream flow is choked and the fluid evolves asymmetrically about the origin, such that
$h(-x)=1-h(x)$
and
$u(-x)=u(x)h(x)/(1-h(x))$
.
The complete solution then comprises a rarefaction of
${\it\lambda}_{-}$
characteristics centred at the origin, with a bore propagating upstream if
$h_{0}>1/2$
and potentially choked flow at the front. First, we construct the continuous solutions
$(h_{0}<1/2)$
: within the rarefaction of
${\it\lambda}_{-}$
characteristics, the Riemann invariant
$R_{+}$
is constant determined by the lock conditions. Thus
$$\begin{eqnarray}-\!\sin ^{-1}(1-2h_{0})=\sin ^{-1}\left(\frac{u}{1-h}\right)-\sin ^{-1}(1-2h),\end{eqnarray}$$
and this applies at the front so that
$R_{+}(0,h_{0})=R_{+}(\mathscr{F}(h_{N}),h_{N})$
. The characteristics are straight lines in the
$(x,t)$
plane, given by
$$\begin{eqnarray}y=\frac{x}{t}=u-{\it\alpha}-({\it\alpha}^{2}+(1-{\it\beta})h)^{1/2}\quad \text{for }y_{a}<y<y_{b},\end{eqnarray}$$
where
$y_{a}={\it\lambda}_{-}(0,h_{0})=-\sqrt{h_{0}(1-h_{0})}$
and
$y_{b}={\it\lambda}_{-}(\mathscr{F}(h_{N}),h_{N})$
and
${\it\alpha}$
and
${\it\beta}$
are given by (2.12) with
$\hat{Q}=0$
. The solution is completed by a uniform region with
$h=h_{N}$
and
$u=\mathscr{F}(h_{N})$
for
$y_{b}<y<y_{N}\equiv u_{N}$
.
If
$h_{0}>1/2$
then the solution includes a discontinuity. Given
$h_{0}$
, (B 2) can be solved to determine
$h_{1}$
. Then replacing (B 4), we have
$$\begin{eqnarray}R_{+}(u_{1},h_{1})=R_{+}(u,h).\end{eqnarray}$$
The flow first becomes choked when
$R_{+}(u_{1},h_{1})=R_{+}(u_{c},h_{c})$
and this determines a critical lock depth
$h_{0}=h_{0m}=0.7965$
. Thus if
$h_{0m}>h_{0}>1/2$
, the flow features a bore moving upstream with velocity
$V_{s}$
; a rarefaction of
${\it\lambda}^{-}$
characteristics centred at the origin for
$V_{s}<x/t<y_{b}$
, within which the flow fields satisfy (B 5) and
$$\begin{eqnarray}R_{+}(u_{1},h_{1})=R_{+}(\mathscr{F}(h_{N}),h_{N})\end{eqnarray}$$
and
$y_{b}={\it\lambda}_{+}(\mathscr{F}(h_{N}),h_{N})$
; and a uniform region at the front of the current within which
$u=u_{N}$
and
$h=h_{N}$
for
$y_{b}<x/t<y_{N}\equiv u_{N}$
.
Figure 17. The height,
$h$
and velocity,
$u$
as functions of
$x/t$
of the lower layer in a dambreak flow for lock depth
$h_{0}=0.2,0.4,0.6,0.8,1.0$
.
Finally, if
$1>h_{0}>h_{0m}$
then the flow is choked at the front and there is a rarefaction of
${\it\lambda}_{+}$
characteristics within which the flow accelerates to match the front condition. Between the rarefactions at the front and rear of the currents, there is a uniform region within which the flow field is spatially uniform and given by
$u=u_{m}$
and
$h=h_{m}$
. Thus we find
$$\begin{eqnarray}R_{+}(u_{1},h_{1})=R_{+}(u_{m},h_{m})\quad \text{and}\quad R_{-}(u_{c},h_{c})=R_{-}(u_{m},h_{m}).\end{eqnarray}$$
It is then straightforward to simultaneously solve (B 8) to determine
$u_{m}$
and
$h_{m}$
and then the solution is given by
$$\begin{eqnarray}\displaystyle x/t={\it\lambda}_{-}(u,h)\quad \text{for }V_{s}<y<y_{b}, & & \displaystyle\end{eqnarray}$$
and
$$\begin{eqnarray}\displaystyle x/t={\it\lambda}_{+}(u,h)\quad \text{for }y_{c}<y<y_{N}, & & \displaystyle\end{eqnarray}$$
where
$y_{b}={\it\lambda}_{-}(u_{m},h_{m})$
,
$y_{c}={\it\lambda}_{-}(u_{m},h_{m})$
and
$y_{N}=u_{c}$
.
The case of a full-depth lock release
$(h_{0}=1)$
is particularly simple. By symmetry
$h_{m}=1/2$
and then from (B 8)
$$\begin{eqnarray}\sin ^{-1}(2u_{m})=\sin ^{-1}\left(\frac{u_{c}}{1-h_{c}}\right)+\sin ^{-1}(1-2h_{c}).\end{eqnarray}$$
This allows the evaluation of
$u_{m}=h_{c}^{1/2}(3-4h_{c})/2=0.4746$
. Also, the edge of the expansion fan may be evaluated as
$$\begin{eqnarray}y_{c}=-y_{b}={\it\lambda}_{+}(u_{m},1/2)=\left({\textstyle \frac{1}{4}}-u_{m}^{2}\right)^{1/2}={\textstyle \frac{1}{2}}(4h_{c}-1)(1-h_{c})^{1/2}=0.1572.\end{eqnarray}$$
The solutions for the height and velocity fields are plotted in figure 17 for a range of lock depths.
It is noteworthy that the jump at the front of the current and the internal bore (when present) are dissipative and so all dambreak currents are always dissipative. It is also worth reiterating that these solutions apply to reservoirs of infinite extent; if the reservoir were finite, the motion would be modified when the leading upstream characteristic (or bore) reached a back wall and is reflected, ultimately catching up with and modifying the motion of the front. These is no counterpart of this phenomenon for gravity currents driven by a sustained source, as analysed in the main body of this paper.
