# Theoretical model of continuous inertial gravity currents including a jump condition

Safir Haddad, Samuel Vaux, Kevin Varrall, Olivier Vauquelin

► To cite this version:

Safir Haddad, Samuel Vaux, Kevin Varrall, Olivier Vauquelin. Theoretical model of continuous inertial gravity currents including a jump condition. Physical Review Fluids, 2022, 7 (8), pp.084802. 10.1103/PhysRevFluids.7.084802. irsn-04111869

HAL Id: irsn-04111869

https://irsn.hal.science/irsn-04111869

Submitted on 31 May 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Theoretical model of continuous inertial gravity currents
including a jump condition

Safir Haddad ◦, 1, ∗ Samuel Vaux ◦, 1, † Kevin Varrall ◦, 2, ‡ and Olivier Vauquelin ◦, 2, §

¹Institut de Radioprotection et de Sûreté Nucléaire (IRSN), PSN-RES, SA2I, LIE, Cadarache,
13115 Saint-Paul-lez-Durance, France

²Aix-Marseille Université (AMU), Laboratoire IUSTI, CNRS UMR 7343, 5 Rue Enrico Fermi,
13453 Marseille Cedex, France

<figure><img src="image_2.png" /></figure>

(Received 26 December 2021; accepted 25 July 2022; published 19 August 2022)

This paper examines the theoretical modeling of a steady horizontal gravity current involving miscible fluids. The main objective is to determine the longitudinal evolution of the current characteristic quantities: its mean velocity *U*, its height *h*, its mean density *ρ*, and the local Richardson number *Ri* = *g*Δρh/ρ*U*², which characterizes the flow regime. The theory developed by Ellison and Turner [J. Fluid Mech. 6, 423 (1959)] for Boussinesq gravity currents is first extended to the general non-Boussinesq case. In this theoretical approach, the differential equations derived from the conservation equations reveal a mathematical singularity which no longer allows them to be solved when the current passes from a supercritical to a subcritical regime, i.e., when the Richardson number reaches unity. To circumvent this problem, we propose to introduce a jump condition into the model which leads to a sudden transition from a supercritical to a subcritical regime when needed. The jump location is set to satisfy the boundary condition at the exit. Numerical simulations are carried out using a large-eddy simulation code in order to obtain reference results. These results are first used to select a suitable entrainment model among those proposed in the literature. The simulations are then compared with the theoretical model. In the case of a supercritical current without regime change, the agreement between simulation and theory is good. In the case of a supercritical current turning subcritical before the exit, the discontinuity imposed in the model by the jump is clearly abrupt in comparison with the physical reality, but it allows us to reproduce the nonmonotonic evolution of the velocity, height, and Richardson number of the current and to obtain an acceptable estimation of these quantities.

DOI: 10.1103/PhysRevFluids.7.084802

# I. INTRODUCTION

A gravity (or density) current forms when a fluid flows predominantly longitudinally into an ambient fluid of different density. It can develop naturally due to a difference in salinity, temperature, or concentration between the moving fluid and the ambient fluid. In the case of miscible fluids, the current develops by engulfing the surrounding ambient fluid. This mixing process is called entrainment.

Gravity currents occur in many environmental situations such as katabatic winds [1], pyroclastic flows [2], snow avalanches [3], or turbidity currents [4], to name but a few. Gravity currents are also

<footnote>* safir.haddad@irsn.fr</footnote>

<footnote>†samuel.vaux@irsn.fr</footnote>

<footnote>‡kevin.varrall@univ-amu.fr</footnote>

<footnote>§olivier.vauquelin@univ-amu.fr</footnote>
<table><caption>TABLE I. Entrainment laws from the literature.</caption><thead><tr><th>Reference</th><th>Entrainment law</th><th>Range of validity</th></tr></thead><tbody><tr><td>Ellison and Turner [17]</td><td><sup>0.08 &minus; 0.1 Ri<sub>B</sub></sup>&frasl;<sub>1 + 5 Ri<sub>B</sub></sub></td><td>Ri<sub>B</sub> &lt; 0.8</td></tr><tr><td>Lofquist [19]</td><td><sup>0.001</sup>&frasl;<sub>Ri<sub>B</sub></sub></td><td>10 &lt; Ri<sub>B</sub> &lt; 110</td></tr><tr><td>Alpert [23]</td><td>0.12 e<sup>[3.9(Ri<sub>B</sub> inj &minus; Ri<sub>B</sub>)]</sup></td><td>No limit</td></tr><tr><td>Jirka [24]</td><td><sup>k(1 &minus; <sup>Ri<sub>B</sub></sup>&frasl;<sub>√Ri<sub>B</sub><sup>2</sup>+0.25<sup>2</sup></sub>)</sup>&frasl;<sub>(1 + Ri<sub>B</sub>)<sup>-1</sup></sub></td><td>Ri<sub>B</sub> &lt; 1</td></tr><tr><td>Parker et al. [20]</td><td>0.0028 Ri<sub>B</sub><sup>-1.2</sup></td><td>Ri<sub>B</sub> &gt; 0.2</td></tr><tr><td>van Kessel and Kranenburg [21]</td><td><sup>5.5 &times; 10<sup>-3</sup></sup>&frasl;<sub>3.6Ri<sub>B</sub>-1 + √(3.6Ri<sub>B</sub>-1)<sup>2</sup>+0.15</sub></sup></td><td>No limit</td></tr><tr><td>Dallimore et al. [22]</td><td><sup>C<sub>k</sub>c<sub>d</sub><sup>3/2</sup>+C<sub>s</sub><sup>4</sup></sup>&frasl;<sub>Ri<sub>B</sub>+10(C<sub>k</sub>C<sub>d</sub><sup>3/2</sup>+C<sub>s</sub><sup>4</sup>)</sub></td><td>No limit</td></tr><tr><td>Princevac et al. [25]</td><td>0.05 Ri<sub>B</sub><sup>-0.75</sup></td><td>0.15 &gt; Ri<sub>B</sub> &gt; 1.5</td></tr></tbody></table>

flows of interest in safety issues (accidental releases of hazardous materials, fire-induced smoke propagation, etc.). In the particular case of a stratified environment (such as the oceans and the atmosphere), a current can also develop between two layers of different densities. It is then called an intrusion (see, for instance, Refs. [5,6]).

From the pioneering work of von Kármán [7], many authors have investigated this flow by developing theoretical and experimental research [8–14]. Most of these contributions dealt with the transient phase and focused on the advance and shape of the current. An extensive literature can be found in the book by Ungarish [15] and in the recent review of Chowdhury and Testik [16].

