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

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

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

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

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

One final equation is obtained from geometrical considerations:

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:

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

where $U_{i}(y)$ has the linear form

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

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),

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,

In the following, we limit ourselves to cases where the gravity current travels from right to left in the laboratory frame, so that $U_{g}\geqslant 0$ . This imposes an upper bound on the shear magnitude in the form

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

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,

in conjunction with the hydrostatic pressure relations far up- and downstream of the gravity current front,

We thus obtain for the pressure drop along the top wall This allows us to calculate the head loss ${\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

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

We recall that for the SL configuration, the velocities in (2.23) are given as

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

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

Using the hydrostatic pressure relations at the inlet and outlet, we arrive at

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

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.

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

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

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).

Figure 16.

Comparison of the steady-state depth-integrated current height $\bar{h}^{\ast }$ (see (4.7)) produced by DNS simulations shown in figure 15. The shear magnitude does not influence the slope of the tangent line right at the stagnation point $\text{O}$ .

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

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

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 }$ ),

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

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:

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. 

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

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.

Figure 26.

Temporally averaged current height $\langle \bar{h}^{\ast }\rangle$ (see (4.17)) for the gravity currents shown in figure 25. As the Reynolds number increases, Kelvin–Helmholtz instabilities and the associated mixing reduce the current height behind the frontal 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$ .