1. Introduction
Gravity currents represent a well-known class of predominantly horizontal flows that form as a result of horizontal density differences and the associated hydrostatic pressure gradients (Benjamin Reference Benjamin1968; Simpson Reference Simpson1997). They encompass a variety of atmospheric and oceanic flows such as turbidity currents, powder snow avalanches and thunderstorm outflows (Moncrieff Reference Moncrieff1978; Hopfinger Reference Hopfinger1983; Meiburg & Kneller Reference Meiburg and Kneller2010; Ungarish Reference Ungarish2010; Linden Reference Linden, Chassignet, Cenedese and Verron2012). Frequently, such gravity currents interact with ambient shear. For example, close to the ground severe thunderstorms often exhibit layers of relatively cold air due to the evaporation and melting of hydro-meteors (Bryan & Rotunno Reference Bryan and Rotunno2014a ). The resulting cold pool of air can form a gravity current which then interacts with the ambient shear to produce complex flow structures. Several detailed observations suggest that the coupling between the cold pool of a thunderstorm outflow and environmental shear may produce squall lines (Moncrieff Reference Moncrieff1978, Reference Moncrieff1992; Rotunno, Klemp & Weisman Reference Rotunno, Klemp and Weisman1988; Xu & Moncrieff Reference Xu and Moncrieff1994). This indicates that gravity currents interacting with background shear play an important role in the organization of vertical convection (Weckwerth & Wakimoto Reference Weckwerth and Wakimoto1992). There have been several efforts to model and predict the dynamics of gravity currents interacting with background shear, e.g. Xu (Reference Xu1992), Xu & Moncrieff (Reference Xu and Moncrieff1994), Xu, Xue & Droegemeer (Reference Xu, Xue and Droegemeer1996), Xue, Xu & Droegemeier (Reference Xue, Xu and Droegemeier1997), Xue (Reference Xue2000, Reference Xue2002) and Bryan & Rotunno (Reference Bryan and Rotunno2014a ,Reference Bryan and Rotunno b ). These investigations focus on predicting the front velocity of gravity currents intruding into sheared environments as a function of the flow configuration and gravity current height. The ambient shear is frequently modelled either as a two-layer stream with a velocity jump or as a single layer with constant shear, and empirical assumptions are commonly required to close the resulting system of equations. Often, these closure assumptions are based on Bernoulli’s equation.
The current investigation extends the vorticity-based modelling approach that we recently introduced for unsheared gravity currents and internal bores (Borden & Meiburg Reference Borden and Meiburg2013a ,Reference Borden and Meiburg b ) to gravity currents with shear. The key new feature of these vorticity-based models is that, in contrast to earlier modelling approaches (von Karman Reference von Karman1940; Benjamin Reference Benjamin1968; Shin, Dalziel & Linden Reference Shin, Dalziel and Linden2004; Nokes et al. Reference Nokes, Davidson, Stepien, Veale and Oliver2008), they do not require any empirical energy-based closure assumptions. This is accomplished by enforcing the conservation of vertical momentum, along with the usual conservation of mass and horizontal momentum. Towards this end, the integral form of the vorticity conservation equation is employed, i.e. a linear combination of the conservation equations for horizontal and vertical momentum. For traditional gravity currents and internal bores without shear, this approach was seen to yield better agreement with direct numerical simulation (DNS) results than the earlier models with energy-based closures.
Section 2 defines the physical problem set-up and derives the governing equations for both the two-layer and the single-layer configurations. Two alternative approaches are outlined for assessing the feasibility of the theoretical solutions, either via determining the loss of mechanical energy along specific streamlines or by analysing the depth-integrated loss of energy from the inflow to the outflow boundary. Furthermore, the relationship between these two approaches is discussed. In § 3, we analyse the gravity current dynamics as a function of the fractional current height and the imposed shear. Interestingly, the theory allows for the existence of gravity currents thicker than half the channel height. Section 4 compares the theoretical model predictions with DNS simulation results and confirms the existence of such thick gravity currents. For high Reynolds numbers, we furthermore investigate the effects of Kelvin–Helmholtz instabilities and mixing on the current height and its propagation velocity. Section 5 summarizes the findings and presents the main conclusions.
2. Physical problem
We will analyse the dynamics of gravity currents propagating into sheared environments for two limiting conditions. Initially, we will consider a two-layer (TL) configuration, in which the ambient fluid consists of two constant-velocity streams separated by a shear layer. Subsequently, we will focus on the related situation of a gravity current propagating into an ambient fluid of constant vorticity, the so-called single-layer (SL) configuration.
2.1. Two-layer configuration
Consider a gravity current of height
$h$
and density
${\it\rho}_{1}$
propagating in a horizontal channel of height
$H$
filled with ambient fluid of density
${\it\rho}_{0}$
, as shown in figure 1. In the laboratory reference frame, the gravity current travels from right to left with constant velocity
$U_{g}$
, into an environment consisting of two uniform streams whose velocities are
${\rm\Delta}U/2$
and
$-{\rm\Delta}U/2$
respectively. This defines the laboratory reference frame. Far ahead of the current front, the inlet conditions
$h_{1i}$
,
$h_{2i}$
and
${\rm\Delta}U$
are given.
Figure 1. The two-layer configuration: a gravity current (grey region) propagates into an ambient consisting of two uniform streams separated by a vortex sheet of strength
${\rm\Delta}U$
.
To facilitate the analysis, we shift to a reference frame moving with the gravity current front, so that the flow field is steady and the fluid inside the current is at rest. The goal of our analysis is to determine the gravity current velocity
$U_{g}$
, as well as the conditions far downstream of the current front in terms of the layer heights (
$h_{1o}$
,
$h_{2o}$
) and velocities (
$U_{1o}$
,
$U_{2o}$
). For the steady-state configuration in the reference frame moving with the gravity current, the continuity equations for the two ambient fluid layers yield
$$\begin{eqnarray}\displaystyle & \displaystyle \left(U_{g}+\frac{{\rm\Delta}U}{2}\right)h_{1i}=U_{1o}h_{1o}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \left(U_{g}-\frac{{\rm\Delta}U}{2}\right)h_{2i}=U_{2o}h_{2o}. & \displaystyle\end{eqnarray}$$
In order to obtain additional equations, we now invoke the conservation of momentum. Rather than applying the
$x$
- and
$y$
-momentum equations separately, at this stage we focus on the vorticity conservation equation, i.e. a linear combination of the
$x$
- and
$y$
-momentum equations. This approach has the advantage of eliminating the pressure variable from the analysis. It should be noted that, if information on the pressure is desired, we can solve the
$x$
-momentum equation for the pressure at a later stage, after all other flow variables have been determined. In this way, by using the vorticity equation we essentially decouple the pressure problem, cf. also Borden & Meiburg (Reference Borden and Meiburg2013a
,Reference Borden and Meiburg
b
).
For two-dimensional steady-state Boussinesq flow, the vorticity field is governed by
$$\begin{eqnarray}\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\bf\omega}=-g^{\prime }\frac{\partial {\it\rho}^{\ast }}{\partial x}+{\it\nu}{\rm\nabla}^{2}{\it\omega}.\end{eqnarray}$$
Here,
${\it\omega}$
,
$g^{\prime }$
,
${\it\rho}^{\ast }$
and
${\it\nu}$
represent the out-of-plane vorticity component
$\partial v/\partial x-\partial u/\partial y$
, the reduced gravity
$g({\it\rho}_{1}-{\it\rho}_{0})/{\it\rho}_{0}$
, the non-dimensional density
$({\it\rho}-{\it\rho}_{0})/({\it\rho}_{1}-{\it\rho}_{0})$
and the kinematic viscosity respectively. It should be noted that we do not make any assumptions regarding the flow being hydrostatic.
For the purpose of determining the gravity current velocity, we assume that the viscosity is sufficiently small for the vorticity to remain confined to narrow layers along the interfaces, so that there is no diffusive interaction between these vorticity layers and no diffusion of vorticity across the top and bottom boundaries. Equation (2.3) can then be integrated for the control volume
$\text{AO}^{\prime }\text{C}^{\prime }\text{C}$
in figure 1. Using the divergence theorem we obtain
$$\begin{eqnarray}\oint {\it\omega}\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{n}\,\text{d}{\it\Gamma}=-\iint _{A}g^{\prime }\frac{\partial {\it\rho}^{\ast }}{\partial x}\,\text{d}A.\end{eqnarray}$$
Here,
$A$
,
${\it\Gamma}$
and
$\boldsymbol{n}$
denote the control volume area, the control volume boundary and the unit outer normal vector respectively. We assume that the normal derivative of the shear stress along the top and bottom boundaries vanishes. At the in- and outflow boundaries, the flow is approximately unidirectional and normal to the boundaries, so that
$\boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\omega}=0$
. Consequently, the diffusion of vorticity across the in- and outflow boundaries can be neglected as well. Equation (2.4) therefore states that the rate at which vorticity is convected out of the control volume at the downstream boundary equals the rate at which it is convected into the control volume at the upstream boundary, plus the rate at which it is generated inside the control volume as a result of baroclinic vorticity production.
We can evaluate the vorticity flux associated with a vortex sheet as the vortex sheet strength multiplied by its principal velocity, cf. Saffman (Reference Saffman1992) and Borden & Meiburg (Reference Borden and Meiburg2013a
). For the vortex sheet separating the two layers of constant-density ambient fluid (line
$\text{BB}^{\prime }$
in figure 1), we obtain
$$\begin{eqnarray}U_{g}{\rm\Delta}U=\frac{U_{1o}+U_{2o}}{2}(U_{1o}-U_{2o}).\end{eqnarray}$$
For the vortex sheet along the interface between the gravity current and ambient fluid, the vorticity flux is balanced by the baroclinic vorticity generation. Since the fluid inside the current is at rest, this yields
$$\begin{eqnarray}{\textstyle \frac{1}{2}}U_{1o}^{2}=g^{\prime }h.\end{eqnarray}$$
One final equation is obtained from geometrical considerations:
$$\begin{eqnarray}h_{1o}+h_{2o}+h=H.\end{eqnarray}$$
The two continuity equations for the ambient fluid layers (2.1) and (2.2), along with the two vorticity equations for the interfaces (2.5) and (2.6) as well as the geometrical constraint (2.7), provide five equations for the five unknowns
$U_{g}$
,
$U_{1o}$
,
$U_{2o}$
,
$h_{1o}$
and
$h_{2o}$
. For any given parameter set (
${\rm\Delta}U$
,
$h$
,
$h_{1i}$
,
$h_{2i}$
,
$H$
and
$g^{\prime }$
), this set of nonlinear equations can be solved iteratively. It should be noted that we did not have to invoke any assumptions regarding energy conservation along any of the streamlines in the flow, in contrast to previous studies, e.g. Benjamin (Reference Benjamin1968), Shin et al. (Reference Shin, Dalziel and Linden2004) and Nokes et al. (Reference Nokes, Davidson, Stepien, Veale and Oliver2008). In the present investigation, we focus primarily on cases in which both ambient currents flow from left to right in the reference frame moving with the gravity current, cf. figure 1, although in § 4.5 we will briefly comment on situations in which one of the ambient streams flows from right to left. Thus, the following conditions will always hold:
$$\begin{eqnarray}\displaystyle & \displaystyle U_{g}\pm \frac{{\rm\Delta}U}{2}\geqslant 0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & U_{1o},U_{2o}\geqslant 0. & \displaystyle\end{eqnarray}$$
2.2. Single-layer configuration
Corresponding conservation arguments can now be invoked for the single-layer case SL, cf. figure 2, in which the ambient flow has constant shear across the channel height. For a given gravity current height
$h$
, reduced gravity
$g^{\prime }$
and imposed shear at the inlet
${\rm\Delta}U=U_{i}(y=0)-U_{i}(y=H)$
, we need to determine the gravity current velocity
$U_{g}$
along with the velocity profile at the outlet
$U_{o}(y)$
.
Figure 2. The single-layer configuration: a gravity current (grey region) propagates into ambient fluid with constant shear.
The continuity equation for the ambient takes the form
$$\begin{eqnarray}\int _{h}^{H}U_{o}(y)\,\text{d}y=\int _{0}^{H}U_{i}(y)\,\text{d}y=U_{g}H,\end{eqnarray}$$
where
$U_{i}(y)$
has the linear form
$$\begin{eqnarray}U_{i}(y)=U_{g}+\frac{{\rm\Delta}U}{2}-\frac{{\rm\Delta}U}{H}y.\end{eqnarray}$$
Since the ambient fluid has constant density, no baroclinic vorticity is produced within the ambient stream, and the slope of the outflow velocity profile equals that of the inflow profile, i.e.
$\text{d}U_{i}/\text{d}y=\text{d}U_{o}/\text{d}y$
. Thus, the outflow velocity profile has the form
$$\begin{eqnarray}U_{o}(y)=U_{o}|_{y=h}-\frac{{\rm\Delta}U}{H}(y-h)\quad y\geqslant h.\end{eqnarray}$$
The outflow velocity at
$U_{o}|_{y=h}$
can be evaluated from the vorticity conservation equation along the interface (
$\text{OA}^{\prime }$
in figure 2),
$$\begin{eqnarray}\left.{\textstyle \frac{1}{2}}U_{o}^{2}\right|_{y=h}=g^{\prime }h.\end{eqnarray}$$
Unlike for the TL case, we can obtain explicit relations to compute the values of the two unknowns
$U_{g}$
and
$U_{o}(y=h)$
for the SL configuration,
$$\begin{eqnarray}\displaystyle & U_{o}|_{y=h}=\sqrt{2g^{\prime }h}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle U_{g}=\left(1-\frac{h}{H}\right)\sqrt{2g^{\prime }h}-\frac{1}{2}\left(1-\frac{h}{H}\right)^{2}{\rm\Delta}U. & \displaystyle\end{eqnarray}$$
$U_{g}\geqslant 0$
. This imposes an upper bound on the shear magnitude in the form 
$$\begin{eqnarray}{\rm\Delta}U\leqslant \frac{2\sqrt{2g^{\prime }h}}{1-h/H}.\end{eqnarray}$$
2.3. Assessment of physically feasible solutions
In order to assess whether the solutions obtained in §§ 2.1 and 2.2 are physically feasible, we need to check whether or not they require an external energy input. Towards this end, we now perform a head loss analysis along selected streamlines.
Applying Bernoulli’s equation with a head loss between the two points
$\text{C}$
and
$\text{C}^{\prime }$
along the top wall in figures 1 and 2 yields
$$\begin{eqnarray}p_{C}+{\textstyle \frac{1}{2}}{\it\rho}_{0}U_{C}^{2}=p_{C^{\prime }}+{\textstyle \frac{1}{2}}{\it\rho}_{0}U_{C^{\prime }}^{2}+{\it\delta}_{CC^{\prime }}.\end{eqnarray}$$
For the system to be physically feasible, the head loss term
${\it\delta}_{CC^{\prime }}$
should be non-negative. From the analysis in §§ 2.1 and 2.2 we know
$U_{C^{\prime }}$
. To obtain
${\it\delta}_{CC^{\prime }}$
, however, we also need to evaluate the pressure difference between these two points. This can be accomplished via the conservation equation for horizontal momentum,
$$\begin{eqnarray}\int _{0}^{H}(p_{i}+{\it\rho}_{0}U_{i}^{2})\,\text{d}y=\int _{0}^{h}p_{o}\,\text{d}y+\int _{h}^{H}(p_{o}+{\it\rho}_{0}U_{o}^{2})\,\text{d}y,\end{eqnarray}$$
in conjunction with the hydrostatic pressure relations far up- and downstream of the gravity current front,
$$\begin{eqnarray}\displaystyle & p_{i}(y)=p_{C}+{\it\rho}_{0}g(H-y), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & p_{o}(y)=\left\{\begin{array}{@{}ll@{}}p_{C^{\prime }}+{\it\rho}_{0}g(H-y)\quad & \text{if }y>h,\\ p_{C^{\prime }}+{\it\rho}_{0}g(H-h)+{\it\rho}_{1}g(h-y)\quad & \text{if }y\leqslant h.\end{array}\right. & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle {\rm\Delta}p|^{TL} & = & \displaystyle p_{C^{\prime }}-p_{C}=-\frac{1}{H}\left[\vphantom{\mathop{\sum }_{1}^{1}}({\it\rho}_{1}-{\it\rho}_{0})g\frac{h^{2}}{2}+{\it\rho}_{0}h_{1o}U_{1o}^{2}+{\it\rho}_{0}h_{2o}U_{2o}^{2}\right.\nonumber\\ \displaystyle & & \displaystyle \left.-\,{\it\rho}_{0}h_{1i}\left(U_{g}+\frac{{\rm\Delta}U}{2}\right)^{2}-{\it\rho}_{0}h_{2i}\left(U_{g}-\frac{{\rm\Delta}U}{2}\right)^{2}\vphantom{\mathop{\sum }_{1}^{1}}\right],\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle {\rm\Delta}p|^{SL} & = & \displaystyle p_{C^{\prime }}-p_{C}=-\frac{1}{H}\left[\vphantom{\mathop{\sum }_{1}^{1}}({\it\rho}_{1}-{\it\rho}_{0})g\frac{h^{2}}{2}-{\it\rho}_{0}\frac{H}{12}({\rm\Delta}U^{2}+12U_{g}^{2})\right.\nonumber\\ \displaystyle & & \displaystyle \left.+\,{\it\rho}_{0}(H-h)\left(\frac{1}{2}g^{\prime }h+\frac{1}{3}\left\{{\rm\Delta}U\left(1-\frac{h}{H}\right)-\frac{3}{2}\sqrt{2g^{\prime }h}\right\}^{2}\right)\vphantom{\mathop{\sum }_{1}^{1}}\right].\end{eqnarray}$$
${\it\delta}_{CC^{\prime }}$
along the top wall via (2.17).
It is straightforward to show that the head loss along any streamline in the ambient is identical to the head loss along the top wall, e.g.
${\it\delta}_{CC^{\prime }}={\it\delta}_{AA^{\prime }}$
(for both SL and TL cases), and
${\it\delta}_{CC^{\prime }}={\it\delta}_{BB^{\prime }}$
in the TL case. In the following, this is demonstrated for the SL case. Using Bernoulli’s equation along
$\text{A}{-}\text{A}^{\prime }$
, the head loss is obtained as
$$\begin{eqnarray}{\it\delta}_{AA^{\prime }}=p_{A}-p_{A^{\prime }}+{\textstyle \frac{1}{2}}{\it\rho}_{0}(U_{A}^{2}-U_{A^{\prime }}^{2})-{\it\rho}_{0}gh.\end{eqnarray}$$
Using the hydrostatic pressure relations (2.19) and (2.20) far up- and downstream of the gravity current front, we obtain the pressure drop between points
$\text{A}$
and
$\text{A}^{\prime }$
as
$$\begin{eqnarray}p_{A}-p_{A^{\prime }}=p_{C}-p_{C^{\prime }}+{\it\rho}_{0}gh.\end{eqnarray}$$
We recall that for the SL configuration, the velocities in (2.23) are given as
$$\begin{eqnarray}U_{A}=U_{g}+\frac{{\rm\Delta}U}{2},\quad U_{A^{\prime }}=\sqrt{2g^{\prime }h},\end{eqnarray}$$
$$\begin{eqnarray}U_{C}=U_{g}-\frac{{\rm\Delta}U}{2},\quad U_{C^{\prime }}=\sqrt{2g^{\prime }h}-\left(1-\frac{h}{H}\right){\rm\Delta}U,\end{eqnarray}$$
where the gravity current velocity
$U_{g}$
for the SL case can be computed via (2.15). By substituting (2.24)–(2.26) into (2.23) and (2.17) and simplifying, we obtain
$$\begin{eqnarray}{\it\delta}_{AA^{\prime }}={\it\delta}_{CC^{\prime }}.\end{eqnarray}$$
Following similar arguments, the same result can be obtained for the TL configuration.
An alternate strategy for assessing the physical feasibility of the flow is to compare the depth-integrated energy loss from the inflow to the outflow. This approach can be useful when there exists a reverse flow over a small portion of the in- or outflow boundaries. We define the net energy loss in the control volume as
$$\begin{eqnarray}\displaystyle {\rm\Delta}{\dot{E}} & = & \displaystyle {\dot{E}}_{i}-{\dot{E}}_{o}=\int _{0}^{H}\left(p+\frac{1}{2}{\it\rho}_{0}U_{i}^{2}(y)+{\it\rho}_{0}gy\right)U_{i}(y)\,\text{d}y\nonumber\\ \displaystyle & & \displaystyle -\,\int _{h}^{H}\left(p+\frac{1}{2}{\it\rho}_{0}U_{o}^{2}(y)+{\it\rho}_{0}gy\right)U_{o}(y)\,\text{d}y.\end{eqnarray}$$
Using the hydrostatic pressure relations at the inlet and outlet, we arrive at
$$\begin{eqnarray}{\rm\Delta}{\dot{E}}=-(p_{C^{\prime }}-p_{C})U_{g}H+\frac{1}{2}{\it\rho}_{0}\int _{0}^{H}U_{i}^{3}(y)\,\text{d}y-\frac{1}{2}{\it\rho}_{0}\int _{h}^{H}U_{o}^{3}(y)\,\text{d}y.\end{eqnarray}$$
It is then straightforward to show for the SL configuration (and similarly for the TL configuration) that the depth-integrated energy loss is related to the head loss via
$$\begin{eqnarray}{\rm\Delta}{\dot{E}}=-(p_{C^{\prime }}-p_{C})U_{g}H+{\textstyle \frac{1}{2}}{\it\rho}_{0}(U_{C}^{2}-U_{C^{\prime }}^{2})U_{g}H=U_{g}H{\it\delta}_{CC^{\prime }}.\end{eqnarray}$$
Equation (2.30) suggests that the head loss and depth-integrated energy loss arguments yield the same result, as far as the feasibility of the solution is concerned.
3. Results
We will now employ the above theoretical framework in order to analyse single- and two-layer configurations. Subsequently, we will compare these theoretical results against DNS simulations. We will pay particular attention to the influence of the inflow conditions on the velocity and head loss of the current, and we will demonstrate the existence of current heights larger than half the channel height. Before doing so, we render the governing equations dimensionless by using as characteristic scales the channel height
$H$
, the buoyancy velocity
$u_{b}=\sqrt{g^{\prime }H}$
and the density difference
${\it\rho}_{1}-{\it\rho}_{0}$
:
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{x}^{\ast }=\frac{\boldsymbol{x}}{H}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{u}^{\ast }=\frac{\boldsymbol{u}}{\sqrt{g^{\prime }H}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle t^{\ast }=\frac{t}{H/\sqrt{g^{\prime }H}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle p^{\ast }=\frac{p}{{\it\rho}_{0}g^{\prime }H}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\dot{E}}^{\ast }=\frac{{\dot{E}}}{{\it\rho}_{0}g^{\prime 3/2}H^{5/2}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\it\rho}^{\ast }=\frac{{\it\rho}-{\it\rho}_{0}}{{\it\rho}_{1}-{\it\rho}_{0}}. & \displaystyle\end{eqnarray}$$
$\text{}^{\ast }$
symbol is employed to indicate a dimensionless quantity.
3.1. Two-layer configuration
We begin by considering case TL with
$h_{1i}^{\ast }=h_{2i}^{\ast }=0.5$
, cf. figure 1. Figure 3 shows the gravity current velocity
$U_{g}^{\ast }$
and the head loss along the upper wall
${\it\delta}_{CC^{\prime }}^{\ast }$
as functions of the gravity current height
$h^{\ast }$
and the inflow shear magnitude
${\rm\Delta}U^{\ast }$
. Physically feasible flow conditions, i.e. those resulting in a positive head loss, are shaded grey. For vanishing shear, we recover Benjamin’s classical result which states that a current occupying half the channel height has a front velocity
$U_{g}^{\ast }=0.5$
and conserves energy, cf. case
$\text{T}_{0}$
in figure 3. For a detailed derivation of the energy analysis we refer to § 3.2.1, which provides such an analysis for the single-layer configuration. In the absence of shear, as discussed by Benjamin (Reference Benjamin1968) and Borden & Meiburg (Reference Borden and Meiburg2013a
), only gravity currents of less than half the channel height are physically possible. However, for positive values of
${\rm\Delta}U^{\ast }$
, when the lower ambient stream moves faster than the upper one, gravity currents with significantly larger heights become feasible. Benjamin (Reference Benjamin1968) points out that such flows require external energy input into the ambient. In light of (2.30) it becomes clear that the presence of shear modifies the pressure jump
$p_{C^{\prime }}-p_{C}$
over the current front, which in turn affects the kinetic energy of the ambient current. If the upper stream is the faster one, on the other hand, currents with thickness greater than half the channel height cannot form. The region for which results are shown is bounded by lines determined from the condition
$U_{g}\pm {\rm\Delta}U/2\geqslant 0$
, so that in the moving reference frame both ambient streams flow from left to right.
Figure 3. Configuration TL with inlet layer heights of
$h_{1i}^{\ast }=h_{2i}^{\ast }=0.5$
: phase-space diagram for (a) the gravity current velocity
$U_{g}^{\ast }$
and (b) the head loss along the upper wall
${\it\delta}_{CC^{\prime }}^{\ast }$
, as functions of
${\rm\Delta}U^{\ast }$
and
$h^{\ast }$
. The grey region indicates physically possible solutions, i.e. solutions with a positive head loss. For positive shear, when the lower ambient stream moves faster than the upper one, currents with thickness larger than half the channel height can form.
To investigate the influence of the inlet properties, we also show the gravity current velocities for
$h_{1i}^{\ast }=0.25$
and
$h_{1i}^{\ast }=0.75$
, cf. figure 4. As the thickness of the lower ambient current decreases, even thicker gravity currents become physically feasible. Conversely, the thinner the lower ambient stream, the less shear is required to obtain a given gravity current thickness
$h^{\ast }>0.5$
.
Figure 4. Phase-space diagram for the gravity current velocity
$U_{g}^{\ast }$
as function of
${\rm\Delta}U^{\ast }$
and
$h^{\ast }$
for (a)
$h_{1i}^{\ast }=0.25$
and (b)
$h_{1i}^{\ast }=0.75$
. Thinner lower ambient streams allow for a wider range of physically feasible gravity current thicknesses.
The existence of a saddle point in all of the velocity diagrams of figures 3(a) and 4 is a result of two competing effects. On one hand, for constant shear the gravity current velocity exhibits a maximum at an intermediate height
$h^{\ast }$
. For vanishing shear specifically, the analyses by Benjamin (Reference Benjamin1968) and Borden & Meiburg (Reference Borden and Meiburg2013a
) place this maximum at
$h^{\ast }=0.347$
and
$h^{\ast }=0.333$
respectively (figure 5
a). Away from this value, the gravity current velocity starts to decrease. On the other hand, for fixed gravity current heights the current velocity displays a minimum at intermediate shear. Away from this minimum, an increase in the shear magnitude leads to faster gravity currents (figure 5
b). As suggested by figures 3(a) and 4, the location of the saddle point is strongly dependent on the flow parameters as well as the inlet heights.
Figure 5. Gravity current velocity for
$h_{1i}^{\ast }=h_{2i}^{\ast }=0.5$
computed for two different scenarios. (a) For a vanishing shear
${\rm\Delta}U^{\ast }=0$
, the current velocity has a maximum at intermediate current heights. (b) For a fixed value of
$h^{\ast }=0.5$
, the gravity current velocity reaches a minimum at an intermediate shear value. Together, these trends account for the existence of the saddle point in figure 3(a).
Figure 6 shows the minimum shear required for a gravity current to have a certain thickness, for different inlet height values
$h_{1i}^{\ast }$
. This minimum shear value is given by the velocity jump
${\rm\Delta}U^{\ast }$
for zero head loss. We observe that the required minimum shear is a strong function of
$h_{1i}^{\ast }$
. For a given current height, a thinner lower ambient stream translates into a significantly smaller value of
${\rm\Delta}U^{\ast }$
.
Figure 6. The minimum shear required to produce a gravity current of a certain height larger than half the channel height, for different lower ambient stream thicknesses
$h_{1i}^{\ast }$
.
Figure 7. Ratio of the outflow shear to the inflow shear of the ambient streams for the TL configuration with
$h_{1i}^{\ast }=h_{2i}^{\ast }=0.5$
.
Figure 7 displays the ratio of outflow to inflow shear,
${\rm\Delta}U_{o}^{\ast }/{\rm\Delta}U^{\ast }=(U_{1o}^{\ast }-U_{2o}^{\ast })/{\rm\Delta}U^{\ast }$
, as a function of the inflow shear
${\rm\Delta}U^{\ast }$
for various gravity current heights. We observe that the outflow shear
${\rm\Delta}U_{o}^{\ast }$
is always less than
${\rm\Delta}U^{\ast }$
. This is a consequence of the higher average velocity of the ambient streams at the outflow in comparison with the inflow. Vorticity conservation dictates that the product of this average velocity and the velocity difference between the two streams remain constant, so that a larger average velocity necessarily leads to a reduced velocity difference. The continuity equation requires that the overall acceleration of the two ambient streams is stronger for larger gravity current heights, so that the reduction in shear becomes more pronounced for larger values of
$h^{\ast }$
.
Interestingly, the ratio of outflow to inflow shear is larger for a given positive shear magnitude (lower ambient stream faster) than for the same negative shear magnitude (upper ambient stream faster). This can be understood in the following way. The outflow velocity of the lower ambient stream is a function of the gravity current height only, cf. (2.6). Consequently, if the lower ambient stream is slower at the inflow than the upper one, it will have to be accelerated more than if it is faster at the inflow. Continuity then tells us that if the upper stream is faster at the inflow, it will have to be accelerated less than if it is slower. As a result, the velocity difference between the ambient streams decreases more strongly if the lower ambient stream is slower at the inflow than the upper one, i.e. if
${\rm\Delta}U^{\ast }<0$
.
This scenario is confirmed by figure 8, which displays the influence of
${\rm\Delta}U^{\ast }$
on the outflow height
$h_{1o}^{\ast }$
of the lower ambient stream, and on the ratio
$h_{2o}^{\ast }/h_{1o}^{\ast }$
for different gravity current heights. If the lower ambient stream is the faster one, i.e. if
${\rm\Delta}U^{\ast }>0$
, it experiences weaker acceleration than if
${\rm\Delta}U^{\ast }<0$
, so that it retains a larger thickness. Geometry then dictates that
$h_{2o}^{\ast }$
has to decrease, so that the ratio
$h_{2o}^{\ast }/h_{1o}^{\ast }$
also decreases.
Figure 8. Influence of inflow shear on the outflow height
$h_{1o}^{\ast }$
of the lower layer (a) and on the ratio
$h_{2o}^{\ast }/h_{1o}^{\ast }$
(b) for several gravity current heights, for the TL configuration with
$h_{1i}^{\ast }=h_{2i}^{\ast }=0.5$
.
3.2. Single-layer configuration
Figure 9 displays the gravity current velocity and the head loss along the upper wall for the SL case with a linear inflow velocity profile, cf. figure 2. Similarly to the TL configuration, we observe that with a constant shear imposed at the inlet, gravity currents exist whose thickness is more than half the channel height. For a detailed derivation of the energy analysis, we refer to § 3.2.1 below. The range of parameters for which this occurs is quite similar to the corresponding TL configuration with
$h_{1i}^{\ast }=h_{2i}^{\ast }=0.5$
, cf. figure 3. Unlike the TL configuration (figure 3), the SL case does not exhibit a saddle point in the velocity phase-space diagram, figure 9(a), since the gravity current velocity in the SL configuration varies monotonically with the shear magnitude
${\rm\Delta}U$
in the ambient layer, cf. (2.15).
Figure 9. Configuration SL: phase-space diagram for (a) the gravity current velocity
$U_{g}^{\ast }$
and (b) the head loss
${\it\delta}_{CC^{\prime }}^{\ast }$
along the upper wall, as functions of shear
${\rm\Delta}U^{\ast }$
and current height
$h^{\ast }$
. The grey region indicates the physically feasible region based on positive values of the head loss
${\it\delta}_{CC^{\prime }}^{\ast }$
.
Figure 10 shows the depth-integrated energy loss, (2.30), for the SL configuration. As discussed earlier, (2.30) indicates that the head loss and depth-integrated energy loss arguments yield the same grey regions, as far as the feasibility of the solution is concerned. While the requirement of a
${\rm\Delta}{\dot{E}}\geqslant 0$
results in the same grey region as obtained via a head loss analysis (see figure 9), it further identifies the parameter regime of physically feasible solutions in which ambient stream flow demonstrates mixed inflow–outflow at the inlet region (dark region in figure 10). This will be further discussed in § 4.5.
Figure 10. Phase-space diagram for net energy loss (
${\rm\Delta}{\dot{E}}^{\ast }$
, see (2.29)) for the SL configuration. The light grey region indicates the physically feasible region according to both the energy difference and positive inlet/outlet velocities while the dark grey region indicates physically feasible solutions when one of the ambient streams flows from right to left.
3.2.1. Gravity currents thicker than half the channel height
For gravity currents propagating into quiescent ambient fluid, Benjamin (Reference Benjamin1968) discusses that without external energy input only current heights
$h/H\leqslant 1/2$
are possible, while currents with
$h/H>1/2$
require external energy input. On the other hand, our theoretical predictions in figures 3 and 9 indicate that, for positive shear values
${\rm\Delta}U>0$
, gravity currents thicker than half the channel height can form in both the TL and SL configurations. In order to investigate this issue further, figure 11 compares the pressure difference
$p_{C^{\prime }}-p_{C}$
, the head loss
${\it\delta}_{CC^{\prime }}$
and the energy loss
${\rm\Delta}{\dot{E}}$
for the classical case without shear with the SL case with shear. We remark that the analysis can be extended to the TL configuration in a similar fashion.
Figure 11. (a) Gravity current velocity
$U_{g}^{\ast }$
and pressure difference
${\rm\Delta}p^{\ast }$
along the top wall (see (3.7) and (3.9)). (b) Head loss
${\it\delta}_{CC^{\prime }}^{\ast }$
and energy loss
${\rm\Delta}{\dot{E}}^{\ast }$
. The solid, dotted and dash-dotted lines correspond to the cases of gravity currents without shear, with negative shear
${\rm\Delta}U^{\ast }=-0.4$
and with positive shear
${\rm\Delta}U^{\ast }=0.4$
respectively in the SL configuration. Only for positive shear can the head loss (and energy loss) be positive for
$h^{\ast }>1/2$
.
The pressure difference along the top wall for the zero-shear case reads
$$\begin{eqnarray}{\rm\Delta}p=p_{C^{\prime }}-p_{C}=-\frac{1}{H}\left[({\it\rho}_{1}-{\it\rho}_{0})g\frac{h^{2}}{2}+{\it\rho}_{0}(H-h)2g^{\prime }h-{\it\rho}_{0}HU_{g0}^{2}\right].\end{eqnarray}$$
Here,
$U_{g0}$
denotes the gravity current velocity given by (see Borden & Meiburg Reference Borden and Meiburg2013a
)
$$\begin{eqnarray}U_{g0}=\left(1-\frac{h}{H}\right)\sqrt{2g^{\prime }h}.\end{eqnarray}$$
For the sake of comparison, we rewrite the pressure drop obtained in the SL configuration (see (2.22)) as
$$\begin{eqnarray}\displaystyle {\rm\Delta}p|^{SL} & = & \displaystyle p_{C^{\prime }}-p_{C}=-\frac{1}{H}\left[\vphantom{\mathop{\sum }_{1}^{1}}({\it\rho}_{1}-{\it\rho}_{0})g\frac{h^{2}}{2}+{\it\rho}_{0}(H-h)2g^{\prime }h-{\it\rho}_{0}HU_{g0}^{2}\right.\nonumber\\ \displaystyle & & \displaystyle +\,{\it\rho}_{0}H\left(1-\frac{h}{H}\right)\left(-\frac{1}{2}g^{\prime }h+\frac{1}{3}\left\{{\rm\Delta}U\left(1-\frac{h}{H}\right)-\frac{3}{2}\sqrt{2g^{\prime }h}\right\}^{2}\right)-{\it\rho}_{0}H\frac{{\rm\Delta}U^{2}}{12}\nonumber\\ \displaystyle & & \displaystyle \left.-\,{\it\rho}_{0}H{\rm\Delta}U\left(1-\frac{h}{H}\right)^{2}\left\{\frac{1}{4}\left(1-\frac{h}{H}\right)^{2}{\rm\Delta}U-\left(1-\frac{h}{H}\right)\sqrt{2g^{\prime }h}\right\}\vphantom{\mathop{\sum }_{1}^{1}}\right].\end{eqnarray}$$
Figure 11 demonstrates that for positive shear values the pressure difference
${\rm\Delta}p$
increases. Recalling the definition of head loss and energy drop from (2.17) and (2.30), we find that this modified pressure difference term allows gravity currents with heights
$h/H>1/2$
to have a positive head loss and energy loss, suggesting a physically possible flow.
4. Comparison of model predictions with DNS results
In order to assess the validity of the above theoretical model predictions, we now present comparisons with two-dimensional DNS simulation results. These are obtained by our simulation 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 are provided in Nasr-Azadani & Meiburg (Reference Nasr-Azadani and Meiburg2011) and Nasr-Azadani, Hall & Meiburg (Reference Nasr-Azadani, Hall and Meiburg2013), so that a brief outline of the simulation approach will suffice here.
The fluid motion is described via the incompressible Navier–Stokes equations in the Boussinesq approximation,
$$\begin{eqnarray}\displaystyle & \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}^{\ast }=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial \boldsymbol{u}^{\ast }}{\partial t^{\ast }}+\boldsymbol{u}^{\ast }\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{u}^{\ast }=-\boldsymbol{{\rm\nabla}}p^{\ast }+\frac{1}{\mathit{Re}^{\ast }}{\rm\nabla}^{2}\boldsymbol{u}^{\ast }+{\it\rho}^{\ast }\boldsymbol{e}^{g}, & \displaystyle\end{eqnarray}$$
$\boldsymbol{e}^{g}$
and
$\mathit{Re}^{\ast }$
denote the unit vector acting in the direction of gravity and the Reynolds number 
$$\begin{eqnarray}\mathit{Re}^{\ast }=\frac{u_{b}H}{{\it\nu}}\end{eqnarray}$$
respectively. In (4.3),
$u_{b}$
,
$H$
and
${\it\nu}$
indicate the buoyancy velocity, the channel height and the kinematic viscosity of the fluid.
To describe the density field
${\it\rho}^{\ast }(\boldsymbol{x}^{\ast },t^{\ast })$
, we employ a continuum description of the density field
${\it\rho}^{\ast }$
and evolve it in an Eulerian manner by
$$\begin{eqnarray}\frac{\partial {\it\rho}^{\ast }}{\partial t^{\ast }}+\boldsymbol{u}^{\ast }\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\rho}^{\ast }=\frac{1}{\mathit{Re}^{\ast }\mathit{Sc}^{\ast }}{\rm\nabla}^{2}{\it\rho}^{\ast }.\end{eqnarray}$$
Here,
$\mathit{Sc}^{\ast }$
represents the Schmidt number associated with the diffusion of the density field
${\it\rho}^{\ast }$
. In most applications
$\mathit{Sc}^{\ast }\gg 1$
, although test simulations indicate that the precise value of
$\mathit{Sc}^{\ast }$
influences the results only weakly as long as
$\mathit{Sc}^{\ast }\geqslant O(1)$
(Härtel, Meiburg & Necker Reference Härtel, Meiburg and Necker2000). With this in mind, we employ
$\mathit{Sc}^{\ast }=6$
in all of our simulations.
The computational set-up employs a channel of size
$L_{x}^{\ast }\times L_{y}^{\ast }=4.5\times 1$
, with in/outflow boundaries in the horizontal direction, cf. figure 12. The Cartesian grid in the
$x^{\ast }$
- and
$y^{\ast }$
-directions is uniformly spaced with
${\rm\Delta}x^{\ast }=0.0056$
and
${\rm\Delta}y^{\ast }=0.0033$
. In the TL configuration, free-slip conditions are applied at the top and bottom walls. In the SL case, we impose the less restrictive condition
$\partial ^{2}u^{\ast }/\partial {y^{\ast }}^{2}=0$
at these walls, so that the slope of the velocity profile is free to adjust itself. At the outlet, all flow variables
$q^{\ast }$
are convected out of the domain via the outflow boundary condition
$$\begin{eqnarray}\frac{\partial q^{\ast }}{\partial t^{\ast }}+\bar{U}^{\ast }\frac{\partial q^{\ast }}{\partial x^{\ast }}=0,\end{eqnarray}$$
where
$\bar{U}^{\ast }$
represents the maximum
$u^{\ast }$
-velocity value in the domain. For the density field, we impose no-flux conditions at the top and bottom walls.
Figure 12. Initial set-up of the density field (a) and the horizontal velocity component (b), for the representative TL case of
$h^{\ast }=0.6$
,
${\rm\Delta}U^{\ast }=0.44$
and
$h_{1i}^{\ast }=0.5$
. The dashed lines in (b) indicate the interface locations.
Figure 13. The two-step initialization of the gravity current for
$h^{\ast }=0.6$
(dashed line),
${\rm\Delta}U^{\ast }=0.44$
,
$h_{1i}^{\ast }=0.5$
and
$\mathit{Re}^{\ast }=2600$
. The preliminary simulation run from the conditions shown in (a–d) ((a)
$t^{\ast }=0$
, (b)
$t^{\ast }=3$
, (c)
$t^{\ast }=6$
) is stopped at
$t^{\ast }=12$
(d). The
$u^{\ast }$
-velocity component in the domain
$x^{\ast }\geqslant 3$
is then set to
$u^{\ast }(x^{\ast }=3,y^{\ast })$
(vertical dotted line), while the
$v^{\ast }$
-velocity is set to zero in this region. The gravity current height is set to the model value
$h^{\ast }=0.6$
, so that
${\it\rho}^{\ast }$
is set to unity below this height and to zero elsewhere. (a–d) Initialization step I, (e) initialization step II.
The initial conditions for the simulation are generated in a two-step procedure, as explained in the following. As a first step, the gravity current front is placed at
$x^{\ast }=1.25$
, and the shape of its interface is approximated by an error function of the form
$$\begin{eqnarray}h_{s}^{\ast }(x^{\ast })=\left\{\begin{array}{@{}ll@{}}0,\quad & x^{\ast }<1.25,\\ \displaystyle h^{\ast }\text{erf}\left(\frac{x^{\ast }-1.25}{0.5}\right),\quad & x^{\ast }\geqslant 1.25.\end{array}\right.\end{eqnarray}$$
The initial ambient layer thickness values are chosen according to the model predictions, with a sharp jump at the location of the gravity current front. The horizontal ambient stream velocities upstream and downstream of the gravity current front are set to the model values. Inside the gravity current, the horizontal velocity is assumed to be zero. The vertical velocity is zero throughout the domain. The density
${\it\rho}^{\ast }$
is set to unity inside the gravity current, and to zero elsewhere.
Upon release, the velocity field quickly adjusts itself to account for the transition in the frontal region of the gravity current, cf. figure 13. Initially, Kelvin–Helmholtz vortices form downstream of the current front. Depending on the Reynolds number of the flow, these Kelvin–Helmholtz waves can be convectively or absolutely unstable (Huerre & Monkewitz Reference Huerre and Monkewitz1990), with larger
$\mathit{Re}^{\ast }$
-values favouring absolute instability. This limits the value of
$\mathit{Re}^{\ast }$
in our simulations, as we would like the vortices to be washed out of the domain, so that a nearly steady flow field is obtained that can be compared with the predictions of the steady-vorticity model. On the other hand, we would like to choose
$\mathit{Re}^{\ast }$
as large as possible, in order to minimize the effects of viscosity and diffusion, which are not contained in the vorticity model. These considerations lead us to employ intermediate
$\mathit{Re}^{\ast }$
-values ranging from 1800 to 3200, depending on the gravity current height. The influence of the precise
$\mathit{Re}^{\ast }$
-value on the results will be discussed in detail below.
In the second step of the initialization procedure, we now obtain improved initial conditions for the simulation in the following manner, cf. figure 13. After the above current reaches nearly steady-state conditions, we stop the simulation and determine a target location, where the gravity current interface first approaches a horizontal slope, cf. the vertical dotted line located at
$x^{\ast }=3$
in figure 13. The distribution of the horizontal velocity component across the channel height at this target location is then imposed everywhere downstream of this target location, and the vertical velocity is set to zero in this region. The gravity current height downstream of the target location is taken as the model value. Subsequently, the actual simulation is started with this set of initial conditions. We emphasize that the precise details of this two-step procedure for determining initial conditions do not affect the final nearly steady-state flow in any significant fashion (cf. Xu et al.
Reference Xu, Xue and Droegemeer1996). Table 1 summarizes the parameter combinations for which we have carried out DNS simulations.
Table 1. Parameter combinations of the simulations discussed in the present study, cf. the symbols in figures 3 and 9.
4.1. Influence of Reynolds number
The vorticity model assumes inviscid flow conditions and vanishing diffusion. The DNS simulations, on the other hand, employ finite Reynolds and Schmidt numbers, thereby giving rise to viscous dissipation and diffuse interfaces. Consequently, we do not expect the DNS results to match the model predictions exactly. This is confirmed by figure 14(a), which compares the quasisteady
$u^{\ast }$
-velocity profiles at a distance of three gravity current heights behind the front, for two different
$\mathit{Re}^{\ast }$
-values. However, we note that for larger Reynolds numbers the DNS velocity profiles match the model prediction increasingly well. Similarly, figure 14(b) compares the front location as a function of time for two different scenarios. While the model predicts a zero front velocity in the moving reference frame, the DNS results for
$\mathit{Re}^{\ast }=1100$
and
$\mathit{Re}^{\ast }=2600$
produce positive quasisteady front velocities of 0.018 and 0.013 respectively. Considering that the model prediction for the front velocity in the laboratory frame is 0.418, these small observed DNS front velocities nevertheless indicate good agreement between the model predictions and DNS results, with a relative error of only 4 % and 3 % respectively. Positive front velocities in the moving reference frame indicate that the gravity current moves somewhat slower than predicted by the vorticity model. We will now discuss the reasons for this reduced gravity current velocity.
Figure 14. Influence of Reynolds number on the DNS results, for the representative case T2 (
$h^{\ast }=0.6$
,
${\rm\Delta}U^{\ast }=0.44$
and
$h_{1i}^{\ast }=0.5$
). (a) Quasisteady horizontal velocity profile some distance downstream of the current front, for two different values of
$\mathit{Re}^{\ast }$
. As the Reynolds number increases, the DNS results approach the vorticity model prediction. (b) Front location as a function of time for two different
$\mathit{Re}^{\ast }$
-values. Front locations are recorded in the reference frame moving with the analytical value of the gravity current velocity
$U_{g}^{\ast }=0.418$
.
To explain this reduction in
$U_{g}^{\ast }$
as a result of viscous effects, we need to consider the resulting diffusion of momentum across shear layers. For case T2, momentum is transferred from the lower ambient layer (the fastest stream) to both the gravity current and the upper layer, both of which are slower. In this way, viscous effects will cause an overall head gain in the upper layer. Since the inflow shear remains the same, this head gain in the upper layer shifts point T2 in figure 3 slightly to the right, i.e. towards slightly lower gravity current velocities. This will be further investigated in § 4.2. Figure 14(b) indicates that for an increased Reynolds number
$u_{f}^{\ast }$
diminishes, so that the difference between the inviscid model prediction and the DNS simulation decreases. Figure 15 compares the depth-integrated height
$$\begin{eqnarray}\bar{h}^{\ast }(x^{\ast })=\int _{0}^{1}{\it\rho}^{\ast }(x^{\ast },y^{\ast })\,\text{d}y^{\ast }\end{eqnarray}$$
at the final simulation time
$t^{\ast }=60$
for several shear magnitudes.
Figure 15. Comparison of the steady-state depth-integrated current height
$\bar{h}^{\ast }$
(see (4.7)) produced by DNS simulations (solid lines) with analytical predictions made by our theoretical model (dotted lines). The values of
${\rm\Delta}U^{\ast }$
are 0, 0.21, 0.44 and 0.68, corresponding to gravity current heights of
$h^{\ast }=0.5$
, 0.55, 0.6 and 0.65 respectively. The predictions made by the circulation model are all for the zero-head-loss cases and are produced in the TL configuration with
$h_{1i}^{\ast }=h_{2i}^{\ast }=0.5$
.
Upon reaching steady-state conditions, we observe that the gravity current approaches the predicted thickness approximately two to three current heights downstream of the front. For increased shear, the current front is steeper, cf. figure 16. von Karman (Reference von Karman1940) and Benjamin (Reference Benjamin1968) employed inviscid assumptions to compute the tangent at the stagnation point for the classical gravity current without shear propagating over a slip wall. They found the angle at the stagnation point to be
${\it\theta}={\rm\pi}/3$
, which has been confirmed both numerically and experimentally. Applying inviscid conditions, the shape of the interface between the gravity current and the ambient fluid in the vicinity of the stagnation point is constrained by
$$\begin{eqnarray}[|\partial _{n}{\it\psi}|^{2}]=2y,\end{eqnarray}$$
with
$[~]$
,
$\partial _{n}$
and
${\it\psi}$
denoting the jump across the interface, the partial derivative normal to the interface and the streamfunction respectively (cf. Xu & Moncrieff Reference Xu and Moncrieff1994). While a non-zero shear in the free stream layer does modify the shape of the interface away from the stagnation point (point
$\text{O}$
in figure 16), the slope right at the stagnation point (shown by the solid line in figure 16) does not appear to depend strongly on the value of the shear in the free stream (or even inside the gravity current) and is very close to the classical value
${\it\theta}={\rm\pi}/3$
(cf. the discussions by Xu Reference Xu1992 and Xu & Moncrieff Reference Xu and Moncrieff1994). While positive shear values steepen the interface away from the stagnation point, for the present set of simulations the slope angle remains below
${\rm\pi}/3$
everywhere. As we will see in § 4.5, for very large shear values, when one of the oncoming free stream flows exhibits a negative velocity, an inflection point on the gravity current interface emerges such that the slope can locally increase beyond
${\it\theta}={\rm\pi}/3$
(Xu & Moncrieff Reference Xu and Moncrieff1994).
For the representative single-layer case S1, figure 17 compares the DNS results with the vorticity model prediction some distance downstream of the gravity current front. Good agreement is observed, especially for the shape of the velocity profile in the ambient stream.
Figure 17. Vertical distribution of the horizontal velocity component
$u^{\ast }$
and the density difference
${\it\rho}^{\ast }$
, plotted three gravity current heights downstream of the gravity current front for case S1. This case corresponds to the SL configuration with
$h^{\ast }=0.6$
and
${\rm\Delta}U^{\ast }=0.42$
. The circles (
$u^{\ast }$
) and triangles (
${\it\rho}^{\ast }$
) represent the DNS results. These are in very good agreement with the predictions of the circulation model, indicated by the dashed and solid lines (see table 1).
Figure 18. Head loss along streamlines as a function of the
$y^{\ast }$
-coordinate of the streamline at the inflow boundary, computed from DNS simulations. The cases T0, T2 and S1 are shown, cf. table 1. A discussion of the results is provided in the text.
4.2. Head loss analysis
We now proceed to analyse the head loss in the DNS simulations. Towards this end, we employ the steady-state DNS solution to evaluate Bernoulli’s equation along all streamlines originating at the inlet (
$x^{\ast }=0$
) and ending at the outlet (
$x^{\ast }=4.5$
). Figure 18 shows the head loss as a function of the
$y^{\ast }$
-coordinate at the inflow boundary, for three different flows. Case T0 refers to a gravity current occupying half the channel height in the absence of shear, cf. table 1, while S1 and T2 correspond to currents of
$h^{\ast }=0.6$
with single-layer and two-layer shear respectively. All of the flows experience positive head loss along the streamlines near the interface between the ambient and the gravity current, since momentum diffuses from the ambient to the gravity current.
For both S1 and T0, the head loss drops to zero as we approach the top wall. For T0, this result is consistent with Benjamin (Reference Benjamin1968) finding that gravity currents with half the channel height conserve energy. For T2, we observe a head loss in the lower faster ambient stream, and a corresponding head gain in the upper slower ambient stream. This is again consistent with the viscous diffusion of momentum from the lower into the upper ambient stream.
4.3. Energy analysis
A DNS-based analysis of the various terms in the kinetic energy equation provides additional useful information. In the Boussinesq approximation, the steady-state kinetic energy equation can be written in the form
$$\begin{eqnarray}\boldsymbol{u}^{\ast }\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}k^{\ast }={\it\rho}^{\ast }\boldsymbol{e}_{g}\boldsymbol{\cdot }\boldsymbol{u}^{\ast }-\boldsymbol{u}^{\ast }\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}P^{\ast }+\frac{\partial }{\partial x_{j}^{\ast }}\left(\frac{2}{\mathit{Re}^{\ast }}u_{i}^{\ast }\unicode[STIX]{x1D61A}_{ij}^{\ast }\right)-\frac{2}{\mathit{Re}^{\ast }}\unicode[STIX]{x1D61A}_{ij}^{\ast }\unicode[STIX]{x1D61A}_{ij}^{\ast }.\end{eqnarray}$$
Here,
$k^{\ast }$
and
$\unicode[STIX]{x1D61A}_{ij}^{\ast }$
denote the kinetic energy
$u_{i}^{\ast }u_{i}^{\ast }/2$
and the strain-rate tensor
$(\partial u_{i}^{\ast }/\partial x_{j}^{\ast }+\partial u_{j}^{\ast }/\partial x_{i}^{\ast })/2$
respectively. It should be noted that
$P^{\ast }$
excludes the hydrostatic pressure due to the ambient fluid, i.e.
$P^{\ast }=(p-{\it\rho}_{0}g(H-y))/{\it\rho}_{0}u_{b}^{2}$
. For a given control volume, application of the Gauss divergence theorem converts (4.9) into
$$\begin{eqnarray}\oint \left(P^{\ast }+\frac{1}{2}u_{i}^{\ast }u_{i}^{\ast }+{\it\rho}^{\ast }y^{\ast }\right)\mathop{u}_{j}^{\ast }n_{j}\,\text{d}{\it\Gamma}^{\ast }=-\underbrace{\oint \frac{-2}{\mathit{Re}^{\ast }}(\mathop{u}_{i}^{\ast }\unicode[STIX]{x1D61A}_{ij}^{\ast }n_{j})\,\text{d}{\it\Gamma}^{\ast }}_{\dot{{\it\Psi}}^{\ast }}-\underbrace{\iint \frac{2}{\mathit{Re}^{\ast }}\unicode[STIX]{x1D61A}_{ij}^{\ast }\unicode[STIX]{x1D61A}_{ij}^{\ast }\,\text{d}A^{\ast }}_{\dot{{\it\Phi}}^{\ast }}.\end{eqnarray}$$
The left-hand side surface integral refers to the net flux of energy across the boundary
${\it\Gamma}^{\ast }$
into and out of the control volume
$A^{\ast }$
. The two terms on the right-hand side of (4.10),
$\dot{{\it\Psi}}^{\ast }$
and
$\dot{{\it\Phi}}^{\ast }$
, denote the work performed by shear forces on the surface of the control volume, i.e. the viscous diffusion of momentum and the dissipation of kinetic energy into heat.
Figure 19 presents the viscous dissipation field for configurations T2 and S1. Both cases demonstrate comparable dissipation rates along the interface between the gravity current and the ambient. Case TL experiences additional significant dissipation along the interface between the two ambient streams, as a result of the local velocity jump
${\rm\Delta}U^{\ast }$
.
Figure 19. Viscous dissipation fields for the cases T2 (a) and S1 (b). Strong dissipation occurs along all interfaces with shear. White corresponds to zero and black corresponds to 0.1. Dashed lines identify the analytical gravity current height
$h^{\ast }=0.6$
.
For the control volume around the ambient stream(s) (
$\text{AO}\text{A}^{\prime }\text{C}^{\prime }\text{C}\text{A}$
in figures 1 and 2), we evaluate the energy influx
${\dot{E}}_{in}^{\ast }$
(integrated at
$\text{AC}$
) and outflux
${\dot{E}}_{out}^{\ast }$
(integrated at
$\text{A}^{\prime }\text{C}^{\prime }$
). We furthermore compute the integral over the viscous dissipation fields separately for the gravity current (
$\dot{{\it\Phi}}_{1}^{\ast }$
) and the ambient (
$\dot{{\it\Phi}}_{0}^{\ast }$
),
$$\begin{eqnarray}\displaystyle & \displaystyle \dot{{\it\Phi}}_{0}^{\ast }=\iint \frac{2}{\mathit{Re}^{\ast }}\unicode[STIX]{x1D61A}_{ij}^{\ast }\unicode[STIX]{x1D61A}_{ij}^{\ast }(1-{\it\rho}^{\ast })\,\text{d}A^{\ast }, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \dot{{\it\Phi}}_{1}^{\ast }=\iint \frac{2}{\mathit{Re}^{\ast }}\unicode[STIX]{x1D61A}_{ij}^{\ast }\unicode[STIX]{x1D61A}_{ij}^{\ast }{\it\rho}^{\ast }\,\text{d}A^{\ast }. & \displaystyle\end{eqnarray}$$
Table 2 reports the various energy budget components for the cases T0, T2 and S1. We compute the energy influx
${\dot{E}}_{in}^{\ast }$
and energy outflux
${\dot{E}}_{out}^{\ast }$
in the free stream layer via
$$\begin{eqnarray}\displaystyle & \displaystyle {\dot{E}}_{in}^{\ast }=\int _{0}^{1}u^{\ast }\left(P^{\ast }+\frac{1}{2}u_{i}^{\ast }u_{i}^{\ast }\right)\,\text{d}y^{\ast }\quad \text{at}~x^{\ast }=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\dot{E}}_{out}^{\ast }=\int _{0}^{1}(1-{\it\rho}^{\ast })u^{\ast }\left(P^{\ast }+\frac{1}{2}u_{i}^{\ast }u_{i}^{\ast }+{\it\rho}^{\ast }y^{\ast }\right)\,\text{d}y^{\ast }\quad \text{at}~x^{\ast }=L_{x}^{\ast }. & \displaystyle\end{eqnarray}$$
Table 2. Different energy components computed for three cases:
$h^{\ast }=0.6$
(TL),
$h^{\ast }=0.6$
(SL) and
$h^{\ast }=0.5$
.
For a gravity current height
$h^{\ast }=0.6$
, approximately 20 % (T2) and 26 % (S1) of the energy influx is dissipated by viscosity, whereas 9 % (T2) and 13 % (S1) is used to perform work on the gravity current. The diffusion of momentum from the ambient to the gravity current causes a recirculation region to form inside the gravity current, cf. figure 19. This recirculation region lifts incoming gravity current fluid from near the bottom wall and alters its velocity, thereby changing the energy of this fluid
${\dot{E}}_{g}^{\ast }$
. To obtain
${\dot{E}}_{g}^{\ast }$
, we compute the following integral at the outlet region:
$$\begin{eqnarray}{\dot{E}}_{g}^{\ast }=\int _{0}^{1}{\it\rho}^{\ast }u^{\ast }\left(P^{\ast }+{\textstyle \frac{1}{2}}u_{i}^{\ast }u_{i}^{\ast }+{\it\rho}^{\ast }y^{\ast }\right)\,\text{d}y^{\ast }\quad \text{at}~x^{\ast }=L_{x}^{\ast }.\end{eqnarray}$$
In steady-state conditions, this energy change
${\dot{E}}_{g}^{\ast }$
and the viscous dissipation inside the gravity current
$\dot{{\it\Phi}}_{1}^{\ast }$
are balanced by the work
$\dot{{\it\Psi}}^{\ast }$
performed by the ambient on the gravity current. We apply this concept to define the relative error
${\it\epsilon}$
for each simulation, i.e.
$$\begin{eqnarray}{\it\epsilon}=\left|\frac{\dot{{\it\Psi}}^{\ast }-(\dot{{\it\Phi}}_{1}^{\ast }+{\dot{E}}_{g}^{\ast })}{{\dot{E}}_{in}^{\ast }}\right|.\end{eqnarray}$$
The values reported in table 2 indicate that the components in the energy budget balance each other with excellent accuracy.
4.4. Non-zero head loss
In the following, we investigate currents that are initialized with conditions far away from the zero-head-loss contour in figure 3. Towards this end, we will focus on the TL configuration with
$h_{1i}^{\ast }=0.5$
. Since the formation of a quasisteady state near the current front can take a significant amount of time, we employ longer computational domains than before.
As a representative case with initial conditions corresponding to a positive head loss, we choose a shear parameter value of
${\rm\Delta}U^{\ast }=0.44$
and a gravity current height of
$h^{\ast }=0.45$
. In the phase-space diagram of figure 3, these initial conditions correspond to point T1, with a gravity current velocity of
$U_{g}^{\ast }=0.475$
and a positive head loss of
${\it\delta}_{CC^{\prime }}^{\ast }=0.0369$
. Figure 20 depicts the temporal evolution of the current in the DNS simulation, in a reference frame moving with
$U_{g}^{\ast }=0.475$
. After the initial formation of Kelvin–Helmholtz instabilities, the flow undergoes a pronounced transition to a quasisteady state, characterized by a thicker current of
$h^{\ast }\approx 0.6$
(shown by a dashed line in figure 20) and a constant front velocity
$u_{f}^{\ast }\approx 0.068$
, cf. figure 21. This indicates that initial conditions corresponding to a large positive head loss value do not persist, and instead a current develops with a head loss near zero for the imposed shear value. This is confirmed by the fact that the DNS value of
$U_{g}^{\ast ^{\prime }}=0.407$
for the front velocity in the laboratory frame is close to the theoretical value of
$U_{g}^{\ast }=0.418$
for the zero-head-loss case.
Figure 20. Temporal evolution of the current initialized with
$h^{\ast }=0.45$
and
${\rm\Delta}U^{\ast }=0.44$
. This corresponds to a positive head loss of
${\it\delta}_{CC^{\prime }}^{\ast }=0.037$
. The current reaches a quasisteady state with
$h^{\ast }\approx 0.6$
(dashed line) and a constant front velocity
$u_{f}^{\ast }\approx 0.068$
: (a)
$t^{\ast }=0$
; (b) 6; (c) 20; (d) 40; (e) 60.
Figure 21. Temporal evolution of the front location for positive and negative head losses. To explore the positive-head-loss scenario, the flow was initialized with
$h^{\ast }=0.45$
and
${\rm\Delta}U^{\ast }=0.44$
. The current reaches a quasisteady state with
$h^{\ast }\approx 0.6$
, which closely corresponds to vanishing head loss. The quasisteady computational front velocity is approximately
$u_{f}^{\ast }\approx 0.070$
, which is close to the analytical value of
$0.475-0.418=0.057$
for zero head loss. To investigate negative-head-loss scenarios, we initiated a current with
$h^{\ast }=0.75$
. Here, a quasisteady computational front velocity
$u_{f}^{\ast }\approx 0.102$
evolves, which is again close to the analytical value of
$0.302-0.418=0.116$
for vanishing head loss.
It is equally interesting to investigate initial conditions corresponding to a strongly negative head loss. We choose
$h^{\ast }=0.75$
and
${\rm\Delta}U^{\ast }=0.44$
, with a resulting head loss of
${\it\delta}_{CC^{\prime }}^{\ast }=-0.107$
and gravity current velocity of
$U_{g}^{\ast }=0.302$
, indicated as T3 in figure 3. Figure 22 presents the temporal evolution of the flow. After a pronounced thinning of the current from the initial conditions, the current starts to thicken again until it reaches a quasisteady state with a current height of
$h^{\ast }\approx 0.6$
and a front velocity of
$u_{f}^{\ast }\approx -0.103$
, cf. figure 21. Similarly to the positive-head-loss scenario, this suggests that initial conditions corresponding to a strongly negative head loss cannot be maintained, and instead a quasisteady current evolves with conditions close to zero head loss.
Figure 22. Temporal evolution of the current initialized with
$h^{\ast }=0.75$
and
${\rm\Delta}U^{\ast }=0.44$
, corresponding to a negative head loss. The current height approaches
$h^{\ast }\approx 0.6$
, which is close to the vanishing-head-loss state: (a)
$t^{\ast }=0$
; (b) 6; (c) 20; (d) 40; (e) 60.
4.5. Mixed inflow–outflow in the moving reference frame
We now discuss a case in which the ambient flow displays mixed in- and outflow in the reference frame moving with the current front, due to high shear magnitude. Such high-shear scenarios can occur under certain atmospheric conditions, for example during the interaction of environmental shear with a cool thunderstorm outflow when producing long-lived squall lines (Xu Reference Xu1992; Xu & Moncrieff Reference Xu and Moncrieff1994; Xue Reference Xue2002; Bryan & Rotunno Reference Bryan and Rotunno2014a
,Reference Bryan and Rotunno
b
). We limit the discussion to the SL configuration, although a corresponding analysis can be applied to the TL case as well. We choose the representative values of
${\rm\Delta}U^{\ast }=2.256$
and
$h^{\ast }=0.9$
, for which (2.15) returns a gravity current velocity equal to
$U_{g}^{\ast }=0.12$
.
Figure 23. Evolution of a gravity current with mixed in- and outflow at the left boundary, case SL. The height
$h^{\ast }=0.9$
and shear magnitude
${\rm\Delta}U^{\ast }=2.256$
result in a gravity current velocity of
$U_{g}^{\ast }=0.123$
. The strong shear causes a recirculation region in the ambient stream, which deforms the gravity current front: (a)
$t^{\ast }=0$
; (b) 2; (c) 6; (d) 60.
As shown in figure 24, the linear velocity profile at the left boundary displays both in- and outflow regions, although the depth-averaged flow still is from left to right. A depth-integrated energy analysis along the lines discussed in § 2.3 indicates zero energy loss for this parameter combination, suggesting that such a case with mixed in- and outflow conditions should be physically feasible.
Figure 23 depicts the corresponding DNS simulation results. After an initial transient period, the current approaches the predicted height of
$h^{\ast }=0.9$
with a stationary nose, i.e.
$u_{f}^{\ast }\approx 0$
. The frontal region of the gravity current displays a very abrupt increase in height, which is linked to the existence of a recirculation region in the ambient current, as can be seen in the steady-state streamline pattern of figure 24. Moreover, unlike the cases shown in figure 16, an inflection point emerges along the interface which allows for interfacial slopes larger than
${\rm\pi}/3$
away from the stagnation point (where the tangent slope remains
${\rm\pi}/3$
). Corresponding analytical and numerical findings have been reported in previous studies (Xue Reference Xue2002; Bryan & Rotunno Reference Bryan and Rotunno2014a
,Reference Bryan and Rotunno
b
).
Figure 24. Steady-state streamlines for the SL case shown in figure 23. The grey region indicates the gravity current with height
$h^{\ast }=0.9$
(shown by the dashed line). The recirculation region in the ambient stream is visible.
4.6. High Reynolds numbers and influence of mixing
Up to now, we have limited ourselves to moderate-Reynolds-number DNS simulations. Here, Kelvin–Helmholtz instabilities at the interface between the gravity current and the ambient are convectively unstable, so that they are swept out of the control volume and a nearly steady-state flow field is established for long times (see table 1 and figure 13). This enabled us to compare steady-state DNS velocity and density profiles with analytical model predictions. In natural settings, however, gravity currents frequently occur at Reynolds numbers orders of magnitude higher than the values discussed in table 1. Under such circumstances, the layer separating the gravity current from the ambient is strongly three-dimensional, unsteady and turbulent, and it exhibits significant mixing. These effects are difficult to capture in simplified models such as the present one, or similar ones by earlier authors (Benjamin Reference Benjamin1968; Xu Reference Xu1992; Shin et al.
Reference Shin, Dalziel and Linden2004), and frequently have to be addressed via numerical simulations (e.g. Härtel et al.
Reference Härtel, Meiburg and Necker2000; Cantero et al.
Reference Cantero, Lee, Balachandar and Garcia2007; Nasr-Azadani & Meiburg Reference Nasr-Azadani and Meiburg2014). In the following, we will employ higher-Reynolds-number simulations in order to explore some of the differences due to unsteadiness and mixing. Towards this end, we restrict ourselves to case S1 in the SL configuration, with a gravity current height of
$h^{\ast }=0.6$
(see table 1).
We discuss five DNS simulations for Reynolds numbers of 3000, 3500, 4000, 5000 and
$10\,000$
, where a nearly steady flow field no longer develops. In order to minimize any influence of the outflow boundary on the flow, we choose a larger computational domain of size
$L_{x}^{\ast }\times L_{y}^{\ast }=9\times 1$
, cf. figure 2.
Figure 25 displays the density fields for different
$\mathit{Re}^{\ast }$
-values, after they reach statistically steady conditions. As
$\mathit{Re}^{\ast }$
increases, the effect of mixing at the interface becomes more pronounced. For
$\mathit{Re}^{\ast }=3000$
, despite the existence of Kelvin–Helmholtz instabilities far downstream, the upstream sections of the current still reach the predicted height, as indicated by the horizontal solid line. As
$\mathit{Re}^{\ast }$
further increases, the onset of Kelvin–Helmholtz instabilities moves closer to the current front, and eventually prevents even the sections next to the current front from reaching the predicted height. For comparison, the numerically evaluated gravity current heights
$\bar{h}^{\ast }$
(4.7) are shown by dashed white lines.
Figure 25. Density fields at
$t^{\ast }=50$
for configuration S1 and different Reynolds numbers
$\mathit{Re}^{\ast }$
: (a) 3000; (b) 3500; (c) 4000; (d) 5000; (e) 10 000. As the Reynolds number increases, Kelvin–Helmholtz instabilities and the associated unsteady mixing prevent the gravity current from approaching the height predicted by the vortex sheet model (shown by the horizontal solid lines). The instantaneous current height
$\bar{h}^{\ast }$
(4.7) is shown by the dashed white line.
For different
$\mathit{Re}^{\ast }$
-values, figure 26 compares the temporally averaged current height
$\langle \bar{h}^{\ast }(x^{\ast })\rangle$
, which is defined as
$$\begin{eqnarray}\langle \bar{h}^{\ast }(x^{\ast })\rangle =\frac{1}{T_{1}^{\ast }-T_{0}^{\ast }}\int _{T_{0}^{\ast }}^{T_{1}^{\ast }}\bar{h}^{\ast }({\it\tau}^{\ast };x^{\ast })\,\text{d}{\it\tau}^{\ast }.\end{eqnarray}$$
Here,
$T_{0}^{\ast }$
and
$T_{1}^{\ast }$
denote the beginning and end times of the interval over which the average is taken, which are set to 30 and 60 respectively. The figure indicates that away from the current front
$\langle \bar{h}^{\ast }(x^{\ast })\rangle$
approaches a value that is significantly lower than the predicted height of
$h^{\ast }=0.6$
and does not appear to depend strongly on the Reynolds number. This drop in height is associated with the vortical motion due to the Kelvin–Helmholtz instabilities, which alter the nature of the flow in the interfacial region.
The above reduction in the current height has also been observed by other authors for numerically simulated gravity currents with and without shear (Xu et al.
Reference Xu, Xue and Droegemeer1996; Bryan & Rotunno Reference Bryan and Rotunno2014b
). Furthermore, Benjamin (Reference Benjamin1968) suggests that for the classical case of a gravity current without shear, heights
$h/H\geqslant 0.347$
should not be expected, due to turbulent dissipation and the resulting diffuse interface between the current and the ambient. Klemp, Rotunno & Skamarock (Reference Klemp, Rotunno and Skamarock1994) employ shallow-water relations along with hydraulic constraints to argue that for the flow conditions to stay subcritical behind the gravity current front (in the reference frame moving with the gravity current head), the gravity current height must always satisfy the
$h/H\leqslant 0.347$
condition.
Xu (Reference Xu1992) considers the SL configuration and extends the approach by Benjamin (Reference Benjamin1968) to gravity currents propagating into shear. Assuming that energy is conserved, he derives in a first step a relationship between the current height and the ambient shear that is identical to the solution obtained in the present investigation in the limit of vanishing energy loss. Based on the observation that numerical simulations usually yield currents of lower height than predicted by the above idealized model, Xu (Reference Xu1992) then develops an augmented model which is designed to allow for a head loss that varies across the streamlines. In this way, this augmented model offers an additional degree of flexibility compared with Benjamin’s model and the present one. However, it needs to quantify this height-dependent head loss from empirical considerations before it can obtain the front velocity.
Bryan & Rotunno (Reference Bryan and Rotunno2014b
) employ Xu’s (Reference Xu1992) modelling approach and study the cases with the maximum energy dissipation possible. Bryan & Rotunno (Reference Bryan and Rotunno2014b
) demonstrate that the gravity current velocity obtained from the inviscid model with imposed energy dissipation agrees well with three-dimensional simulation results (see figure 15 in Bryan & Rotunno Reference Bryan and Rotunno2014b
). However, the gravity current height does not track the model predictions closely. We note that in the investigations of Xu et al. (Reference Xu, Xue and Droegemeer1996) and Bryan & Rotunno (Reference Bryan and Rotunno2014b
), the influence of the
$\mathit{Re}^{\ast }$
-value could not be analysed, since in their simulations the physical viscosity was set to zero, so that they contained only numerical dissipation.
In order to further explore the effect of high Reynolds numbers on the gravity current height, we carry out several additional simulations that are initiated with conditions on the zero-head-loss line for different values of
$h^{\ast }$
and
${\rm\Delta}U^{\ast }$
in the SL configuration. Positive, zero and negative shear magnitudes are considered, and for each case we evaluate the gravity current velocity
$U_{g}^{\ast }$
via (2.15). The simulation parameters are given in the first three columns of table 3.
Figure 27. (a) Temporally averaged gravity current height
$\langle \bar{h}^{\ast }\rangle$
at
$x^{\ast }=7$
for the cases discussed in table 3. (b) Gravity current velocity computed for the same cases. The solid lines indicate the model predictions for zero head loss, while the dashed lines show the model predictions for maximum head loss. The black symbols represent DNS results.
Table 3. Parameters used for DNS simulations chosen on the zero-head-loss line. All simulations are conducted in the SL configuration and with the Reynolds number (formed with the gravity current height
$h$
) equal to 2300. With the presence of Kelvin–Helmholtz instabilities, the averaged current height obtained from the DNS runs
$\langle \bar{h}^{\ast }\rangle$
seems to follow a reduced height corresponding to the case with maximum energy loss (head loss),
$h_{ml}^{\ast }$
.
We conduct each DNS simulation with a Reynolds number (formed with the gravity current height
$h$
) of 2300. For this
$\mathit{Re}^{\ast }$
-value, all simulations exhibit absolutely unstable Kelvin–Helmholtz instabilities in the interfacial region. Upon reaching statistically steady conditions, we record the gravity current height
$\bar{h}^{\ast }$
downstream of the current front at
$x^{\ast }=7$
in order to obtain
$\langle \bar{h}^{\ast }\rangle$
from (4.17) at this location (fifth column of table 3).
Figure 27(a) compares the temporally averaged gravity current height
$\langle \bar{h}^{\ast }\rangle$
recorded from the DNS simulations (black symbols) with the analytical predictions for zero head loss (solid line). The error bars indicate the maximum and minimum gravity current heights recorded in the time interval
$30\leqslant t^{\ast }\leqslant 60$
. The DNS values clearly do not follow the zero-head-loss predictions. For comparison, the analytical predictions for maximum head loss are shown by the dashed line. This curve exhibits good agreement with the DNS data, which suggests that for high
$\mathit{Re}^{\ast }$
-values the gravity current height adjusts itself so as to produce maximum head loss. A corresponding observation was reported by Xu (Reference Xu1992) based on the maximum energy dissipation function of his model. On the other hand, the gravity current velocity, shown in figure 27(b), does not seem to be strongly affected by the effects of unsteadiness and mixing.
For three different Reynolds numbers, figure 28 compares the temporal evolution of the front location
$x_{f}^{\ast }$
in the reference frame moving with the velocity predicted by the model. As
$\mathit{Re}^{\ast }$
increases, the front velocity in this moving reference frame increases, indicating deteriorating agreement with the model prediction. Figure 29 displays temporally averaged profiles of the horizontal velocity
$\langle u^{\ast }\rangle$
and the density field
$\langle {\it\rho}^{\ast }\rangle$
downstream of the gravity current front at
$x^{\ast }=6$
, for three different
$\mathit{Re}^{\ast }$
-values. Compared with the earlier steady case in figure 17, where the velocity and density profiles approach the values in the ambient layer rather sharply, the present unsteady flows demonstrate a much smoother transition from current to ambient. This behaviour is investigated by Borden & Meiburg (Reference Borden and Meiburg2013a
), who develop a diffuse mixing layer model with smooth velocity and density profiles connecting the gravity current and the ambient (see figure 8 in Borden & Meiburg Reference Borden and Meiburg2013a
). By employing profiles based on the DNS data, they obtain better agreement between model predictions and DNS results at high
$\mathit{Re}^{\ast }$
-values.
Figure 28. Influence of the Reynolds number on the front velocity. As the Reynolds number increases, the Kelvin–Helmholtz instabilities and mixing at the interface cause the agreement with the model predictions to deteriorate. Front velocities are recorded in the reference frame moving with the analytical value of the gravity current velocity
$U_{g}^{\ast }=0.404$
.
Figure 29. Temporally averaged horizontal velocity and density profiles across the channel height at
$x^{\ast }=6$
, for the gravity currents shown in figure 25. As the Reynolds number increases, strong Kelvin–Helmholtz instabilities smear out the averaged velocity and density profiles, leading to deteriorating agreement with predictions by the vortex sheet model. The horizontal solid line indicates the predicted gravity current height
$h^{\ast }=0.6$
. The temporal averaging extends from
$t^{\ast }=30$
to
$t^{\ast }=60$
.
5. Concluding remarks
We have introduced an analytical vorticity-based model for steady-state inviscid Boussinesq gravity currents in sheared ambients. The model is based on the conservation of mass and horizontal and vertical momentum, and it does not require any empirical closure assumptions. Among other things, the model predicts the existence of gravity currents that are thicker than half the channel height, due to the energy input provided by the shear in the ambient stream. The feasibility of such currents is confirmed by moderate-Reynolds-number DNS results, as well as by an analysis of the energy loss in the flow. The general agreement between the model predictions and the DNS results is excellent. Specifically, the DNS results demonstrate the formation of gravity currents characterized by a relatively small amount of energy loss predominantly caused by viscous and diffusion effects on the flow field. Differences between the model predictions and the DNS results can be understood via an analysis of the head loss/head gain as a result of the diffusion of momentum across sheared interfaces.
If the DNS simulations prescribe initial conditions corresponding to a large energy loss, or to an energy gain, the flow undergoes a transient readjustment until it reaches a quasisteady state that corresponds to a small energy loss.
We further investigated the influence of the ambient layer thicknesses and shear magnitude on the range of physically feasible solutions. Specifically, the results demonstrate that in the TL configuration the thickness of the lower ambient layer plays a crucial role with regard to the range of shear magnitudes for which gravity currents can exist with thicknesses larger than half the channel height.
Consistent with the findings of earlier authors (Benjamin Reference Benjamin1968; Xu Reference Xu1992; Borden & Meiburg Reference Borden and Meiburg2013a ), we find that for larger Reynolds numbers the current exhibits a qualitatively different behaviour, characterized by Kelvin–Helmholtz vortices as well as strong mixing and dissipation along the interface between the current and the ambient. While the current velocity does not change drastically, the temporally averaged velocity and density profiles of such flows vary smoothly between the current and the ambient, resulting in a current height that is substantially lower than the one predicted by the vortex sheet model. Our DNS simulations suggest that for a given shear, such gravity currents acquire the height that corresponds to the maximum energy dissipation.
Acknowledgement
We gratefully acknowledge support through NSF grant CBET-1335148.