In the case of a continuous injection, when the current reaches the exit of the propagation domain, the flow tends to a steady state. It is this configuration that was addressed by Ellison and Turner [17] in their seminal work. These authors developed a theoretical model for a Boussinesq (small density differences involved), weakly turbulent, steady density current on a horizontal or inclined wall. Their equations, based on the conservation of mass, momentum, and buoyancy fluxes, allow the variations of the height, density, and velocity of the current to be calculated along the flow main direction. One of the key parameters of their model is the entrainment coefficient, $E$, which is introduced into the mass conservation equation. This coefficient represents the proportionality between the velocity of the inflow of ambient fluid (into the density current) and the mean velocity of the current. In contrast to plume-type flows [18], in the case of a density current, the entrainment coefficient is not a universal constant but depends on the local stability of the flow, which can be quantified by the Boussinesq Richardson number defined by $Ri_B = \Delta\rho gh/\rho_a U^2$, with $\Delta\rho = |\rho_a - \rho|$, $\rho_a$ the density of the ambient fluid, $\rho$ the density of the current, $g$ the acceleration due to gravity, $h$ the height of the current, and $U$ the velocity of the current. Note that when the inertial forces are predominant in the current ($Ri_B < 1$) the flow is supercritical, when the buoyancy forces dominate the flow is subcritical ($Ri_B > 1$), and when these forces are balanced ($Ri_B = 1$) the flow is considered to be critical. From experiments involving salt water injected into fresh water, Ellison and Turner [17] found that $E \propto Ri_B^{-1}$ up to the value $Ri_B = 0.8$, beyond which the entrainment vanishes. This supposes that only a supercritical current is expected to entrain ambient fluid.

Subsequently, many of the experimental studies carried out have aimed at refining this entrainment submodel. We can quote among others Refs. [19–22] in hydraulics, [23,24] for fire safety issues, and [25] for atmospheric flows. The corresponding entrainment laws and their range of validity are summarized in Table I. We refer the reader to the papers of Christodoulou [26] and Fernando [27] for a more extensive review.
<figure><img src="image_3.png" /><figcaption>FIG. 1. Schematic diagram of the density current.</figcaption></figure>

While research into the measuring and modeling of the entrainment is extensive, we note that the resolution of the governing equations proposed by Ellison and Turner [17] has received much less attention. Indeed, as we shall see in the next section, depending on the length of the studied domain, these governing equations may reveal a mathematical singularity which no longer allows them to be solved when the flow reaches the critical state. Recently, Guo et al. [28] tackled the problem in the case of an initially supercritical gravity current. They solved the equations of Ellison and Turner [17] up to the singular value ($Ri_B = 1$) and then considered that the flow remains critical up to the exit of the domain. Their hypothesis requires the artificial modification of one variable of the current (i.e., one of either the height, velocity, or density) and unfortunately fails to conserve simultaneously buoyancy, mass, and momentum fluxes.

In this paper, we propose a resolution of the equations of Ellison and Turner [17] for an initially supercritical gravity current by introducing a discontinuity similar to a jump as done in Refs. [29–31] for free-surface hydraulic flows. This approach allows the fluxes to be conserved and the nonmonotonic behavior of the velocity and height of the current to be reproduced. In order to tackle problems involving large density contrasts, calculations will be done in the general non-Boussinesq case.

The paper is organized as follows: Sec. II deals with non-Boussinesq theoretical equations for the steady gravity current. Section III focuses on the large-eddy simulation setup. Section IV presents observations from the numerical simulations and comparisons between theoretical and numerical results. Conclusions are drawn in Sec. V.

## II. THEORY

### A. Governing equations

We consider a steady two-dimensional gravity current developing along a rigid horizontal wall of length $L$ as represented in Fig. 1. The current is released with a horizontal velocity $U_i$ and a density $\rho_i$ from a two-dimensional opening of height $h_i$. The ambient fluid, of density $\rho_a > \rho_i$, is at rest. Along the longitudinal $x$ axis, the gravity current cross-section scales of velocity, height, and density are, respectively, denoted $U(x)$, $h(x)$, and $\rho(x)$ and are obtained using the following integral formulations for the mass, volume, and momentum fluxes per unit width:

$$ \rho(x) U(x) h(x) = \int_{0}^{\infty} \rho(x, z) u(x, z) dz, \qquad (1) $$

$$ U(x) h(x) = \int_{0}^{\infty} u(x, z) dz, \qquad (2) $$

$$ \rho(x) U(x)^2 h(x) = \int_{0}^{\infty} \rho(x, z) u(x, z)^2 dz, \qquad (3) $$

where $u(x, z)$ is the local velocity of the layer and $\rho(x, z)$ its local density. For the sake of clarity, in the following, $\rho(x)$, $U(x)$, and $h(x)$ will be denoted $\rho$, $U$, and $h$, respectively. Note that a momentum
correction factor should multiply the left-hand side of Eq. (3) [32,33]. Nevertheless, as mentioned by Gu and Lawrence [34], for the type of flow we are dealing with, the correction is small and can be neglected. Consequently, for the sake of simplicity this correction factor will be set to unity.

The governing equations are established in the general non-Boussinesq case, by considering the mass, momentum, and buoyancy conservation over a horizontal infinitesimal element of length $dx$ of the gravity current. These equations are as follows:

$$ \frac{d(\rho U h)}{dx} = E U \rho_a, \qquad (4) $$

$$ \frac{d(\rho U^2 h)}{dx} = -C_D U^2 \rho - \frac{1}{2} \frac{d(\Delta\rho g h^2)}{dx}, \qquad (5) $$

$$ \frac{d}{dx} \left( \frac{\Delta \rho}{\rho_a} g U h \right) = 0, \qquad (6) $$

where $E$ is the entrainment coefficient and $C_D$ is the drag coefficient. The terms on the right-hand side of Eq. (5) represent, respectively, the turbulent frictional drag on the wall and the pressure force associated with the change of height and density of the layer. Equation (6) supposes a constant buoyancy flux. This is no longer accurate when heat transfers occur between the wall and the current (a heat sink term would appear in the right-hand side). We will consider only the case of an adiabatic flow in the following.

Equations (4), (5), and (6) can be combined so that the first derivatives $dh/dx$, $dU/dx$, and $d\rho/dx$ can be expressed as

$$ \frac{dh}{dx} = \frac{(1 + \frac{\rho_a}{\rho} - \frac{1}{2} \mathrm{Ri})E + C_D}{(1 - \mathrm{Ri})}, \qquad (7) $$

$$ \frac{dU}{dx} = -\frac{U \left( \frac{\rho_a}{\rho} + \frac{1}{2} \mathrm{Ri} \right) E + C_D}{h (1 - \mathrm{Ri})}, \qquad (8) $$

$$ \frac{d\rho}{dx} = \frac{\Delta\rho E}{h}, \qquad (9) $$

