2.1 Partially blocked flows

Figure 4 illustrates the temporal evolution of the density field for simulation S1. Transient, convectively unstable Kelvin–Helmholtz waves along the gravity current interface are soon washed out of the control volume, and a quasisteady state evolves. We remark that for some simulations involving larger values of $Re$ (not shown) the KH instabilities remained within the computational domain, so that a strictly quasisteady state did not emerge.

In order to compute the gravity current thickness, we mark the gravity current fluid initially by a passive dye of concentration $c_{d}^{\ast }=1$ . This dye has the same diffusivity as the density field. The local dimensionless height $d^{\ast }$ of the gravity current is then obtained as

(2.8) $$\begin{eqnarray}d^{\ast }(x^{\ast },t^{\ast })=\int _{0}^{1}c_{d}^{\ast }(x^{\ast },y^{\ast },t^{\ast })\,\text{d}y^{\ast }.\end{eqnarray}$$

The gravity current front location $x_{f}^{\ast }$ is defined as the leftmost point at which $d^{\ast }$ exceeds 0.01. The local thicknesses $d_{1}^{\ast }$ and $d_{2}^{\ast }$ of the ambient layers can be obtained analogously from corresponding passive dye fields. We furthermore identify the right and left edges of the ambient interfacial wave, $x_{r}^{\ast }$ and $x_{l}^{\ast }$ , in the following fashion: to obtain $x_{r}^{\ast }$ as a function of time, we sweep the ambient interface from the current tip towards the outflow boundary and evaluate $x_{r}^{\ast }$ as the first point at which $|\unicode[STIX]{x2202}(d^{\ast }+d_{1}^{\ast })/\unicode[STIX]{x2202}x^{\ast }|<0.01$ . $x_{l}^{\ast }$ is taken as the location of the leftmost point of the ambient interface where $d_{1}^{\ast }-h_{1i}^{\ast }>0.01$ . For those cases in which the flow gives rise to an internal bore propagating ahead of the expansion fan, this procedure yields the front location $x_{b}^{\ast }$ of the bore. In the presence of a bore, the location $x_{l}^{\ast }$ of the left edge of the expansion wave is then defined as the first point to the right of the bore at which $\unicode[STIX]{x2202}d_{1}^{\ast }/\unicode[STIX]{x2202}x^{\ast }|_{x^{\ast }}-\unicode[STIX]{x2202}d_{1}^{\ast }/\unicode[STIX]{x2202}x^{\ast }|_{x^{\ast }-\unicode[STIX]{x0394}x^{\ast }}<0$ . This definition is based on the rapid fall of the ambient interface slope as it transitions from the bore to the expansion fan region. In some cases this decrease occurs only very close to the gravity current tip, so that the disturbance takes the form of an expansion fan without a leading bore.

Within the moving reference frame of the simulation, figure 5(a) depicts the displacement with time of the gravity current tip $x_{f}^{\ast }$ , relative to its initial location $x_{f,i}^{\ast }$ , for all simulations of table 1. In all cases the front velocity $\text{d}x_{f}^{\ast }/\text{d}t^{\ast }$ initially varies with time, but eventually settles down to a constant value $u_{f}^{\ast }$ . The front velocity in the laboratory frame is then obtained as

(2.9) $$\begin{eqnarray}U_{g,lab}^{\ast }=U_{g}^{\ast }-u_{f}^{\ast }.\end{eqnarray}$$

Figure 5.

DNS results for the temporal evolution of (a) the displacement of the gravity current tip $x_{f}^{\ast }$ with respect to its initial location $x_{f,i}^{\ast }$ , (b) the displacement of the expansion wave’s left edge $x_{l}^{\ast }$ relative to the gravity current front location $x_{f}^{\ast }$ and (c) the distance of the right edge $x_{r}^{\ast }$ of the expansion fan from the current tip, evaluated for the simulations of table 1. The straight line segments indicate the quasisteady velocities of the corresponding features in the moving reference frame, which have been obtained by linear fits of the DNS results between $t^{\ast }=20$ and $t^{\ast }=40$ , when the flow properties have become nearly independent of time. Note that the quasisteady velocities vanish for simulations S1, S4 and S5 in panel (b), and for all simulations in panel (c).

Figure 5(b) displays the displacement of the left edge of the expansion wave $x_{l}^{\ast }$ with respect to the gravity current front location $x_{f}^{\ast }$ as a function of time, for all simulations. After a brief transitional phase, $x_{l}^{\ast }-x_{f}^{\ast }$ remains nearly constant with time for S1, S4 and S5, indicating that the left edge of the wave moves at the same speed as the current tip. On the other hand, for cases S2 and S3 the left edge of the wave travels faster than the current tip, so that a quasisteady flow configuration is not achieved, as will be discussed in the context of figures 6 and 7. Here the relative velocity of the left edge of the wave with respect to the front approaches a constant value $u_{l}^{\ast }$ for long times. In the laboratory frame, the left edge of the wave thus travels with velocity

(2.10) $$\begin{eqnarray}U_{l,lab}^{\ast }=U_{g,lab}^{\ast }+u_{l}^{\ast }.\end{eqnarray}$$

If the flow gives rise to a bore, as will be discussed in the following, the bore speed $u_{b}^{\ast }$ in the moving reference frame can be determined analogously, by tracking the temporal evolution of $x_{b}^{\ast }-x_{f}^{\ast }$ . The bore velocity in the laboratory frame is then obtained as

(2.11) $$\begin{eqnarray}U_{b,lab}^{\ast }=U_{g,lab}^{\ast }+u_{b}^{\ast }.\end{eqnarray}$$

The interface height $h_{l}^{\ast }$ at the left edge of the wave is evaluated as

(2.12) $$\begin{eqnarray}h_{l}^{\ast }=d_{1}^{\ast }(x_{l}^{\ast },t_{s}^{\ast }),\end{eqnarray}$$

where we take $t_{s}^{\ast }=35$ , which falls into the quasisteady range.

The distance of the expansion wave’s right edge from the gravity current tip $x_{r}^{\ast }-x_{f}^{\ast }$ is shown in figure 5(c) for all simulations. In all cases, this distance acquires a quasisteady value for sufficiently large times, which indicates that the right edge of the wave is anchored to the gravity current front.

In order to avoid any potential contamination from the outflow boundary, we evaluate the downstream values of the layer thicknesses some distance upstream of this boundary as

(2.13) $$\begin{eqnarray}\displaystyle & \displaystyle h^{\ast }=d^{\ast }(0.9L_{x},t_{s}^{\ast }), & \displaystyle\end{eqnarray}$$

(2.14) $$\begin{eqnarray}\displaystyle & \displaystyle h_{1o}^{\ast }=d_{1}^{\ast }(0.9L_{x},t_{s}^{\ast }), & \displaystyle\end{eqnarray}$$

(2.15) $$\begin{eqnarray}\displaystyle & \displaystyle h_{2o}^{\ast }=d_{2}^{\ast }(0.9L_{x},t_{s}^{\ast }), & \displaystyle\end{eqnarray}$$

at $t_{s}^{\ast }=35$ . The downstream velocities of the ambient fluid layers are obtained as the average horizontal velocity across the corresponding layer

(2.16) $$\begin{eqnarray}\displaystyle & \displaystyle U_{1o}^{\ast }=\frac{\displaystyle \int _{h^{\ast }}^{h^{\ast }+h_{1o}^{\ast }}u^{\ast }(0.9L_{x},y^{\ast },t_{s}^{\ast })\,\text{d}y^{\ast }}{h_{1o}^{\ast }}, & \displaystyle\end{eqnarray}$$