where $\mathrm{Ri}$ is the non-Boussinesq Richardson number defined as

$$ \mathrm{Ri} = \frac{\Delta\rho gh}{\rho U^2}. \qquad (10) $$

By combining (7), (8), (9), and (10), we obtain the expression for the first derivative of $\mathrm{Ri}$:

$$ \frac{d\mathrm{Ri}}{dx} = \frac{\mathrm{Ri} \left(1 + 2\frac{\rho_a}{\rho}\right) \left(1 + \frac{\mathrm{Ri}}{2}\right) E + 3C_D}{h (1 - \mathrm{Ri})}. \qquad (11) $$

Note that Eqs. (7), (8), (9), and (11) are an extension of those of Ellison and Turner [17] in the general non-Boussinesq case.

This model requires the knowledge of two parameters: the entrainment coefficient $E$ and the drag coefficient $C_D$. The choice of the law for the entrainment coefficient will be discussed in Sec. IV B. Concerning the drag coefficient, since there is no consensus in the literature, we have decided to choose the value used by Kunsch [35] for a model of smoke backflow in a tunnel: $C_D = 0.0065$ (as in Refs. [36–38]). In order to estimate the influence of this parameter, a sensitivity study has been performed in the case of a supercritical current. It appears that doubling the value of $C_D$ leads to relative variations less than 5% for the longitudinal evolution of the primary variables ($U$, $h$, and $\rho$).

Through Eqs. (7), (8), and (11), the theoretical model exhibits a mathematical singularity when $\mathrm{Ri} = 1$. For an initially supercritical gravity current ($\mathrm{Ri}(x = 0) < 1$), the right-hand side in Eq. (11) is always positive, which implies that $\mathrm{Ri}(x)$ is a monotonic increasing function. From there, two cases are to be considered: (i) the flow remains supercritical ($\mathrm{Ri} < 1$) over the whole domain and
<figure><img src="image_4.png" /><figcaption>FIG. 2. Schematic diagram of a density jump.</figcaption></figure>

Eqs. (7), (8), (9), and (11) can be solved directly and (ii) the gravity current transitions from a supercritical flow ($Ri < 1$) to a subcritical flow ($Ri > 1$) through the critical state ($Ri = 1$), which is theoretically problematic. In order to circumvent this problem, we then introduce a mathematical discontinuity similar to a jump.

## B. Density jump

A hydraulic jump is a well-known phenomenon in free-surface flows in hydraulics. It is characterized by a sudden drop in velocity and an abrupt increase in the thickness of the fluid layer flowing along a wall. In the case of miscible fluids, this phenomenon is referred to as a density jump [39].

As shown in Fig. 2, using respectively the subscripts 1 and 2 for the quantities upstream and downstream of the density jump, we write the conservation equations for the mass, momentum, and buoyancy, similarly to Yih and Guha [40]:

$$ \rho_1 U_1 h_1 + \rho_a \epsilon U_1 h_1 = \rho_2 U_2 h_2, \quad (12) $$

$$ \rho_1 U_1^2 h_1 + \frac{\Delta \rho_1 g h_1}{2} = \rho_2 U_2^2 h_2 + \frac{\Delta \rho_2 g h_2}{2}, \quad (13) $$

$$ \Delta \rho_1 U_1 h_1 = \Delta \rho_2 U_2 h_2, \quad (14) $$

with $\epsilon$ being the ratio of entraining mass flux to the upstream mass flux [41–43]), $\Delta\rho_1 = \rho_a - \rho_1$ and $\Delta\rho_2 = \rho_a - \rho_2$.

Combining (12), (13), and (14), we obtain an equation for the ratio $h_2/h_1$ allowing the downstream height $h_2$ to be expressed as a function of the upstream conditions:

$$ \left(\frac{h_2}{h_1}\right)^3 - \left(\frac{h_2}{h_1}\right) \frac{(2 + \mathrm{Ri}_1)}{\gamma \mathrm{Ri}_1} + \frac{2(1 + \epsilon)}{\gamma^2 \mathrm{Ri}_1} = 0, \quad (15) $$

with $\gamma = \frac{\rho_a}{\rho_a + \epsilon\rho_1}$. Solutions of (15) are given in Regev et al. [42]. In the particular case where $\epsilon = 0$ (no entrainment into the jump), the resolution of (15) gives the well-known equation of Bélanger:

$$ \frac{h_2}{h_1} = \frac{\sqrt{1 + \frac{8}{\mathrm{Ri}_1}} - 1}{2}. \quad (16) $$

Given that $\epsilon = 0$, $\Delta\rho_1 = \Delta\rho_2$ and Eq. (12) or (14) gives trivially $U_2/U_1 = h_1/h_2$. Note that the jump being a complex phenomenon, it is unlikely that no entrainment occurs in it, and considering $\epsilon = 0$ is a strong assumption. Unfortunately, this issue is poorly documented in the literature. Consequently, in the absence of reliable laws or correlations, and to provide an uncluttered model,
we have chosen to consider no entrainment within the jump in the theoretical model. If this had not been the case, it would have been physically meaningful to also consider a finite length for the jump, rather than a simple discontinuity. These issues will be discussed in Sec. V.

Finally, by using the ratios of height, velocity, and density, the ratio of Richardson numbers between the upstream and downstream positions is given by

$$ \frac{Ri_2}{Ri_1} = \frac{1}{8} \left( \sqrt{1 + \frac{8}{Ri_1}} - 1 \right)^3. \qquad (17) $$

This relation shows that for a supercritical upstream flow ($Ri_1 < 1$), the flow downstream of the jump is subcritical ($Ri_2 > 1$), and the evolution of the Richardson number is then monotonically decreasing [see Eq. (11)]. In addition, at the exit of domain (similar to a step of infinite height), as explained by Henderson [33], the flow must be close to the critical condition, i.e., $Ri(x = L) \approx 1$. This condition is justified in the Appendix. In our theoretical approach, this boundary condition (at the exit of the domain) allows the mathematical discontinuity (i.e., the jump) to be located on the x axis.

In summary, the theoretical model developed above allows us to calculate the characteristics (velocity, height, and density) of a non-Boussinesq density current flowing along a horizontal wall of length L. In the event of a transition from supercritical to subcritical flow, this model includes a mathematical discontinuity (jump) in order to respect not only the boundary condition at the exit, but also the conservation of all the fluxes.

### C. Theoretical resolution

To solve the set of coupled ordinary differential Eqs. (7), (8), and (9) for initially supercritical gravity currents, we use the method of Dormand and Prince [44].

As previously mentioned, two situations can occur. In the first, the Richardson number remains less than one until the exit of the domain, and Eqs. (7), (8), and (9) are numerically solved without any difficulty. In this case, the flow remains supercritical.

In the second situation, the Richardson number reaches the critical value ($Ri = 1$) before $x = L$ (point 1 in Fig. 3). The equations can no longer be solved because of the singularity. Practically, from here we stop the resolution of (7), (8), and (9) before this singularity and then apply a jump according to Eqs. (15) and (17). Downstream of the jump, we obtain values which are the initial conditions to be used in (7), (8), and (9) for the subcritical flow. If the jump is located too far upstream (point 2 in Fig. 3), the Richardson number at the exit is higher than unity, which is not in agreement with the boundary condition. Conversely, if the jump is not located sufficiently upstream (point 3 in Fig. 3), the Richardson number reaches unity before the exit of the domain, which is unacceptable physically. In practice, we proceed by successive iterations until we find the abscissa (point 4 in Fig. 3), which allows us to catch the critical regime ($Ri = 1$) at the exit of the domain.

## III. NUMERICAL SETUP

We consider horizontal, isothermal, and continuous releases of light gases to simulate non-Boussinesq turbulent miscible gravity currents. Simulations have been carried out for a density ratio $\rho_i/\rho_a$ ranging from 0.5 to 0.83. A total of four simulations were carried out over the range 0.003 $\lesssim$ $Ri \lesssim$ 0.112 corresponding originally to the supercritical non-Boussinesq regime. In all simulations, the Reynolds number based on the source height is set to $Re > 20000$, where $Re = U_i h_i / \nu$ (see Table II for the injection parameters of the continuous gravity current simulations) and $\nu$ is the kinematic viscosity.

In order to perform the numerical simulation of an inertial gravity current, large-eddy simulations (LES) are used to solve the Favre-filtered Navier-Stokes equations (mass and momentum balance) along with species transport equations. We use the numerical computational code CALIF³S-ISIS [soft-
<figure><img src="image_5.png" /><figcaption>FIG. 3. Graphical representation of the calculation method in the case of a density jump. The solid line represents the first calculation until Ri reaches the unity (point 1). The two black dashed lines represent two calculations (points 2 and 3) not in agreement with the hypothesis of a critical flow at the exit of the domain. The red dashed line (point 4) represents the final solution giving Ri = 1 at the exit of the domain.</figcaption></figure>

ware developed at the French Institut de Radioprotection et de Sûreté Nucléaire (IRSN)], dedicated
to three-dimensional simulations of turbulent and slightly compressible flows (low-Mach-number
approach). In Cartesian coordinates, the three-dimensional filtered Navier-Stokes equations are as
follows:

$$
\frac{\partial \bar{\rho}}{\partial t} + \frac{\partial (\bar{\rho} \tilde{u}_i)}{\partial x_i} = 0, \quad (18)
$$

$$
\frac{\partial (\bar{\rho} \tilde{u}_i)}{\partial t} + \frac{\partial (\bar{\rho} \tilde{u}_i \tilde{u}_j)}{\partial x_j} = - \frac{\partial \bar{p}}{\partial x_i} + \frac{\partial \bar{S}_{ij}}{\partial x_j} + (\rho_a - \bar{\rho}) g_i - \frac{\partial \tau_{ij}}{\partial x_j}, \quad (19)
$$

where $\tilde{u}_i$ is the Favre-filtered velocity and $\bar{p}$ is the dynamic pressure. The density $\bar{\rho}$ is the filtered
density of the fluid and is determined by the ideal gas law in combination with the mass frac-
tion of the different species of the gas mixture. In Eq. (19), $g_i$ is the gravitational acceleration,
$\tau_{ij} = \bar{\rho}u_i u_j - \bar{\rho}\tilde{u}_i \tilde{u}_j$ denotes the subgrid-scale Reynolds stress, and $\bar{S}_{ij} = -(2/3)\mu(\partial\tilde{u}_k/\partial x_k)\delta_{ij} +$
$\mu(\partial\tilde{u}_i/\partial x_j + \partial\tilde{u}_j/\partial x_i)$ represents the filtered strain rate tensor, where $\mu$ is the dynamic molecular
viscosity calculated as a function of the individual viscosities and molar masses as well as the
corresponding mass fractions.

<table><caption>TABLE II. Source parameters of continuous gravity current simulations and length of the domain.</caption>
  <thead>
    <tr>
      <th></th>
      <th><i>U</i><sub>i</sub> (m/s)</th>
      <th><i>h</i><sub>i</sub> (m)</th>
      <th><i>ρ</i><sub>i</sub> (kg/m<sup>3</sup>)</th>
      <th><i>Ri</i><sub>i</sub></th>
      <th><i>Re</i><sub>i</sub></th>
      <th><i>L</i> (m)</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <td>Case 1</td>
      <td>8</td>
      <td>0.1</td>
      <td>1</td>
      <td>0.003</td>
      <td>43 260</td>
      <td>10</td>
    </tr>
    <tr>
      <td>Case 2</td>
      <td>4.2</td>
      <td>0.2</td>
      <td>0.6</td>
      <td>0.112</td>
      <td>22 401</td>
      <td>10</td>
    </tr>
    <tr>
      <td>Case 3</td>
      <td>5</td>
      <td>0.5</td>
      <td>0.75</td>
      <td>0.073</td>
      <td>100 931</td>
      <td>50</td>
    </tr>
    <tr>
      <td>Case 4</td>
      <td>10</td>
      <td>0.5</td>
      <td>0.75</td>
      <td>0.018</td>
      <td>201 862</td>
      <td>50</td>
    </tr>
  </tbody>
</table>
For each species k, the mass fraction $y_k$ is governed by a species transport equation which reads as

$$ \frac{\partial \bar{\rho} \tilde{y}_k}{\partial t} + \frac{\partial (\bar{\rho} \tilde{y}_k \tilde{u}_i)}{\partial x_i} = \frac{\partial}{\partial x_i} \left( \bar{\rho} D \frac{\partial \tilde{y}_k}{\partial x_i} + \frac{\mu_t}{Sc_t} \frac{\partial \tilde{y}_k}{\partial x_i} \right), \quad (20) $$

where $\tilde{y}_k$ represents the Favre-filtered mass fraction of the kth component of the mixture and D stands for the molecular diffusivity of the mixture. In Eq. (20), we employ the simple gradient diffusion hypothesis to close the problem with a turbulent Schmidt number $Sc_t$ set to 0.7.