(2.17) $$\begin{eqnarray}\displaystyle & \displaystyle U_{2o}^{\ast }=\frac{\displaystyle \int _{h^{\ast }+h_{1o}^{\ast }}^{1}u^{\ast }(0.9L_{x},y^{\ast },t_{s}^{\ast })\,\text{d}y^{\ast }}{h_{2o}^{\ast }}, & \displaystyle\end{eqnarray}$$

for $t_{s}^{\ast }=35$ . The DNS results from this section will be compared to predictions by the vorticity model in § 5.

We now proceed to discuss simulations S2, which gives rise to an expansion fan, and S3, which results in the formation of an expansion fan and a bore. These observations are consistent with the findings of Tan et al. (Reference Tan, Nobes, Fleck and Flynn2010) and White & Helfrich (Reference White and Helfrich2012). As shown for simulation S2 in figure 6(c), in the absence of a bore the interfacial slope remains approximately constant over much of the wave region. If a bore is present, on the other hand, the slope of this interface undergoes a noticeable change at some intermediate location, as indicated by the white square in figure 7. For $h_{1i}^{\ast }$ close to one half, the bore becomes weak and increasingly difficult to observe. This is in line with the observations of Rottman & Simpson (Reference Rottman and Simpson1983), who noticed that the bore emerging ahead of the expansion wave in their partial lock-release experiments becomes difficult to identify when the ratio of the lock height to the tank depth is below 0.7. Similarly, the simulations of Khodkar et al. (Reference Khodkar, Nasr-Azadani and Meiburg2017) exhibit a more distinct bore when this ratio is 0.8 or greater.

Figure 6.

Temporal evolution of the density field for simulation S2, at $t^{\ast }=5$ (a), 15 (b) and 35 (c). The simulation parameters are $h_{1i}^{\ast }=0.5$ , $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}^{\ast }=1.4$ , $Re=5000$ and $Pe=20\,000$ . The initial conditions are prescribed according to the simplified model for configuration C2 derived in § 3. The simulation gives rise to an upstream-propagating disturbance in the form of an expansion wave. The white horizontal solid and dashed lines represent the quasisteady gravity current height and the interface location between the ambient layers at the outflow, respectively. The white solid line tangent to the ambient interface shows that the interface height varies approximately linearly along the expansion wave. The white circle, square and cross indicate $x_{f}^{\ast }$ , $x_{l}^{\ast }$ and $x_{r}^{\ast }$ , respectively.

Figure 7.

Snapshot of the density field for simulation S3 at $t^{\ast }=35$ . The simulation parameters are $h_{1i}^{\ast }=0.35$ , $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}^{\ast }=1.4$ , $Re=5000$ and $Pe=20\,000$ . The initial conditions are prescribed according to the simplified model for configuration C3 derived in § 3. The simulation gives rise to an upstream-propagating disturbance in the form of a bore, followed by an expansion wave. The white horizontal solid and dashed lines represent the quasisteady gravity current height and the interface location between the ambient layers at the outflow, respectively. The white solid line tangent to the ambient interface is intended to show that the bore region has a larger slope than the expansion fan region. The white circle, triangle, square and cross indicate $x_{f}^{\ast }$ , $x_{b}^{\ast }$ , $x_{l}^{\ast }$ and $x_{r}^{\ast }$ , respectively.

Figure 8.

Simulation S2: DNS results for the spatial variation of various flow properties inside the wave region. (a) Interface height $d^{\ast }+d_{1}^{\ast }$ , (b) velocity $u_{1}^{\ast }$ of the lower ambient layer, (c) velocity $u_{2}^{\ast }$ of the upper ambient layer and (d) the interfacial vortex sheet strength $\unicode[STIX]{x1D6FE}^{\ast }=u_{2}^{\ast }-u_{1}^{\ast }$ . The results have been obtained for $t^{\ast }=25,30$ and 35, indicated by the dash-dotted, dashed and solid lines, respectively. The vertical dotted lines denote the location of the gravity current tip, with the leftmost line corresponding to $t^{\ast }=25$ and the rightmost one corresponding to $t^{\ast }=35$ . Over much of the expansion wave region, the respective flow properties are seen to vary approximately linearly.

We now analyse the variation of the ambient layer properties within the wave region, in preparation for developing simplified models in § 3. Panels (a–d) of figure 8 display the ambient interface height $d^{\ast }+d_{1}^{\ast }$ , the horizontal flow velocities within the lower and upper ambient layers $u_{1}^{\ast }$ and $u_{2}^{\ast }$ , and the vortex sheet strength $\unicode[STIX]{x1D6FE}^{\ast }=u_{2}^{\ast }-u_{1}^{\ast }$ , within the expansion wave region as functions of the scaled horizontal location $x_{s}^{\ast }=(x^{\ast }-x_{l}^{\ast })/(x_{r}^{\ast }-x_{l}^{\ast })$ , for simulation S2. Here, $u_{1}^{\ast }$ and $u_{2}^{\ast }$ are defined as

(2.18) $$\begin{eqnarray}\displaystyle & \displaystyle u_{1}^{\ast }(x^{\ast },t^{\ast })=\frac{\displaystyle \int _{d^{\ast }}^{d^{\ast }+d_{1}^{\ast }}u^{\ast }(x^{\ast },y^{\ast },t^{\ast })\,\text{d}y^{\ast }}{d_{1}^{\ast }}, & \displaystyle\end{eqnarray}$$

(2.19) $$\begin{eqnarray}\displaystyle & \displaystyle u_{2}^{\ast }(x^{\ast },t^{\ast })=\frac{\displaystyle \int _{d^{\ast }+d_{1}^{\ast }}^{1}u^{\ast }(x^{\ast },y^{\ast },t^{\ast })\,\text{d}y^{\ast }}{d_{2}^{\ast }}. & \displaystyle\end{eqnarray}$$

The results are shown for the three times $t^{\ast }=25$ , 30 and 35, when the Kelvin–Helmholtz instabilities have left the domain and the flow variables far up- and downstream have become time-independent. While the ambient interface height, as well as $u_{1}^{\ast }$ , $u_{2}^{\ast }$ and $\unicode[STIX]{x1D6FE}^{\ast }$ , vary approximately linearly along much of the wave region, they experience more abrupt changes in the vicinity of the gravity current front (vertical dotted lines in figure 8), due to the strong upward acceleration of the ambient fluid caused by the current nose.

We note that in terms of their dimensionless parameters, all three of the above simulations fall into the partially blocked regime identified by Rottman & Simpson (Reference Rottman and Simpson1989) for two-layer stratified flows over solid obstacles: cf. figure 3. However, as mentioned earlier, it is important to keep in mind the subtle differences between gravity current flows and those over solid obstacles; cf. the discussion in Flynn et al. (Reference Flynn, Ungarish and Tan2012). Consequently, the classification introduced by Rottman & Simpson (Reference Rottman and Simpson1989) does not transfer directly to the gravity current problem, which explains why only S2 and S3, but not S1, give rise to upstream-propagating waves.

2.2 Supercritical flow

In addition to the partially blocked flows with $U_{g,lab}^{\ast }>U_{l,lab}^{\ast }$ (or $u_{l}^{\ast }=0$ ), there is another flow regime that does not give rise to upstream waves and can be modelled by C1. This supercritical region was introduced by Baines & Davies (Reference Baines and Davies1980) and is sketched in figure 3. Figure 9(a), which corresponds to S4 in table 1, reflects this flow regime.

Figure 9.