The present large-eddy simulations implicitly apply a box filter in each direction and the wall adapting local eddy subgrid-scale model for the subgrid Reynolds stress [45] is adopted. We use a staggered grid with a cell-centered piecewise constant representation of the scalar variables and a marker-and-cell-type finite volume approximation for the velocity. For the time discretization, we employ a fractional step algorithm decoupling balance equations for the transport of species and Navier-Stokes equations, which are solved by a pressure correction technique. Since we consider gravity currents in an infinite (open) environment, the computational domain must be bounded by artificial boundary conditions which perturb as little as possible the flow in the interior of the domain. In our simulations, the boundary conditions used are based on the usual control of the kinetic energy and allow us to distinguish between the flow that enters the domain and the flow that leaves. This type of boundary condition was originally established for the incompressible case in Bruneau and Fabrie [46,47], and its extension to compressible flows was tackled in Bruneau [48]. Periodic boundary conditions are imposed in the spanwise direction.

The three-dimensional computational domain $\Omega$ is a rectangular box of dimensions $L_x \times L_y \times L_z$. For each simulation, the value set to the horizontal streamwise length is $L_x = L + L_i$ (where $L_i$ is the horizontal length from the left boundary of the domain up to the planar injection source) and the vertical length $L_z$ depends on the physical parameters of the flow under study. The planar injection source is positioned vertically (of height $h_i$) slightly downstream of the left boundary of the domain and in contact with the top solid boundary of the computational domain $\Omega$. The injection source is taken as the origin of the longitudinal $x$ axis. On this inflow boundary, the flow emerges horizontally with a uniform velocity profile $U_i$. We use a refined Cartesian grid with a uniform rectangular mesh ($\Delta x \times \Delta y$) from the source horizontal position $x_i$ to the outlet located at $L_x$. From the source horizontal position $x_i$ to the left boundary of the domain ($x = -1$ m), the grid is horizontally stretched. In the vertical direction $z$, the grid spacing ($\Delta z$) is kept uniform from $z = 0$ to $z = h_i$ (with $\Delta z = \Delta z_1$), still uniform from $z = h_i$ up to a vertical distance $L_{1z}$ (with $\Delta z = \Delta z_2$) and then stretched toward the opposite boundary (i.e. the bottom of the domain). In order to initiate the turbulence at the source, we apply an azimuthal forcing similar to Zhou, Luo, and Williams [49].

For each simulated case, a grid-convergence study was carried out to validate the domain height $L_z$, extent of the vertical subregion $L_{1z}$, spanwise width $L_y$, and grid spacing in each direction. The choice of $L_{1z}$ was made to approximately demarcate the thickness of the gravity current. The spanwise width of the domain $L_y$ ranges from $4h_i$ to $6h_i$. We tested vertical grid spacings $\Delta z_1/h_i$ ranging from 0.067 to 0.04, $\Delta z_2/h_i$ ranging from 0.24 to 0.05, horizontal grid spacings $\Delta x/h_i$ varying from 0.15 to 0.075, and spanwise grid spacings $\Delta y/h_i$ from 0.1 to 0.05. Concerning the time discretization, a CFL (Courant-Friedrichs-Lewy) number close to unity has been imposed for each calculation even if time step sizes for which CFL numbers greater than one are allowed when using implicit schemes.

As an illustration of the flow obtained by these simulations, Fig. 4 presents an instantaneous density field of a gravity current developing along a wall until the exit of the domain. During the flow development, the layer thickens and engulfs surrounding fluid through large eddies.

The duration times of the simulations were set sufficiently large to ensure first that the steady states of the currents are reached and second to guarantee the convergence of the time-averaged values of the current variables (the variations of the mean field fall below 2% of the value of the mean). For instance, in case 1 (see Table II), we set the duration of the numerical simulation to $t = 250$ s. Given the flow initial conditions and the domain length, the transient phase lasts for a
<figure><img src="image_6.png" /><figcaption>FIG. 4. Instantaneous density field (case 1).</figcaption></figure>

short duration of a few seconds. The statistics are then calculated from this time until the end of
the simulation (i.e., for the first case, from 15 to 250 s). To illustrate the temporal convergence,
the longitudinal evolution of the time-averaged cross-sectional mass flow rate per unit width for
different durations is plotted in Fig. 5. We can see that the statistics are converged for $t = 200$ s and
are consequently fully converged at $t = 250$ s.

Furthermore, previous simulations presented in Refs. [50,51] compared the aforementioned LES approach on turbulent miscible Boussinesq and non-Boussinesq flows with experimental data. These papers confirm the suitability of the CALIF³S-ISIS code to properly evaluate the behaviors of turbulent buoyant flows exhibiting large density differences.

## IV. RESULTS

### A. Flow observations

For the four cases presented in Table II, Fig. 6 represents the time-averaged $x$ velocity (left column) and the time-averaged density (right column) fields provided by the large-eddy simulations. In the analysis of these fields, we focus only on the domain between the injection position and the vertical dotted lines, in order to exclude the zone of the flow influenced by the weir at the exit. The location of these dotted lines will be discussed and justified later.

From the $x$-velocity fields, we can see two distinct behaviors. For cases 2, 3, and 4 represented in Figs. 6(c), 6(e), and 6(g), we observe a velocity decrease in the first meters after the release, followed by a velocity increase until the exit of the studied domain. This nonmonotonic behavior of the velocity suggests the existence of the two regimes (supercritical and then, subcritical) and thus of

<figure><img src="image_7.png" /><figcaption>FIG. 5. Longitudinal evolution of the mass flow rate per unit width at different times. Dashed line corresponds to $t = 50$ s, dotted lines corresponds to $t = 200$ s, and solid line corresponds to $t = 250$ s.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>FIG. 6. Time-averaged *x*-velocity field for case 1 (a), case 2 (c), case 3 (e), and case 4 (g) and time averaged density field for case 1 (b), case 2 (d), case 3 (f), and case 4 (h).</figcaption></figure>

a critical transition. In contrast, for case 1 depicted in Fig. 6(a), the velocity decreases monotonically until the exit of the studied domain, which is characteristic of a regime remaining supercritical.

We also observe two distinct behaviors from the density fields. For case 1, Fig. 6(b) indicates that the layer thickens monotonically over the entire studied domain. In contrast, for the other cases (2, 3, and 4) represented in Figs. 6(d), 6(f), and 6(h), we see a nonmonotonic behavior with a current thickening and then thinning after a location, which probably corresponds to a critical transition.

## B. Entrainment law

Before going further with the resolution of Eqs. (7)–(9) and their comparisons with the numerical simulations, we need to choose the entrainment law among those available in the literature. We obtain it by using case 1 for the sake of simplicity since the behavior of the Richardson number
<figure><img src="image_9.png" /><figcaption>FIG. 7. Comparison between different entrainment laws. Circle symbols (○) represent the entrainment determined using Eq. (21). The other symbols represent laws found in the literature with the law of Ellison and Turner [17] [triangle symbols (▷)], the law of van Kessel and Kranenburg [21] [dashed curve with cross symbols (+)], the law of Jirka [24] [star symbols (∗)], and the law of Princevac et al. [25] [square symbols (□)].</figcaption></figure>

seems monotonic without regime change. The velocity and density vertical profiles obtained in the numerical simulation are integrated according to Eqs. (1)–(3), and the entrainment coefficient is calculated along the x axis with the following relation:

$$E = \frac{d(\rho U h)}{dx} \frac{1}{\rho_a U}. \qquad (21)$$

In Fig. 7 the calculated entrainment coefficient is compared with four laws found in the literature. In our case, three laws [17,21,24] are in good agreement with our data. The fourth one [25] gives for its part a much higher value. This difference is caused by the range of Reynolds number used in their work (~10⁷). For this range of Reynolds number, according to Princevac et al. [25], the entrainment coefficient increases considerably. Hereafter, we use the law of van Kessel and Kranenburg [21] since it was established from theoretical analysis and not constrained by a restricted Richardson-number range.

### C. Comparison between numerical results and the theoretical model

In order to compare the theoretical model predictions with the numerical simulations, the numerical data are first integrated along the vertical z axis using relations (1)–(3). Comparisons are carried out for the four cases presented in Table II. For each case, we compare the longitudinal evolution of the height, velocity, density, and Richardson number of the current. As mentioned previously in Sec. IV A, our comparisons are made on the entire domain, except the last portion of the domain located beyond the dotted lines in Fig. 6. The choice of the location of these dotted lines is based on the numerical data and is inherent to the regime of the flow:

(1) For cases requiring a jump, since the theoretical model needs the critical condition to be reached at the end (see the Appendix), we set the comparison domain between the injection and the abscissa where the simulated flow reaches Ri = 1 slightly upstream of the weir.

(2) For cases not requiring a jump, the theoretical model allows only monotonic evolution to be reproduced. In contrast, the simulated flow shows a nonmonotonic behavior due to the weir influence before the exit. Consequently, we have set the comparison domain between the injection
<figure><img src="image_10.png" /><figcaption>FIG. 8. Longitudinal evolution of (a) height, (b) velocity, (c) density, and (d) Richardson number of the current for case 1. Comparisons between theoretical model (solid line) and numerical data from the LES simulation (circle symbols).</figcaption></figure>

and the abscissa where the simulated flow changes its behavior, corresponding to the location where $\frac{dRi}{dx} \le 0$.

These features lead us to consider about 80%–90% of the domain length L to compare theoretical and numerical data.

Case 1 corresponds to a turbulent, Boussinesq, highly inertial current that develops along a 10-m-long horizontal wall. In this case, the current remains supercritical from the source release to the exit of the domain. All the characteristic quantities (height, velocity, and density) evolve in a monotonic way, as well as the Richardson number, which increases while remaining much lower than the critical value ($Ri = 1$). The comparison between the numerical data and the model is presented in Fig. 8. A good agreement is observed with only moderate deviations. For this simulation, since the critical value is not reached, the theoretical model can be solved directly on the whole domain.

Case 2 corresponds to a turbulent, non-Boussinesq, inertial current that develops along a 10-m-long horizontal wall. Compared to case 1, the source Richardson number for case 2 is much closer to unity in order to reach the critical regime before the exit of the domain. As can be seen in Fig. 9, the simulation shows a nonmonotonic behavior for the velocity, the height and the Richardson number, which reaches the critical value ($Ri = 1$) at about 20% of the total domain length ($0.2L$).
<figure><img src="image_11.png" /><figcaption>FIG. 9. Longitudinal evolution of (a) height, (b) velocity, (c) density, and (d) Richardson number of the current for case 2. Comparisons between theoretical model (solid line) and numerical data from the LES simulation (circle symbols).</figcaption></figure>

The Richardson number goes through a maximum value in the subcritical regime (at about 0.3L) and then decreases towards unity.

To account for this nonmonotonic behavior, the model triggers a jump that is found to be located at the distance 0.16L from the source, as observed in Fig. 9. This jump is clearly abrupt, but it allows us to match satisfactorily the upstream (supercritical) and downstream (subcritical) evolution of the velocity, height, and the Richardson number. The monotonic behavior of the density is also properly reflected by the model. After the jump, in the subcritical regime, the density increases slowly, which is qualitatively well reproduced by the theoretical model. Nevertheless, we notice a gap for the density between the model and the simulation. One possibility to reduce this gap would be to introduce a specific entrainment into the jump [$\epsilon \neq 0$ in (15)], as will be discussed in Sec. V.

For cases 3 and 4, the length of the domain is increased by a factor of 5 (i.e., 50 m). The current is still turbulent, non-Boussinesq and inertially dominated, but the strategy to compare theory and simulation is different from the two previous cases. Indeed, the theoretical model source conditions are chosen to produce a jump close to the injection (case 3) and a jump much further away (case 4). The numerical simulations are then compared with the results predicted by the model as shown in Figs. 10 and 11. Again, the evolution of the characteristic quantities is relatively well predicted with the theoretical model including a jump.
<figure><img src="image_12.png" /><figcaption>FIG. 10. Longitudinal evolution of (a) height, (b) velocity, (c) density, and (d) Richardson number of the current for case 3. Comparisons between theoretical model (solid line) and numerical data from the LES simulation (circle symbols).</figcaption></figure>

V. CONCLUSIONS

We have investigated the dynamics of non-Boussinesq steady-state turbulent inertial gravity currents by numerical simulations and a theoretical approach similar to that developed in the seminal work of Ellison and Turner [17]. In the theoretical approach, the equations obtained present a mathematical singularity when the Richardson number is equal to 1 (critical condition). This does not allow the transition between a supercritical and a subcritical regime to be represented. In order to overcome this problem, and to reproduce the nonmonotonic behaviors observed in the numerical simulations, we have introduced a mathematical discontinuity in the theoretical model which is similar to a density jump. The amplitude of this jump is given by the Bélanger equation, and its location is determined from the exit condition (weir), which imposes a critical condition to be reached.

In the case where the flow remains supercritical over the entire domain, the equations can be directly solved since the mathematical singularity never appears in this situation. No jump is then required. A good agreement between the theoretical and numerical data is observed for the primary variables of the current (velocity, density, and thickness).

In the case where the flow transitions from the supercritical ($Ri < 1$) to the subcritical ($Ri > 1$) state before the end of the domain, the introduction of the jump in the theoretical model allows
<figure><img src="image_13.png" /><figcaption>FIG. 11. Longitudinal evolution of (a) height, (b) velocity, (c) density, and (d) Richardson number of the current for case 4. Comparisons between theoretical model (solid line) and numerical data from the LES simulation (circle symbols).</figcaption></figure>