Snapshots of the quasisteady density fields for simulations S4 (a) and S5 (b), at $t^{\ast }=20$ . The simulation parameters are $h_{1i}^{\ast }=0.05$ and $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}^{\ast }=1.89$ , and slightly different values of $Re$ and $Pe$ : cf. table 1. S4 and S5, while both corresponding to configuration C1, are classified as supercritical and partially blocked flows, respectively, by Rottman & Simpson (Reference Rottman and Simpson1989). The white horizontal solid and dashed lines represent the quasisteady gravity current height and the ambient interface location at the outflow. The white circles, squares and crosses indicate $x_{f}^{\ast }$ , $x_{l}^{\ast }$ and $x_{r}^{\ast }$ , respectively.

Consistent with figure 3, in order to achieve supercritical behaviour we take the ratio $D_{o}=h^{\ast }/h_{1i}^{\ast }\approx 1$ , while the Froude number $F_{o}=U_{g}^{\ast }/\sqrt{h_{1i}^{\ast }\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}^{\ast }}$ needs to be sufficiently large. To increase $F_{o}$ we select a small value for $h_{1i}^{\ast }$ and focus on cases that yield large gravity current front velocities. According to the analytical model to be discussed in § 3, and corresponding results of § 4, for very small values of $h_{1i}^{\ast }$ , $U_{g}^{\ast }$ becomes largest when $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}^{\ast }\approx 2$ . Additionally, because $D_{o}\approx 1$ , the gravity current height needs to be small as well. For these reasons, in case S4 we choose $h_{1i}^{\ast }=h^{\ast }=0.05$ . Simulating flows of very thin gravity currents whose density is close to that of the lower ambient fluid layer, requires large values of $Re$ and $Pe$ .

For moderate values of $h_{1i}^{\ast }\gtrsim 0.15$ the simulations typically yielded only a single possible gravity current thickness for each $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}^{\ast }$ . On the other hand, for very small values of $h_{1i}^{\ast }$ we observed that gravity currents with two different heights might form for the same ambient density difference. Simulations S4 and S5, shown in figures 9(a) and 9(b), represent an example in this regard. These simulations have the same values for $h_{1i}^{\ast }$ and $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}^{\ast }$ , and slightly different values of $Re$ and $Pe$ , resulting in very different gravity current heights. According to Rottman & Simpson (Reference Rottman and Simpson1989), S5 would fall into the completely blocked flow regime for a two-layer flow over a solid obstacle. On the other hand, in our DNS simulations for gravity currents propagating into two-layer stratified ambients we never observed a completely blocked flow, i.e. a flow for which the ambient interface does not pass over the gravity current, even for very small lower ambient inflow thicknesses. Again, this reflects the differences between the solid obstacle and gravity current flows alluded to earlier. We note that, to the best of our knowledge, supercritical gravity currents with such small values for $h_{1i}^{\ast }$ and $h^{\ast }$ have not been observed experimentally or numerically in previous investigations.

3.1 Flows without upstream-propagating disturbance

The most basic flow configuration C1, as sketched in figure 1(a), is characterized by a smooth, steady ambient interface without upstream wave propagation. For a given gravity current height $h^{\ast }$ , and prescribed ambient densities and upstream layer thicknesses, we have five unknowns: the velocity $U_{g}^{\ast }$ of the gravity current, as well as the ambient velocities ( $U_{1o}^{\ast }$ , $U_{2o}^{\ast }$ ) and layer thicknesses ( $h_{1o}^{\ast }$ , $h_{2o}^{\ast }$ ) at the outflow boundary. In the reference frame moving with the gravity current front, the dimensionless mass conservation equations for the ambient layers read

(3.1) $$\begin{eqnarray}\displaystyle & \displaystyle U_{g}^{\ast }h_{1i}^{\ast }=U_{1o}^{\ast }h_{1o}^{\ast }, & \displaystyle\end{eqnarray}$$

(3.2) $$\begin{eqnarray}\displaystyle & \displaystyle U_{g}^{\ast }h_{2i}^{\ast }=U_{2o}^{\ast }h_{2o}^{\ast }. & \displaystyle\end{eqnarray}$$

Vorticity conservation along the top of the gravity current $O-2o$ and the interface between the ambient fluid layers $2i-3o$ results in

(3.3) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{1}{2}{U_{1o}^{\ast }}^{2}=\frac{g_{1}^{\prime }}{g^{\prime }}h^{\ast }, & \displaystyle\end{eqnarray}$$

(3.4) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{U_{1o}^{\ast }+U_{2o}^{\ast }}{2}(U_{2o}^{\ast }-U_{1o}^{\ast })=\frac{g_{2}^{\prime }}{g^{\prime }}(h_{1o}^{\ast }+h^{\ast }-h_{1i}^{\ast }), & \displaystyle\end{eqnarray}$$

where

(3.5) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{g_{1}^{\prime }}{g^{\prime }}=\frac{\displaystyle \unicode[STIX]{x1D70C}_{1}-\left(\unicode[STIX]{x1D70C}_{2}+\frac{\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}}{2}\right)}{\unicode[STIX]{x1D70C}_{1}-\unicode[STIX]{x1D70C}_{2}}, & \displaystyle\end{eqnarray}$$

(3.6) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{g_{2}^{\prime }}{g^{\prime }}=\frac{\displaystyle \left(\unicode[STIX]{x1D70C}_{2}+\frac{\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}}{2}\right)-\left(\unicode[STIX]{x1D70C}_{2}-\frac{\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}}{2}\right)}{\unicode[STIX]{x1D70C}_{1}-\unicode[STIX]{x1D70C}_{2}}=\frac{\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x1D70C}_{1}-\unicode[STIX]{x1D70C}_{2}}. & \displaystyle\end{eqnarray}$$

A fifth equation is obtained from the geometrical constraint at the outflow boundary

(3.7) $$\begin{eqnarray}h^{\ast }+h_{1o}^{\ast }+h_{2o}^{\ast }=1.\end{eqnarray}$$

The above equations (3.1)–(3.4) and (3.7) represent a closed algebraic system for the five unknowns $U_{g}^{\ast }$ , $U_{1o}^{\ast }$ , $U_{2o}^{\ast }$ , $h_{1o}^{\ast }$ and $h_{2o}^{\ast }$ , so that empirical closure assumptions are not needed.

As an aside, we insert here a brief comment on the relationship between the head losses in the upper and lower layers, for the quasisteady flow of C1 in the moving reference frame. Within each layer, the head loss is of course constant along all streamlines, since the flow is irrotational. We write Bernoulli’s equation with a head loss in dimensionless form, for fluid elements immediately below and above the streamline $2i-3o$ :

(3.8) $$\begin{eqnarray}\displaystyle & \displaystyle p_{2i}^{\ast }+\frac{1}{2}{U_{g}^{\ast }}^{2}+\frac{\displaystyle \left(\unicode[STIX]{x1D70C}_{2}+\frac{\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}}{2}\right)g}{\unicode[STIX]{x1D70C}_{2}g^{\prime }}h_{1i}^{\ast }=p_{3o}^{\ast }+\frac{1}{2}{U_{1o}^{\ast }}^{2}+\frac{\displaystyle \left(\unicode[STIX]{x1D70C}_{2}+\frac{\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}}{2}\right)g}{\unicode[STIX]{x1D70C}_{2}g^{\prime }}(h^{\ast }+h_{1o}^{\ast })+\unicode[STIX]{x1D6E5}_{1}^{\ast },\qquad & \displaystyle\end{eqnarray}$$