the nonmonotonic behaviors of the current thickness and velocity, as well as that of the Richardson number, to be reproduced in a qualitatively satisfactory way. Obviously, close to the jump, large deviations are observed between the theory and the numerical simulations. However, further away from this zone, the agreement becomes acceptable.

As a possible improvement in the theoretical model, we could consider a local entrainment into the jump region [39] as well as a finite length [52] for this region of transition between the supercritical and subcritical regimes. As an example, in Fig. 12 we present the evolution of the Richardson number for the numerical simulation, the initial model including an abrupt jump without entrainment and the model modified to include an entrainment into a jump of finite length. For the improved model, we have set the entrainment coefficient to $\epsilon = 0.1$, which is a realistic order of magnitude considering the work of Regev et al. [41], and a length of $L_j = 6.1h_2$, proposed by Henderson [33]. As can be seen in the figure, the modifications brought to the theoretical model improve its relevance. Nevertheless, the introduced parameters (entrainment coefficient $\epsilon$ and jump length $L_j$) cannot be chosen in an universal way due to a lack of reliable data and theoretical works in the literature.

On the other hand, from an energetic point of view, the vortex dynamics developing into a density current has an influence and should be taken into account in a theoretical approach, especially in
<figure><img src="image_14.png" /><figcaption>FIG. 12. Longitudinal evolution of the Richardson number of the current for case 3 (see Table II). Solid line (—) is the initial model. Dashed line (- -) is the improved model. Circle symbols (o) are the LES simulation.</figcaption></figure>

the momentum equation. Among the most relevant studies that have addressed this issue, we can mention Hornung et al. [53] and Richard and Gavrilyuk [54], who use the entrophy as a source term. Moreover, it would also be interesting to consider the potentialities of the circulation-based model developed in Borden and Meiburg [55] or Khodkar et al. [56] in order to improve the theoretical approach presented in this article.

APPENDIX: WEIR BOUNDARY CONDITION

As illustrated in Fig. 13, we consider a steady two-dimensional gravity current reaching the exit of a horizontal wall. The two points A and B, are located just upstream of the end and separated by a short distance dx. Given that dx is small, we assume that the friction and the mixing between A and B can be neglected. Thus, we write the Bernoulli formula between these two points as follows:

$$P_A + \frac{1}{2} \rho U_A^2 + \rho g h_A = P_B + \frac{1}{2} \rho U_B^2 + \rho g h_B. \quad (A1)$$

The relation between static pressures $P_A$ and $P_B$ is obtained by hydrostatic considerations out of the flow:

$$P_A = P_B + \rho_0 g (h_A - h_B). \quad (A2)$$

By combining Eqs. (A1) and (A2) and assuming that dx is an infinitesimal distance, we obtain

$$\frac{d}{dx} \left( \Delta \rho g h + \frac{1}{2} \rho u^2 \right) = 0. \quad (A3)$$

<figure><img src="image_15.png" /><figcaption>FIG. 13. Sketch of the flow at the exit of the domain.</figcaption></figure>
By introducing the volume flow rate $q = uh$, (A3) can be rewritten as

$$ \frac{dh}{dx} \frac{d}{dh} \left( \Delta \rho gh + \frac{\rho q^2}{2h^2} \right) = 0. \qquad (A4) $$

After some algebraic manipulations it comes that

$$ \frac{dh}{dx}(Ri - 1) = 0. \qquad (A5) $$

Assuming that point B is influenced by the weir boundary condition at the exit, contrary to point A, the flow has a smaller thickness at point B than at point A (i.e. $h_B < h_A$ and so $dh/dx \neq 0$). Therefore, the Richardson must be unitary at point B to comply with Eq. (A5), and consequently $Ri(x = B) = 1$. This was also observed by Bauer and Graf [57] and Graf and Altinakar [58].

[1] D. R. Fitzjarrald, Katabatic wind in opposing flow, J. Atmos. Sci. 41, 1143 (1984).

[2] G. A. Valentine, Stratified flow in pyroclastic surges, Bull. Volcanol. 49, 616 (1987).

[3] E. J. Hopfinger, Snow avalanche motion and related phenomena, Annu. Rev. Fluid Mech. 15, 47 (1983).

[4] E. Meiburg and B. Kneller, Turbidity currents and their deposits, Annu. Rev. Fluid Mech. 42, 135 (2010).

[5] M. R. Flynn and P. F. Linden, Intrusive gravity currents, J. Fluid Mech. 568, 193 (2006).

[6] R. E. Britter and J. E. Simpson, A note on the structure of the head of an intrusive gravity current, J. Fluid Mech. 112, 459 (1981).

[7] T. von Kármán, The engineer grapples with nonlinear problems, Bull. Amer. Math. Soc. 46, 615 (1940).

[8] T. B. Benjamin, Gravity currents and related phenomena, J. Fluid Mech. 31, 209 (1968).

[9] D. Hoult, Oil spreading on the sea, Annu. Rev. Fluid Mech. 4, 341 (1972).

[10] J. E. Simpson and R. E. Britter, The dynamics of the head of a gravity current advancing over a horizontal surface, J. Fluid Mech. 94, 477 (1979).

[11] H. E. Huppert, The propagation of two-dimensional and axisymmetric viscous gravity currents at a fluid interface, J. Fluid Mech. 121, 43 (1982).

[12] A. Dai, C. E. Ozdemir, M. I. Cantero, and S. Balachandar, Gravity currents from instantaneous sources down a slope, J. Hydraul. Eng. 138, 237 (2012).

[13] Z. Borden, E. Meiburg, and G. Constantinescu, Internal bores: An improved model via a detailed analysis of the energy budget, J. Fluid Mech. 703, 279 (2012).

[14] M. A. Khodkar, M. M. Nasr-Azadani, and E. Meiburg, Partial-depth lock-release flows, Phys. Rev. Fluids 2, 064802 (2017).

[15] M. Ungarish, An Introduction to Gravity Currents and Intrusions (CRC Press, New York, 2009).

[16] M. R. Chowdhury and F. Y. Testik, A review of gravity currents formed by submerged single-port discharges in inland and coastal waters, Environ. Fluid Mech. 14, 265 (2014).

[17] T. H. Ellison and J. S. Turner, Turbulent entrainment in stratified flows, J. Fluid Mech. 6, 423 (1959).

[18] B. R. Morton, G. Taylor, and J. S. Turner, Turbulent gravitational convection from maintained and instantaneous sources, Proc. R. Soc. London A 234, 1 (1956).

[19] K. Lofquist, Flow and stress near an interface between stratified liquids, Phys. Fluids 3, 158 (1960).

[20] G. Parker, M. Garcia, Y. Fukushima, and W. Yu, Experiments on turbidity currents over an erodible bed, J. Hydraul. Res. 25, 123 (1987).

[21] T. van Kessel and C. Kranenburg, Gravity current of fluid mud on sloping bed, J. Hydraul. Eng. 122, 710 (1996).

[22] C. J. Dallimore, J. Imberger, and T. Ishikawa, Entrainment and turbulence in saline underflow in Lake Ogawara, J. Hydraul. Eng. 127, 937 (2001).

[23] R. L. Alpert, Turbulent ceiling-jet induced by large-scale fires, Combust. Sci. Technol. 11, 197 (1975).

[24] G. H. Jirka, Turbulent Buoyant Jets and Plumes, edited by W. Rodi, Vol. 6 (Pergamon Press, Oxford, 1982), pp. 69–119.
[25] M. Princevac, H. J. Fernando, and C. D. Whiteman, Turbulent entrainment into natural gravity-driven flows, *J. Fluid Mech.* **533**, 259 (2005).

[26] G. C. Christodoulou, Interfacial mixing in stratified flows, *J. Hydraul. Res.* **24**, 77 (1986).

[27] H. J. Fernando, Turbulent mixing in stratified fluids, *Annu. Rev. Fluid Mech.* **23**, 455 (1991).

[28] Q. Guo, Y. Z. Li, H. Ingason, Z. Yan, and H. Zhu, Theoretical studies on buoyancy-driven ceiling jets of tunnel fires with natural ventilation, *Fire Safety J.* **119**, 103228 (2021).

[29] D. Bonn, A. Andersen, and T. Bohr, Hydraulic jumps in a channel, *J. Fluid Mech.* **618**, 71 (2009).

[30] R. Dasgupta, G. Tomar, and R. Govindarajan, Numerical study of laminar, standing hydraulic jumps in a planar geometry, *Eur. Phys. J. E* **38**, 45 (2015).

[31] M. Dhar, G. Das, and P. K. Das, Planar hydraulic jumps in thin film flow, *J. Fluid Mech.* **884**, A11 (2020).

[32] V. T. Chow, *Open-Channel Hydraulics*, McGraw-Hill Civil Engineering Series (McGraw-Hill, New York, 1959).

[33] F. Henderson, *Open Channel Flow*, Macmillan Series in Civil Engineering (Macmillan, New York, 1966).

[34] L. Gu and G. A. Lawrence, Analytical solution for maximal frictional two-layer exchange flow, *J. Fluid Mech.* **543**, 1 (2005).

[35] J. P. Kunsch, Critical velocity and range of a fire-gas plume in a ventilated tunnel, *Atmos. Environ.* **33**, 13 (1998).

[36] L. Hu, R. Huo, Y. Li, H. Wang, and W. Chow, Full-scale burning tests on studying smoke temperature and velocity along a corridor, *Tunnell. Underground Space Tech.* **20**, 223 (2005).

[37] W. Chow, Y. Gao, J. Zhao, J. Dang, C. Chow, and L. Miao, Smoke movement in tilted tunnel fires with longitudinal ventilation, *Fire Safety J.* **75**, 14 (2015).

[38] W. Chow, Y. Gao, J. Zhao, J. Dang, and N. C. Chow, A study on tilted tunnel fire under natural ventilation, *Fire Safety J.* **81**, 44 (2016).

[39] D. L. Wilkinson and I. R. Wood, A rapidly varied flow phenomenon in a two-layer flow, *J. Fluid Mech.* **47**, 241 (1971).

[40] C.-S. Yih and C. R. Guha, Hydraulic jump in a fluid system of two layers, *Tellus* **7**, 358 (1955).

[41] A. Regev, S. Hassid, and M. Poreh, Density jumps in smoke flow along horizontal ceilings, *Fire Safety J.* **39**, 465 (2004).

[42] A. Regev, S. Hassid, and M. Poreh, Calculation of entrainment in density jumps, *Environ. Fluid Mech.* **6**, 407 (2006).

[43] N. Najafpour, M. Sarnie, B. Firoozabadi, and H. Afshin, Theoretical and experimental investigation of density jump on an inclined surface, *Sci. Iran.* **21**, 1655 (2014).

[44] J. R. Dormand and P. J. Prince, A reconsideration of some embedded Runge-Kutta formulae, *J. Comput. Appl. Math.* **15**, 203 (1986).

[45] F. Nicoud and F. Ducros, Subgrid-scale stress modelling based on the square of the velocity gradient tensor, *Flow Turbul. Combust.* **62**, 183 (1999).

[46] C. H. Bruneau and P. Fabrie, Effective downstream boundary conditions for incompressible Navier-Stokes equations, *Int. J. Numer. Methods Fluids* **19**, 693 (1994).

[47] C. H. Bruneau and P. Fabrie, New efficient boundary conditions for incompressible Navier-Stokes equations: A well-posedness result, *ESAIM: M2AN* **30**, 815 (1996).

[48] C. H. Bruneau, Boundary conditions on artificial frontiers for incompressible and compressible Navier-Stokes equations, *ESAIM: M2AN* **34**, 303 (2000).

[49] X. Zhou, K. H. Luo, and J. J. Williams, Large-eddy simulation of a turbulent forced plume, *Eur. J. Mech. B 20*, 233 (2001).

[50] R. Mehaddi, S. Vaux, F. Candelier, and O. Vauquelin, On the modelling of steady turbulent fountains, *Environ. Fluid Mech.* **15**, 1115 (2015).

[51] S. Vaux, R. Mehaddi, O. Vauquelin, and F. Candelier, Upward versus downward non-Boussinesq turbulent fountains, *J. Fluid Mech.* **867**, 374 (2019).

[52] W. H. Hager, R. Bremen, and N. Kawagoshi, Classical hydraulic jump: Length of roller, *J. Hydraul. Res.* **28**, 591 (1990).

[53] H. G. Hornung, C. Willert, and S. Turner, The flow field downstream of a hydraulic jump, *J. Fluid Mech.* **287**, 299 (1995).
[54] G. L. Richard and S. L. Gavrilyuk, The classical hydraulic jump in a model of shear shallow-water flows, *J. Fluid Mech.* **725**, 492 (2013).

[55] Z. Borden and E. Meiburg, Circulation based models for Boussinesq gravity currents, *Phys. Fluids* **25**, 101301 (2013).

[56] M. A. Khodkar, M. M. Nasr-Azadani, and E. Meiburg, Gravity currents propagating into two-layer stratified fluids: Vorticity-based models, *J. Fluid Mech.* **844**, 994 (2018).

[57] S. W. Bauer and W. H. Graf, Free overall as flow measuring device, *J. Irrig. Drain. Div.* **97**, 73 (1971).

[58] W. H. Graf and M. S. Altinakar, *Hydraulique Fluviale*, Vol. 16 (Presses Polytechniques et Universitaires Romandes, Lausanne, 2000).