(3.9) $$\begin{eqnarray}\displaystyle & \displaystyle p_{2i}^{\ast }+\frac{1}{2}{U_{g}^{\ast }}^{2}+\frac{\displaystyle \left(\unicode[STIX]{x1D70C}_{2}-\frac{\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}}{2}\right)g}{\unicode[STIX]{x1D70C}_{2}g^{\prime }}h_{1i}^{\ast }=p_{3o}^{\ast }+\frac{1}{2}{U_{2o}^{\ast }}^{2}+\frac{\displaystyle \left(\unicode[STIX]{x1D70C}_{2}-\frac{\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}}{2}\right)g}{\unicode[STIX]{x1D70C}_{2}g^{\prime }}(h^{\ast }+h_{1o}^{\ast })+\unicode[STIX]{x1D6E5}_{2}^{\ast }.\qquad & \displaystyle\end{eqnarray}$$

Here $\unicode[STIX]{x1D6E5}_{1}^{\ast }$ and $\unicode[STIX]{x1D6E5}_{2}^{\ast }$ represent the dimensionless head losses within the lower and upper ambient layers, respectively, scaled by $\unicode[STIX]{x1D6E5}_{ref}=\unicode[STIX]{x1D70C}_{2}g^{\prime }H$ . Subtracting (3.9) from (3.8) yields

(3.10) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6E5}_{1}^{\ast }-\unicode[STIX]{x1D6E5}_{2}^{\ast }=\frac{1}{2}({U_{2o}^{\ast }}^{2}-{U_{1o}^{\ast }}^{2})-\frac{g_{2}^{\prime }}{g^{\prime }}(h^{\ast }+h_{1o}^{\ast }-h_{1i}^{\ast }). & & \displaystyle\end{eqnarray}$$

The right-hand side of this equation vanishes according to (3.4), and we conclude that vorticity conservation immediately implies that both ambient layers have the same head loss. We will return to this point when discussing the assumptions underlying previous models of earlier investigators in § 6.

3.2 Propagation of a rarefaction wave towards the inflow boundary

Figure 1(b) displays configuration C2, a flow with a rarefaction wave travelling towards the inflow boundary, which causes the steady-state assumption to be invalid, even in the reference frame moving with the gravity current tip. The DNS results of § 2 show that the right edge $R$ of the wave remains stationary in the reference frame of the gravity current. The speed $U_{l}$ of the left edge $L$ of the wave is merely a function of the inflow conditions

(3.11) $$\begin{eqnarray}U_{l}^{\ast }=\sqrt{\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}^{\ast }h_{1i}^{\ast }h_{2i}^{\ast }}\end{eqnarray}$$

(cf. Khodkar et al. Reference Khodkar, Nasr-Azadani and Meiburg2017). We then need five additional equations to determine the remaining unknowns $U_{g}^{\ast }$ , $U_{1o}^{\ast }$ , $U_{2o}^{\ast }$ , $h_{1o}^{\ast }$ and $h_{2o}^{\ast }$ . These can again be obtained from the mass balances of the two ambient layers, the vorticity balances along the two interfaces and the geometrical constraint provided by (3.7), although the unsteadiness renders the equations slightly more complex. While the vorticity balance along the steady upper boundary of the gravity current is still given by (3.3), the unsteady evolution of the rarefaction wave modifies the other three equations.

In order to quantify the temporal changes of the ambient layer volumes, and of the circulation along the boundary separating these layers, we assume that the ambient interface height as well as the vortex sheet strength vary linearly within the wave region from $L$ to $R$ . The DNS results of § 2 above showed that this approximation holds to a good degree over much of the wave region, although it becomes questionable in the region immediately above the gravity current front: cf. figure 8(a). Hence it will be interesting to compare the predictions obtained on the basis of this approximation with the actual DNS results further below. Consequently, in the reference frame of the gravity current the mass balance for the deforming upper layer within control volume $L-3o-4o-L^{\prime }$ yields

(3.12) $$\begin{eqnarray}(U_{l}^{\ast }-U_{g}^{\ast })\frac{h_{2i}^{\ast }+h_{2o}^{\ast }}{2}=U_{l}^{\ast }h_{2i}^{\ast }-U_{2o}^{\ast }h_{2o}^{\ast }.\end{eqnarray}$$

Mass conservation for the entire channel $1i-1o-4o-3i$ results in

(3.13) $$\begin{eqnarray}\displaystyle U_{g}^{\ast }=U_{1o}^{\ast }h_{1o}^{\ast }+U_{2o}^{\ast }h_{2o}^{\ast }. & & \displaystyle\end{eqnarray}$$

To obtain an equation for the vorticity along the ambient interface $2i-3o$ , we start from the integral form of the inviscid vorticity conservation equation

(3.14) $$\begin{eqnarray}\displaystyle \frac{\text{d}}{\text{d}t^{\ast }}\int _{CV(t^{\ast })}\unicode[STIX]{x1D714}^{\ast }\,\text{d}A^{\ast }+\oint _{CS(t^{\ast })}\unicode[STIX]{x1D714}\boldsymbol{V}_{r}^{\ast }\boldsymbol{\cdot }\boldsymbol{n}\,\text{d}S^{\ast }=-\int _{CV(t^{\ast })}\frac{g_{2}^{\prime }}{g^{\prime }}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}^{\ast }}{\unicode[STIX]{x2202}x^{\ast }}\,\text{d}A^{\ast }, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D714}^{\ast }=\unicode[STIX]{x1D714}/\sqrt{g^{\prime }/H}$ , and $CS$ denotes the boundary of the control volume that contains the ambient interface. $\boldsymbol{V}_{r}^{\ast }$ represents the relative flow velocity at this boundary. The temporal rate of change of the circulation inside the control volume is due to the elongation of the wave region resulting from the rarefaction wave propagation. The rate at which the wave region lengthens is given by $U_{l}^{\ast }-U_{g}^{\ast }$ . Based on the DNS results of § 2 (cf. figure 8 d), we assume that the vortex sheet strength varies linearly from $L$ to $R$ , so that its average value in the wave region is $(U_{2o}^{\ast }-U_{1o}^{\ast })/2$ . Hence we obtain

(3.15) $$\begin{eqnarray}\displaystyle \frac{\text{d}}{\text{d}t^{\ast }}\int _{CV(t^{\ast })}\unicode[STIX]{x1D714}^{\ast }\,\text{d}A^{\ast }=-\frac{1}{2}(U_{l}^{\ast }-U_{g}^{\ast })(U_{2o}^{\ast }-U_{1o}^{\ast }). & & \displaystyle\end{eqnarray}$$

There is no influx of vorticity into the control volume, and the shear at the outflow boundary yields

(3.16) $$\begin{eqnarray}\displaystyle \oint _{CS(t)}\unicode[STIX]{x1D714}^{\ast }\boldsymbol{V}_{r}^{\ast }\boldsymbol{\cdot }\boldsymbol{n}\,\text{d}S^{\ast }=-\frac{1}{2}(U_{2o}^{\ast }-U_{1o}^{\ast })(U_{2o}^{\ast }+U_{1o}^{\ast }). & & \displaystyle\end{eqnarray}$$

The baroclinic vorticity production term can be evaluated as

(3.17) $$\begin{eqnarray}\displaystyle -\!\int _{CV(t^{\ast })}\frac{g_{2}^{\prime }}{g^{\prime }}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}^{\ast }}{\unicode[STIX]{x2202}x^{\ast }}\,\text{d}A^{\ast }=-\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}^{\ast }(h_{1o}^{\ast }+h^{\ast }-h_{1i}^{\ast }). & & \displaystyle\end{eqnarray}$$

By substituting (3.15)–(3.17) into (3.14), we obtain the vorticity equation in the form

(3.18) $$\begin{eqnarray}\displaystyle \frac{1}{2}(U_{l}^{\ast }-U_{g}^{\ast })(U_{2o}^{\ast }-U_{1o}^{\ast })+\frac{{U_{2o}^{\ast }}^{2}-{U_{1o}^{\ast }}^{2}}{2}=\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}^{\ast }(h^{\ast }+h_{1o}^{\ast }-h_{1i}^{\ast }). & & \displaystyle\end{eqnarray}$$

The system of equations (3.3), (3.12), (3.13) and (3.18) is then fully closed by the geometrical constraint (3.7).

3.3 Rarefaction wave following an upstream bore

When the rarefaction wave gives rise to a leading internal bore, we model the flow based on configuration C3, as sketched in figure 1(c). As shown by Khodkar et al. (Reference Khodkar, Nasr-Azadani and Meiburg2017), the mass and vorticity conservation equations along the rarefaction wave and in the laboratory reference frame can be written as

(3.19) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}u_{1}^{\ast }}{\text{d}d_{1}^{\ast }}=-\frac{u^{\ast }+u_{1}^{\ast }-U_{g}^{\ast }}{d_{1}^{\ast }}, & \displaystyle\end{eqnarray}$$

(3.20) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}u_{2}^{\ast }}{\text{d}d_{1}^{\ast }}=\frac{u^{\ast }+u_{2}^{\ast }-U_{g}^{\ast }}{d_{2}^{\ast }}, & \displaystyle\end{eqnarray}$$

(3.21) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{(u^{\ast }+u_{1}^{\ast }-U_{g}^{\ast })^{2}}{d_{1}^{\ast }}+\frac{(u^{\ast }+u_{2}^{\ast }-U_{g}^{\ast })^{2}}{d_{2}^{\ast }}=\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}^{\ast }. & \displaystyle\end{eqnarray}$$

We recall from § 2 that $u_{1}^{\ast }$ and $u_{2}^{\ast }$ represent the dimensionless horizontal velocities of the lower and upper layers, while $d_{1}^{\ast }$ and $d_{2}^{\ast }$ denote the local depths of the layers. $u^{\ast }$ indicates the interfacial velocity along the wave in the laboratory frame, which is positive towards the inflow boundary. Consequently, $u^{\ast }(x_{r}^{\ast })=U_{g}^{\ast }$ and $u^{\ast }(x_{l}^{\ast })=U_{l}^{\ast }$ . The above system of equations, which can be obtained by applying the mass and vorticity conservation principles to an infinitesimal control volume along the expansion fan, yields the flow field within the fan as a function of the interface height $d_{1}^{\ast }(x^{\ast },t^{\ast })$ , without invoking any assumptions regarding the interfacial shape: cf. Khodkar et al. (Reference Khodkar, Nasr-Azadani and Meiburg2017).

If $|u^{\ast }(x^{\ast })|$ decreases monotonically from left to right, only the expansion wave forms. On the other hand, if $|u^{\ast }(x^{\ast })|$ has a maximum somewhere between $L$ and $R$ , the interface steepens locally and a bore emerges as shown in figure 1(c). To investigate the conditions under which a maximum appears, we take the derivative of (3.19)

(3.22) $$\begin{eqnarray}\displaystyle & & \displaystyle \left[\frac{2(u^{\ast }+u_{2}^{\ast }-U_{g}^{\ast })}{1-d_{1}^{\ast }}+\frac{2(u^{\ast }+u_{1}^{\ast }-U_{g}^{\ast })}{d_{1}^{\ast }}\right]\frac{\text{d}u^{\ast }}{\text{d}d_{1}^{\ast }}\nonumber\\ \displaystyle & & \displaystyle \quad -\,\frac{3(u^{\ast }+u_{1}^{\ast }-U_{g}^{\ast })^{2}}{{d_{1}^{\ast }}^{2}}+\frac{3(u^{\ast }+u_{2}^{\ast }-U_{g}^{\ast })^{2}}{(1-d_{1}^{\ast })^{2}}=0,\end{eqnarray}$$

where $\text{d}u_{1}^{\ast }/\text{d}d_{1}^{\ast }$ and $\text{d}u_{2}^{\ast }/\text{d}d_{1}^{\ast }$ have been substituted from (3.19) and (3.20), respectively. By letting $\text{d}u^{\ast }/\text{d}d_{1}^{\ast }$ be zero, we arrive at

(3.23) $$\begin{eqnarray}\frac{u^{\ast }+u_{1}^{\ast }-U_{g}^{\ast }}{d_{1}^{\ast }}=\frac{u^{\ast }+u_{2}^{\ast }-U_{g}^{\ast }}{1-d_{1}^{\ast }}.\end{eqnarray}$$

We expect that the local extremum first appears for the smallest interface height, i.e. at the left edge of the expansion wave, where $d_{1}^{\ast }=h_{1i}^{\ast }$ and $u_{1}^{\ast }-U_{g}^{\ast }=u_{2}^{\ast }-U_{g}^{\ast }=0$ . By substituting these values into (3.23), we find that a bore emerges only if $h_{1i}^{\ast }<0.5$ . Khodkar et al. (Reference Khodkar, Nasr-Azadani and Meiburg2017) show that in the presence of a bore, the flow chooses $h_{l}^{\ast }=0.5$ , which results in a non-dissipative flow across the bore. The equations for mass and vorticity conservation around the internal bore, and in the reference frame moving with the bore, yield

(3.24) $$\begin{eqnarray}\displaystyle & \displaystyle U_{b}^{\ast }h_{1i}^{\ast }=(U_{b}^{\ast }-U_{1b}^{\ast })h_{l}^{\ast }, & \displaystyle\end{eqnarray}$$

(3.25) $$\begin{eqnarray}\displaystyle & \displaystyle U_{b}^{\ast }h_{2i}^{\ast }=(U_{b}^{\ast }+U_{2b}^{\ast })(1-h_{l}^{\ast }), & \displaystyle\end{eqnarray}$$

(3.26) $$\begin{eqnarray}\displaystyle & \displaystyle {\textstyle \frac{1}{2}}(U_{2b}^{\ast }+U_{1b}^{\ast })(U_{2b}^{\ast }-U_{1b}^{\ast }+2U_{b}^{\ast })=\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}^{\ast }(h_{l}^{\ast }-h_{1i}^{\ast }), & \displaystyle\end{eqnarray}$$

where $U_{b}^{\ast }$ is the bore propagation speed, and $U_{1b}^{\ast }$ and $U_{2b}^{\ast }$ represent the velocities of the lower and upper layers downstream of the bore in the laboratory reference frame. Equations (3.24)–(3.26) can be solved to obtain

(3.27) $$\begin{eqnarray}\displaystyle & \displaystyle U_{b}^{\ast }=0.5\sqrt{\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}^{\ast }}, & \displaystyle\end{eqnarray}$$

(3.28) $$\begin{eqnarray}\displaystyle & \displaystyle U_{1b}^{\ast }=U_{2b}^{\ast }=\sqrt{\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}^{\ast }}(0.5-h_{1i}^{\ast }). & \displaystyle\end{eqnarray}$$

$U_{1b}^{\ast }$ and $U_{2b}^{\ast }$ can then be substituted into the vorticity balance equation (3.21) to arrive at a relationship for the propagation velocity of the left edge of the rarefaction wave

(3.29) $$\begin{eqnarray}U_{l}^{\ast }=\sqrt{\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}^{\ast }h_{1i}^{\ast }h_{2i}^{\ast }},\end{eqnarray}$$

which is identical to that found for C2 in (3.11).

After finding the bore speed as well as the velocities of the ambient layers downstream of the bore, the gravity current velocity and the heights and velocities of the ambient layers at the outflow boundary can then be obtained by an analysis corresponding to those of §§ 3.1 and 3.2. We thus obtain

(3.30) $$\begin{eqnarray}\displaystyle & \displaystyle (U_{l}^{\ast }-U_{g}^{\ast })\frac{1-h_{l}^{\ast }+h_{2o}^{\ast }}{2}=(U_{2b}^{\ast }+U_{l}^{\ast })(1-h_{l}^{\ast })-U_{2o}^{\ast }h_{2o}^{\ast }, & \displaystyle\end{eqnarray}$$

(3.31) $$\begin{eqnarray}\displaystyle & & \displaystyle \frac{1}{2}(U_{l}^{\ast }-U_{g}^{\ast })(U_{2o}^{\ast }-U_{1o}^{\ast }+U_{1b}^{\ast }+U_{2b}^{\ast })+\frac{{U_{2o}^{\ast }}^{2}-{U_{1o}^{\ast }}^{2}}{2}\nonumber\\ \displaystyle & & \displaystyle \quad -\,\frac{1}{2}(U_{1b}^{\ast }+U_{2b}^{\ast })(U_{2b}^{\ast }-U_{1b}^{\ast }+2U_{l}^{\ast })=\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}^{\ast }(h^{\ast }+h_{1o}^{\ast }-h_{l}^{\ast }),\end{eqnarray}$$

where the vortex sheet strength again has been assumed to vary linearly. Equations (3.3), (3.7), (3.13), (3.30) and (3.31) can then be solved to obtain the flow quantities for configuration C3. Alternatively, one may consider applying two-layer shallow-water equations to the rarefaction wave region, or developing a model such as that of Khodkar et al. (Reference Khodkar, Nasr-Azadani and Meiburg2017), which concurrently applies mass and vorticity conservation principles to infinitesimal control volumes within the expansion fan.

3.4 Physically possible solutions

Since no energy-based closure assumptions were employed while deriving the above analytical models, the energetics of the flow can be assessed a posteriori, i.e. after determining the flow field. Towards this end, we will now analyse the rate of energy loss for the entire flow in the reference frame of the gravity current. The most general form of the energy equation is

(3.32) $$\begin{eqnarray}\displaystyle {\dot{E}}_{i}^{\ast }-{\dot{E}}_{o}^{\ast }-{\dot{E}}_{d}^{\ast }-{\dot{E}}_{s}^{\ast }=0. & & \displaystyle\end{eqnarray}$$

Here ${\dot{E}}_{d}^{\ast }$ and ${\dot{E}}_{s}^{\ast }$ represent the dimensionless rates of dissipation and storage of energy inside the control volume, while ${\dot{E}}_{i}^{\ast }$ and ${\dot{E}}_{o}^{\ast }$ indicate the rates at which energy is convected into and out of the control volume. These rates have been rendered dimensionless by $\unicode[STIX]{x1D70C}_{2}{g^{\prime }}^{3/2}H^{5/2}$ .

We begin by discussing the steady configuration C1, for which the storage term vanishes. ${\dot{E}}_{i}^{\ast }$ and ${\dot{E}}_{o}^{\ast }$ can be calculated as

(3.33) $$\begin{eqnarray}\displaystyle & \displaystyle {\dot{E}}_{i}^{\ast }=\int _{0}^{1}\left(p^{\ast }(y^{\ast })+\frac{1}{2}{U_{i}^{\ast }}^{2}+\frac{\unicode[STIX]{x1D70C}g}{\unicode[STIX]{x1D70C}_{2}g^{\prime }}y^{\ast }\right)U_{i}^{\ast }\,\text{d}y^{\ast }, & \displaystyle\end{eqnarray}$$

(3.34) $$\begin{eqnarray}\displaystyle & \displaystyle {\dot{E}}_{o}^{\ast }=\int _{0}^{1}\left(p^{\ast }(y^{\ast })+\frac{1}{2}{U_{o}^{\ast }}^{2}+\frac{\unicode[STIX]{x1D70C}g}{\unicode[STIX]{x1D70C}_{2}g^{\prime }}y^{\ast }\right)U_{o}^{\ast }\,\text{d}y^{\ast }. & \displaystyle\end{eqnarray}$$

Since the flow direction at the in- and outflow boundaries is horizontal, we can assume the local pressure distribution to be hydrostatic

(3.35) $$\begin{eqnarray}\displaystyle & \displaystyle p_{i}^{\ast }(y^{\ast })=\left\{\begin{array}{@{}ll@{}}\displaystyle p_{3i}^{\ast }+\frac{\displaystyle \left(\unicode[STIX]{x1D70C}_{2}-\frac{\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}}{2}\right)g}{\unicode[STIX]{x1D70C}_{2}g^{\prime }}(1-y^{\ast })\quad & y^{\ast }\geqslant h_{1i}^{\ast },\\ \displaystyle p_{3i}^{\ast }+\frac{\displaystyle \left(\unicode[STIX]{x1D70C}_{2}-\frac{\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}}{2}\right)g}{\unicode[STIX]{x1D70C}_{2}g^{\prime }}h_{2i}^{\ast }+\frac{\displaystyle \left(\unicode[STIX]{x1D70C}_{2}+\frac{\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}}{2}\right)g}{\unicode[STIX]{x1D70C}_{2}g^{\prime }}(h_{1i}^{\ast }-y^{\ast })\quad & y^{\ast }<h_{1i}^{\ast }.\end{array}\right. & \displaystyle\end{eqnarray}$$

(3.36) $$\begin{eqnarray}\displaystyle & \displaystyle p_{o}^{\ast }(y^{\ast })=\left\{\begin{array}{@{}l@{}}\displaystyle p_{4o}^{\ast }+\frac{\displaystyle \left(\unicode[STIX]{x1D70C}_{2}-\frac{\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}}{2}\right)g}{\unicode[STIX]{x1D70C}_{2}g^{\prime }}(1-y^{\ast })\quad \\ \quad \qquad y\geqslant h^{\ast }+h_{1o}^{\ast },\quad \\ \displaystyle p_{4o}^{\ast }+\frac{\displaystyle \left(\unicode[STIX]{x1D70C}_{2}-\frac{\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}}{2}\right)g}{\unicode[STIX]{x1D70C}_{2}g^{\prime }}h_{2o}^{\ast }+\frac{\displaystyle \left(\unicode[STIX]{x1D70C}_{2}+\frac{\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}}{2}\right)g}{\unicode[STIX]{x1D70C}_{2}g^{\prime }}(h^{\ast }+h_{1o}^{\ast }-y^{\ast })\quad \\ \quad \qquad h^{\ast }\leqslant y^{\ast }<h^{\ast }+h_{1o}^{\ast },\quad \\ \displaystyle p_{4o}^{\ast }+\frac{\displaystyle \left(\unicode[STIX]{x1D70C}_{2}-\frac{\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}}{2}\right)g}{\unicode[STIX]{x1D70C}_{2}g^{\prime }}h_{2o}^{\ast }+\frac{\displaystyle \left(\unicode[STIX]{x1D70C}_{2}+\frac{\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}}{2}\right)g}{\unicode[STIX]{x1D70C}_{2}g^{\prime }}h_{1o}^{\ast }+\frac{\unicode[STIX]{x1D70C}_{1}g}{\unicode[STIX]{x1D70C}_{2}g^{\prime }}(h^{\ast }-y^{\ast })\quad \\ \quad \qquad y^{\ast }<h^{\ast }.\quad \end{array}\right. & \displaystyle\end{eqnarray}$$

Here, $p_{i}^{\ast }(y^{\ast })$ and $p_{o}^{\ast }(y^{\ast })$ denote the dimensionless pressure at the in- and outflow boundaries, respectively. The pressure drop $p_{3i}^{\ast }-p_{4o}^{\ast }$ along the top wall of the channel can be computed in a straightforward fashion, by applying the integral horizontal momentum balance for the entire channel in the reference frame moving with the gravity current

(3.37) $$\begin{eqnarray}{\dot{M}}^{\ast }=\int _{0}^{1}(p_{i}^{\ast }+{U_{i}^{\ast }}^{2})\,\text{d}y^{\ast }-\int _{0}^{1}(p_{o}^{\ast }+{U_{o}^{\ast }}^{2})\,\text{d}y^{\ast }.\end{eqnarray}$$

Here ${\dot{M}}^{\ast }$ has been rendered dimensionless by $\unicode[STIX]{x1D70C}_{2}g^{\prime }H^{2}$ . We note that for C1, ${\dot{M}}^{\ast }=0$ , due to the steady nature of the flow. By substituting the velocities evaluated in §§ 3.1–3.3, and the pressure functions of (3.35) and (3.36) into (3.37), we obtain the following expression for the pressure drop along the top wall:

(3.38) $$\begin{eqnarray}\displaystyle p_{3i}^{\ast }-p_{4o}^{\ast } & = & \displaystyle \frac{\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}^{\ast }}{2}\left(\frac{{h_{1o}^{\ast }}^{2}-{h_{2o}^{\ast }}^{2}}{2}\right)-\frac{\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}^{\ast }}{2}h_{1o}^{\ast }h_{2o}^{\ast }-\frac{\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}^{\ast }}{2}h^{\ast }h_{2o}^{\ast }\nonumber\\ \displaystyle & & \displaystyle +\,\frac{\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}^{\ast }}{2}h^{\ast }h_{1o}^{\ast }+\frac{\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}^{\ast }}{2}\left(\frac{{h_{2i}^{\ast }}^{2}-{h_{1i}^{\ast }}^{2}}{2}\right)+\frac{\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}^{\ast }}{2}h_{1i}^{\ast }h_{2i}^{\ast }\nonumber\\ \displaystyle & & \displaystyle +\,\frac{{h^{\ast }}^{2}}{2}+{U_{1o}^{\ast }}^{2}h_{1o}^{\ast }+{U_{2o}^{\ast }}^{2}h_{2o}^{\ast }-{U_{g}^{\ast }}^{2}+{\dot{M}}^{\ast }.\end{eqnarray}$$

For C2 and C3 the unsteady rarefaction wave results in ${\dot{M}}={\dot{M}}_{1}^{\ast }+{\dot{M}}_{2}^{\ast }\neq 0$ , where ${\dot{M}}_{1}^{\ast }$ and ${\dot{M}}_{2}^{\ast }$ indicate the rates at which momentum accumulates in the lower and upper layers, respectively. These rates can be obtained as

(3.39) $$\begin{eqnarray}\displaystyle & \displaystyle {\dot{M}}_{1}^{\ast }=\left[\frac{\text{d}}{\text{d}t^{\ast }}\int _{x_{l}^{\ast }}^{x_{r}^{\ast }}u_{1}^{\ast }(x^{\ast })d_{1}^{\ast }(x^{\ast })\,\text{d}x^{\ast }\right]-(U_{l}^{\ast }-U_{g}^{\ast })U_{g}^{\ast }h_{1i}^{\ast }, & \displaystyle\end{eqnarray}$$

(3.40) $$\begin{eqnarray}\displaystyle & \displaystyle {\dot{M}}_{2}^{\ast }=\left[\frac{\text{d}}{\text{d}t^{\ast }}\int _{x_{l}^{\ast }}^{x_{r}^{\ast }}u_{2}^{\ast }(x^{\ast })d_{2}^{\ast }(x^{\ast })\,\text{d}x^{\ast }\right]-(U_{l}^{\ast }-U_{g}^{\ast })U_{g}^{\ast }h_{2i}^{\ast }, & \displaystyle\end{eqnarray}$$

where $x_{l}^{\ast }$ and $x_{r}^{\ast }$ denote the horizontal locations of the left and right edges of the wave, as defined in § 2. In order to be able to evaluate the above integrals analytically, we assume that the horizontal velocities of the individual layers vary linearly as a function of $x^{\ast }$ along the wave, which is supported by the DNS data of the previous section. The validity of this assumption in the context of evaluating the energy budget of the flow, will be further investigated in § 5. The summation of ${\dot{M}}_{1}^{\ast }$ and ${\dot{M}}_{2}^{\ast }$ then gives the net accumulation rate ${\dot{M}}^{\ast }$ of momentum within the entire channel. The pressure drop along the top wall can be evaluated by substituting ${\dot{M}}^{\ast }$ into (3.38).

${\dot{E}}_{s}^{\ast }$ is the rate at which the energy within the channel changes with time, as the rarefaction wave travels to the left. To assess the rate at which the kinetic energy of the flow changes, it is again convenient to analyse each layer separately. The rate of change of the kinetic energy of the lower and upper layers can be calculated as

(3.41) $$\begin{eqnarray}\displaystyle & \displaystyle {\dot{E}}_{k,1}^{\ast }=\left[\frac{\text{d}}{\text{d}t^{\ast }}\int _{x_{l}^{\ast }}^{x_{r}^{\ast }}\frac{1}{2}{u_{1}^{\ast }}^{2}(x^{\ast })d_{1}^{\ast }(x^{\ast })\,\text{d}x^{\ast }\right]-\frac{1}{2}(U_{l}^{\ast }-U_{g}^{\ast }){U_{g}^{\ast }}^{2}h_{1i}^{\ast }, & \displaystyle\end{eqnarray}$$

(3.42) $$\begin{eqnarray}\displaystyle & \displaystyle {\dot{E}}_{k,2}^{\ast }=\left[\frac{\text{d}}{\text{d}t^{\ast }}\int _{x_{l}^{\ast }}^{x_{r}^{\ast }}\frac{1}{2}{u_{2}^{\ast }}^{2}(x^{\ast })d_{2}^{\ast }(x^{\ast })\,\text{d}x^{\ast }\right]-\frac{1}{2}(U_{l}^{\ast }-U_{g}^{\ast }){U_{g}^{\ast }}^{2}h_{2i}^{\ast }. & \displaystyle\end{eqnarray}$$

${\dot{E}}_{k,1}^{\ast }$ and ${\dot{E}}_{k,2}^{\ast }$ can then be added to calculate the net rate of change of kinetic energy ${\dot{E}}_{k}^{\ast }$ within the channel.

The propagation of the expansion wave furthermore lifts up the denser fluid of the lower layer, thereby increasing the potential energy. To quantify the rate of growth of potential energy, we analyse the evolution of the constant-slope interface over a small time interval $\unicode[STIX]{x0394}t$ , due to the propagation of the rarefaction wave, as shown in figure 10. The shaded region indicates the volume of light, upper layer fluid replaced by dense, lower layer fluid during $\unicode[STIX]{x0394}t$ . The rate at which the potential energy increases is then proportional to the rate at which the area of this triangle grows times the height of its centre of gravity, which results in

(3.43) $$\begin{eqnarray}{\dot{E}}_{p}^{\ast }=\frac{1}{2}\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}^{\ast }(U_{l}^{\ast }-U_{g}^{\ast })(h^{\ast }+h_{1o}^{\ast }-h_{1i}^{\ast })\left(\frac{h^{\ast }+h_{1o}^{\ast }+2h_{1i}^{\ast }}{3}\right).\end{eqnarray}$$

The potential energy rate of growth can be added to the kinetic energy rate of growth to obtain ${\dot{E}}_{s}^{\ast }$ . We then substitute ${\dot{E}}_{s}^{\ast }$ into (3.32) to compute the rate of energy dissipation ${\dot{E}}_{d}^{\ast }$ .

Figure 10.

Sketch in support of evaluating the unsteady potential energy term, in the reference frame moving with the gravity current. Over a small time interval $\unicode[STIX]{x0394}t$ , the left edge of the expansion wave travels from $L$ to $L_{n}$ with the constant velocity $U_{l}-U_{g}$ , while the right edge $R$ remains stationary. Within the expansion wave region, the interface is assumed to have a constant slope that changes with time. The shaded region indicates the volume of dense fluid of the lower layer that has replaced the light fluid of the upper layer during $\unicode[STIX]{x0394}t$ .

A corresponding analysis can be conducted for configuration C3 in order to obtain the rate of energy dissipation. Towards this end, we first compute the pressure drop along the top wall from (3.38), where the rate of accumulation of momentum ${\dot{M}}^{\ast }$ is now evaluated according to

(3.44) $$\begin{eqnarray}\displaystyle {\dot{M}}_{1}^{\ast } & = & \displaystyle \left[\frac{\text{d}}{\text{d}t^{\ast }}\int _{x_{l}^{\ast }}^{x_{r}^{\ast }}u_{1}^{\ast }(x^{\ast })d_{1}^{\ast }(x^{\ast })\,\text{d}x^{\ast }\right]+(U_{b}^{\ast }-U_{l}^{\ast })(U_{g}^{\ast }-U_{1b}^{\ast })h_{l}^{\ast }\nonumber\\ \displaystyle & & \displaystyle -\,(U_{b}^{\ast }-U_{g}^{\ast })U_{g}^{\ast }h_{1i}^{\ast },\end{eqnarray}$$

(3.45) $$\begin{eqnarray}\displaystyle {\dot{M}}_{2}^{\ast } & = & \displaystyle \left[\frac{\text{d}}{\text{d}t^{\ast }}\int _{x_{l}^{\ast }}^{x_{r}^{\ast }}u_{2}^{\ast }(x^{\ast })d_{2}^{\ast }(x^{\ast })\,\text{d}x^{\ast }\right]+(U_{b}^{\ast }-U_{l}^{\ast })(U_{g}^{\ast }+U_{2b}^{\ast })(1-h_{l}^{\ast })\nonumber\\ \displaystyle & & \displaystyle -\,(U_{b}^{\ast }-U_{g}^{\ast })U_{g}^{\ast }h_{2i}^{\ast }.\end{eqnarray}$$

The results of the above integrals can then be added up to obtain ${\dot{M}}^{\ast }$ for C3, which then enables us to compute $p_{3i}^{\ast }-p_{4o}^{\ast }$ from (3.38). In a corresponding fashion, we obtain

(3.46) $$\begin{eqnarray}\displaystyle {\dot{E}}_{k,1}^{\ast } & = & \displaystyle \left[\frac{\text{d}}{\text{d}t^{\ast }}\int _{x_{l}^{\ast }}^{x_{r}^{\ast }}\frac{1}{2}{u_{1}^{\ast }}^{2}(x^{\ast })d_{1}^{\ast }(x^{\ast })\,\text{d}x^{\ast }\right]+\frac{1}{2}(U_{b}^{\ast }-U_{l}^{\ast })(U_{g}^{\ast }-U_{1b}^{\ast })^{2}h_{l}^{\ast }\nonumber\\ \displaystyle & & \displaystyle -\,\frac{1}{2}(U_{b}^{\ast }-U_{g}^{\ast }){U_{g}^{\ast }}^{2}h_{1i}^{\ast },\end{eqnarray}$$

(3.47) $$\begin{eqnarray}\displaystyle {\dot{E}}_{k,2}^{\ast } & = & \displaystyle \left[\frac{\text{d}}{\text{d}t^{\ast }}\int _{x_{l}^{\ast }}^{x_{r}^{\ast }}\frac{1}{2}{u_{2}^{\ast }}^{2}(x^{\ast })d_{2}^{\ast }(x^{\ast })\,\text{d}x^{\ast }\right]+\frac{1}{2}(U_{b}^{\ast }-U_{l}^{\ast })(U_{g}^{\ast }+U_{2b}^{\ast })^{2}(1-h_{l}^{\ast })\nonumber\\ \displaystyle & & \displaystyle -\,\frac{1}{2}(U_{b}^{\ast }-U_{g}^{\ast }){U_{g}^{\ast }}^{2}h_{2i}^{\ast },\end{eqnarray}$$

and

(3.48) $$\begin{eqnarray}{\dot{E}}_{p}^{\ast }=\frac{1}{2}\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}^{\ast }\left[(U_{l}^{\ast }-U_{g}^{\ast })(h^{\ast }+h_{1o}^{\ast }-h_{l}^{\ast })\frac{h^{\ast }+h_{1o}^{\ast }+2h_{l}^{\ast }}{3}+(U_{b}^{\ast }-U_{g}^{\ast })\frac{{h_{l}^{\ast }}^{2}-{h_{1i}^{\ast }}^{2}}{2}\right],\end{eqnarray}$$

so that we can evaluate the rate of energy dissipation ${\dot{E}}_{d}^{\ast }$ from (3.32).

5.1 Influence of dimensionless parameters

As stated earlier, in the analytical models presented in § 3, the fluids are assumed to be inviscid and the interfaces are considered sharp, so that no diffusive momentum or mass transfer can occur between the different flow layers. On the other hand, the DNS simulations of § 2 employ finite Reynolds and Péclet numbers, so that momentum and scalar interfaces have finite thicknesses. Here we investigate the quantitative influence of the $Re$ - and $Pe$ -values on the DNS results.

Figure 15.

Influence of the Reynolds number on the downstream density (a), and velocity (b) profiles across the channel height, for simulation S1 shown in figure 4. $Pe$ is held fixed at 20 000. The straight dashed lines represent the model predictions of § 3. As $Re$ increases, the agreement between model predictions and simulation results generally improves.

Figure 16.

Influence of the Péclet number on the downstream density (a), and velocity (b) profiles across the channel height, for simulation S1 shown in figure 4. $Re$ is held fixed at 5000. The straight dashed lines represent the model predictions of § 3. As $Pe$ increases, the agreement between model predictions and simulation results generally improves.

For simulation S1, figures 15 and 16 display the vertical density and velocity profiles far downstream of the gravity current front, at $x^{\ast }=0.9L_{x}$ . Note that in all simulations of figure 15, $Pe=20\,000$ , and in those of figure 16, $Re=5000$ . While the model assumes that the velocity inside the gravity current vanishes in the reference frame moving with the gravity current front, the DNS simulations show that the gravity current fluid moves towards the front near the bottom wall, and away from it near the upper interface, as a result of the finite $Re$ - and $Pe$ -values. The agreement between the model predictions and DNS results improves for larger $Re$ - and $Pe$ -values, as viscous and diffusive effects diminish. Similarly, the velocity gradients in the interfacial regions approach the step-like model profiles for larger $Re$ and $Pe$ . Corresponding observations can be made regarding the density profiles, in that they become more step-like and closer to the model predictions as $Re$ and $Pe$ increase. The influence of $Re$ on the density profiles is relatively small, however. We note that as $Re$ is increased beyond approximately 6000, the agreement between model predictions and DNS results will no longer improve, due to the appearance of unsteady Kelvin–Helmholtz vortices.