On the modelling of shallow turbidity flows

TL;DR

This paper introduces a simple hyperbolic model for shallow turbidity flows that incorporates turbulence, sediment entrainment, and water mixing, validated against experimental data and capable of capturing various flow behaviors.

Contribution

A novel, minimal-parameter hyperbolic model for shallow turbidity flows that accounts for turbulence, sediment dynamics, and mixing, validated through experiments and analytical solutions.

Findings

Model accurately fits experimental data

Special solutions include travelling and self-similar flows

Unsteady simulations show long-term asymptotic states

Abstract

In this study we investigate shallow turbidity density currents and underflows from mechanical point of view. We propose a simple hyperbolic model for such flows. On one hand, our model is based on very basic conservation principles. On the other hand, the turbulent nature of the flow is also taken into account through the energy dissipation mechanism. Moreover, the mixing with the pure water along with sediments entrainment and deposition processes are considered, which makes the problem dynamically interesting. One of the main advantages of our model is that it requires the specification of only two modeling parameters - the rate of turbulent dissipation and the rate of the pure water entrainment. Consequently, the resulting model turns out to be very simple and self-consistent. This model is validated against several experimental data and several special classes of solutions (such as travelling, self-similar and steady) are constructed. Unsteady simulations show that some special solutions are realized as asymptotic long time states of dynamic trajectories.

Tables

Table 1. Table 1. Definitions of various parameters used in the approximate models. For the illustration see Figure  1 . A more complete list of employed nomenclature is given in  C .

Variable	Significance
$b_s$	Buoyancy of sediments
$b_{\mathrm{\hspace{0.17em}0}}$	Buoyancy of the lower (bottom) layer
$b(x,t)$	Buoyancy in the mixing layer
${\zeta }_s$	Thickness of the sediment deposit
$\zeta (x,t)$	Thickness of the lower (bottom) layer
$h(x,t)$	Thickness of the mixing layer
$w(x,t)$	Depth-averaged mean velocity in the bottom layer
$u(x,t)$	Depth-averaged mean velocity in the mixing layer
$q(x,t)$	Depth-averaged turbulent mean squared velocity in the mixing layer
$ℰ(x,t)$	Turbulent energy dissipation
$u_{\ast }$	Friction velocity

Table 2. Table 2. Parameters used in numerical simulations of the experiments from Rastello & Hopfinger (2004) [ 70 , Table 2] , schematically depicted in Figure  5 .

Parameter	Experiment 1	Experiment 2
Slope angle, $\phi$	${32}^{\circ }$	${45}^{\circ }$
Container height, $h_{\ell }$	$6.5\mathrm{𝖼𝗆}$	$20\mathrm{𝖼𝗆}$
Heavy fluid buoyancy, $b_{\ell }$	$19\frac{\mathrm{𝖼𝗆}}{s^2}$	$1\frac{\mathrm{𝖼𝗆}}{s^2}$
Sediment deposit height, ${\zeta }_s$	$0.2\mathrm{𝖼𝗆}$	$2.0\mathrm{𝖼𝗆}$
Minimal sediment “mass”, $m_s^{\circ }$	${10}^{-3}\frac{{\mathrm{𝖼𝗆}}^2}{{𝗌}^2}$	${10}^{-3}\frac{{\mathrm{𝖼𝗆}}^2}{{𝗌}^2}$
Sediment “mass”, $m_s$	$10\frac{{\mathrm{𝖼𝗆}}^2}{{𝗌}^2}$; $20\frac{{\mathrm{𝖼𝗆}}^2}{{𝗌}^2}$	$3\frac{{\mathrm{𝖼𝗆}}^2}{{𝗌}^2}$
Final simulation time, $T$	$10𝗌$	$20𝗌$

Equations

$${\zeta }_t+[\zeta w{]}_x=$$

$$h_t+[hu{]}_x=$$

$$w_t+w{w}_x+[b_0\zeta +bh{]}_x=$$

$$b_t+u{b}_x=$$

$$u_t+u{u}_x+b[\zeta +h{]}_x+\frac{1}{2}h{b}_x=$$

$$q_t+u{q}_x=$$

$$\displaystyle+\ \bigl{(}u^{\,2}\ -\ q^{\,2}\ -\ b\,h\bigr{)}\,\chi^{\,+}\ -\ \mathcal{E}\,\Bigr{]}\,.$$

$$\mathcal{E}\stackrel{\mathrm{def}}{:=}\kappa q^3,\text{mbox}where\kappa >0.$$

$${\chi }^{-}\equiv 0,{\chi }^{+}\equiv \sigma q.$$

$$\frac{u_{\ast }^2}{\zeta }\equiv -b_0{d}_x,$$

$$h_t+[hu{]}_x=$$

$$b_t+u{b}_x=$$

$$u_t+u{u}_x+b{h}_x+\frac{1}{2}h{b}_x=$$

$$q_t+u{q}_x=$$

$$\delta \stackrel{\mathrm{def}}{:=}\frac{\kappa }{\sigma }>0.$$

$$h_t+[hu{]}_x=$$

$$(hb{)}_t+[hbu{]}_x=$$

$$\displaystyle(h\,u)_{\,t}\ +\ \bigl{[}\,h\,u^{\,2}\ +\ {\textstyle{1\over 2}}\,b\,h^{\,2}\,\bigr{]}_{\,x}\ =$$

$$\displaystyle\bigl{(}h\,(u^{\,2}\ +\ q^{\,2}\ +h\,b)\bigr{)}_{\,t}\ +\ \bigl{[}\,(u^{\,2}\ +\ q^{\,2}\ +\ 2\,h\,b)\;h\,u\,\bigr{]}_{\,x}\ =$$

$$(h\,q)_{t}\ +\ [\,h\,q\,u\,]_{x}\ =\ {\textstyle{1\over 2}}\,\sigma\;\bigl{(}u^{2}\ +\ (1-\delta)\,q^{2}\ -\ h\,b\bigr{)}\,.$$

$$\displaystyle w\,\zeta_{\,x}\ +\ \zeta\,w_{\,x}\$$

$$\displaystyle u\,h_{\,x}\ +\ h\,u_{\,x}\$$

$$\displaystyle w\,w_{\,x}\ +\ b_{\,0}\,\zeta_{\,x}\ +\ b\,h_{\,x}\ +\ h\,b_{\,x}\$$

$$\displaystyle u\,u_{\,x}\ +\ b\,(\zeta_{\,x}\ +\ h_{\,x})\ +\ {\textstyle{1\over 2}}\,h\,b_{\,x}\$$

$$\displaystyle h\,u\,b_{\,x}\$$

$$\displaystyle 2\,q\,h\,u\,q_{\,x}\$$

$$\displaystyle+\ \bigl{(}u^{\,2}\ -\ q^{\,2}\ -\ b\,h\bigr{)}\,\chi^{\,+}\ -\ \kappa\,q^{\,3}\,.$$

$$\Delta\,(x)\ =\ \bigl{(}u^{\,2}\ -\ b\,h\bigr{)}\cdot\bigl{(}w^{\,2}\ -\ b_{\,0}\,\zeta\bigr{)}\ -\ b^{\,2}\,\zeta\,h\,.$$

$$\displaystyle h_{\,x}\$$

$$\displaystyle u_{\,x}\$$

$$\displaystyle b_{\,x}\$$

$$\displaystyle\zeta_{\,x}\$$

$$\displaystyle w_{\,x}\$$

$$\displaystyle q_{\,x}\$$

$$\displaystyle a_{\,1}\$$

$$\displaystyle a_{\,2}\$$

$${\chi }^{+}=2\sigma q,{\chi }^{-}=-\sigma q.$$

$$w_0^2>b_0{\zeta }_0,$$

$$b_{\,0}\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ \frac{\bigl{(}\rho_{\,0}\ -\ \rho_{\,a}\bigr{)}\,g}{\rho_{\,a}}\,.$$

$$b_L=\frac{1}{2}{b}_0,u_L=\frac{1}{2}{w}_0,q_L=\frac{w_0}{\sqrt{4+2\delta }}.$$

$$b_{\,x}\ =\ \frac{\sigma\,q\,\bigl{(}b_{\,0}\ -\ 2\,b\bigr{)}}{h\,u}\,.$$

$$\begin{gathered} \sigma =0.15,\delta =6,\phi =10.8^{\circ },x_0=8\mathsf{cm}, \\ {\zeta }_0=7\mathsf{cm},w_0=5\frac{\mathsf{cm}}{\mathsf{s}},b_0=1.4\frac{\mathsf{cm}}{{\mathsf{s}}^2}. \end{gathered}$$

Full text

Key words and phrases: turbidity currents; density flows; shallow water flows; conservation laws; finite volumes; travelling waves; self-similar solutions

MSC:

PACS:

∗ Corresponding author

Last modified:

1. 1 Introduction

2. 2 Mathematical model

3. 3 Analytical study of the model

4. 3.1 Mixing layer formation

5. 3.1.1 Problem statement and solution

6. 3.1.2 Sediment layer

7. 3.1.3 Intermediate conclusions

8. 3.2 Steady flows

9. 3.3 Travelling waves

10. 3.3.1 General considerations

11. 3.3.2 Stability of equilibria

12. 3.3.3 A particular class of travelling waves

13. 3.3.4 Further considerations

14. 3.4 Similarity solutions

15. 4 Model validation and unsteady simulations

16. 4.1 Problem formulation

17. 4.2 Experimental set-up

18. 4.3 Numerical results

19. 5 Discussion

20. 5.1 Conclusions

21. 5.2 Perspectives

22. Acknowledgments

23. A An exact solution

24. A.1 Steady flow

25. A.2 Travelling waves

26. B Derivation of the speed-Froude relation

27. C Nomenclature

Introduction

Underwater turbidity currents are sediment-laden underflows that play an important rôle in the morphology of the continental shelves (more generally of ocean bottoms) and in the global sediment cycle going to the formation of hydrocarbon reservoirs. We refer to [80] for a self-contained and comprehensive account of the theory of gravity currents and intrusions. The presence and entrainment of sediments differentiates them from stratified flows due to, *e.g.*temperature or salinity differences. The main physical mechanisms include the deposition, erosion and dispersion of important amounts of heavy sediment particles. Turbidity currents are not to be confused with debris flows, which represent fast-moving masses of poorly sorted heterogeneous material where interactions among the material pieces ($\approx$ particles) are important. Moreover, debris mix little with the ambient fluid. Debris flows have been a mainstream topic in the scientific literature due to their hazard they wreak in mountain regions (and not only).

The driving force is the gravity acceleration acting on dispersed sediment particles along steep and moderate bottom slopes. The initial perturbation is amplified by this acceleration, which in turn destabilizes the flow into shear instabilities that result in turbulent mixing and the transfer of mass and momentum. This gravity force creates the horizontal pressure gradient due to the increase of hydrostatic pressure resulting from the addition of particles. The heavy sediment particles are suspended in the mixing layer by fluid turbulence. The studied here processes are responsible of the transfer of littoral sediments to deep ocean regions. One should not disregard the destructive potential of gravity currents onto underwater structures such as pipelines, cables, *etc.*Turbidity currents in submarine canyons can attain surprisingly high velocities of the order of $8\ \thicksim\ 14\ {\sf m}/{\sf s}$ [43, 65]. These high velocities in the downstream direction result from the self-acceleration (and self-suspension) process from an appropriate initial perturbation, when more and more sediments are entrained by the flow from the bed, thus, increasing the rate of work performed by gravity [65]. This process is sometimes referred to as the “ignition” [64, 65, 24], which translates the energy imbalance property of such flows. One of important scientific questions is to determine the conditions necessarily to have an igniting flow. However, the self-acceleration stage cannot continue indefinitely. Most often the bed slope drops off (due to the bed morphology) or, simply, the sediment supply ceases. The mechanism of ignition was already described in Pantin (1979) [63]. However, the first laboratory demonstration of self-accelerated turbidity flows took $30$ more years [74].

Turbidity currents is a particular case of (continuously) stratified flows and they are fundamentally different from classical density underflows [17]. The main difference comes from the fact that the source of the density gradient, *i.e.*the suspended sediment, is not conservative. The suspended sediments are free to exchange with the core layer near the sea bed. The ambient still water is also entrained into this process. These exchanges are difficult to quantify and they constitute one of the main difficulties in the modeling of such flows [63, 64, 65]. In this respect turbidity currents are fundamentally non-conservative flows in their nature. Gravity flows may occur in the atmosphere111For instance, downslope windstorms over topography in Colorado (US) were observed and examined in [49, 60]. over topography, sub-aerial (*e.g.*avalanches, pyroclastic flows) and sub-aqueous environments (*e.g.*turbidity currents) over bathymetries. They may result also from anthropogenic activities such as when a dense buoyant industrial effluent or pollutant is released into a lake, river or ocean. In the present study we shall consider mainly sub-aqueous flows due to the abundance of available experimental data. We refer to [66, 55] as general excellent reviews on this topic.

Perhaps, the first serious attempts to observe turbidity currents in natural environments were performed in late 1960’s at Scripps Canyon offshore of La Jolla, California. They were reported in [38]. However, the flows reported in that study were so violent that the instrumentation was lost during these density currents making the detailed analysis extremely difficult [66]. The exact time moment of these underwater events is unpredictable which make them difficult to monitor in natural environments. Most of our physical knowledge on underwater turbidity currents come from small scale laboratory experiments [57, 24, 45, 44]. The experiments are bound to use common liquids for practical reasons. In general, it is not possible to respect all scalings. To give an example, we can mention the issue with particle sizes and their settling velocity. Nevertheless, taking into account the difficulties in obtaining field data, laboratory experiments are the only source of quantitative data about turbidity currents. The mathematical modelling is needed to extrapolate these experimental results to the scales on which these processes occur in nature. Nonetheless, the experiments offer a great opportunity for the verification of numerical results.

The gravity current can be divided geometrically into the flow head, body and tail. The head is shaped as an ellipse and, generally the head is higher than the flow body. In the present study we are mainly interested in the flow head modelling, where the most intensive mixing processes take place. Consequently, it influences the whole flow dynamics. The most advanced point of the flow head is called the front or nose.

The main difficulties in understanding the dynamics of gravity turbidity currents come from their genuinely turbulent nature. Moreover, the phenomenon is nonlinear, heterogeneous and unsteady. The flow complexity increases when the flow entrains more and more sediments in suspension. The literature devoted to the mathematical modeling of the density currents is abundant. First of all, we would like to mention the classical monographs on this subject [79, 77, 48]. The first and simplest models intended to explain the classical lock-exchange configurations was proposed in [36]. These models are referred to as integral, box or 0D models, since all quantities are averaged in space. The modern approaches to the mathematical modeling of such flows were initiated in [63, 64, 65]. A dense cloud 0D model for powder-snow avalanches including non-Boussinesq and sediment entrainment effects along the avalanche path was proposed in [70]. Powder-snow avalanches are large-scale, finite volume release turbidity currents (in the form of large scale suspension clouds) occurring on mountain slopes. These clouds sometimes reach $100\ {\sf m}$ in height and the front velocities of the order of $100\ \frac{{\sf m}}{{\sf s}}\,$. Without sediments (*i.e.*snow in the case of avalanches) distributed over the incline, the density current first accelerates and then decelerates without reaching important velocities. With sediments entrainment, the current can be maintained in the accelerating self-sustaining state during sufficient intervals of time to reach the velocities indicated above. In [34] a fair correlation of the avalanche velocity with the snow cover was demonstrated. The measurements of an avalanche front velocity in the Sion valley, Switzerland demonstrate a constant increase of the front velocity with traveled distance [12] (during the accelerating phase, of course). Thus, we come to the conclusion that the inclusion of sediments entrainment effect is of capital importance to predict the correct density current front velocity.

Some of recent studies devoted to the sediments transport within depth-averaged models include [8, 42, 22, 58, 5]. This list is far from being exhaustive. The shallow water approach assumes that vertical accelerations are negligible, so the pressure being essentially hydrostatic. The sediment concentration is a passive tracer with exchanges among different layers. The flow is fully turbulent, even if pure viscous effects are generally negligible. Moreover, the energy required to keep the sediments in the suspension cloud is a negligible portion of the total turbulent energy production [65]. Thus, the base model has to be first Reynolds-averaged [75] before applying the long wave approximation. Several authors made an effort to take into account the turbulence modeling into the shallow water type models [54, 19, 20]. Our approach to solve this issue will be detailed below. Nowadays, the multi-layer approaches to the density stratified flows become more and more popular [2]. Finally, some researchers chose a more CFD222Computational Fluid Dynamics (CFD).-like approach to the simulation of density flows incorporating eventually the advanced turbulence modeling [6, 18, 61, 62]. Perhaps, the first Direct Numerical Simulation (DNS) of the gravity current dates back to the years of $2\,000$ [32]. These simulations have an advantage of being depth-resolving and, thus, providing very a complete information about the flow structure in two or even three dimensions. However, due to the high computational complexity, only idealized academic configurations can be considered within reasonable CPU-time at the current state of technologies. Recently proposed three-dimensional (3D) turbidity-current models can be found in [37, 41, 35]. Moreover, the 3D DNS computations are often limited in the bulk Reynolds number.

In the present study we adopt a simplified (1.5D) approach along the lines of [73, 11, 47] based on the Eulerian formulation and depth-averaged formulations. A Lagrangian simplified BANG1D model was proposed in [69]. A simple 1.5D model was proposed in [39]. The authors parametrized their model by making the entrainment velocity depending on the dimensionless Richardson number. In the present study we close the model in an alternative way.

Very similar physical processes take place in powder-snow avalanches where the snow particles suspension flows down the mountains and the snow plays the rôle of sediments in underflows [34]. Consequently, very similar mathematical models appear in these two fields and to make the bibliography review more complete we mention some recent results in powder-snow avalanche modeling [59, 7, 1, 18, 78, 56, 13].

The goal here is to propose a simple multi-layer (just two or three layers) shallow water-type model which takes into account mixing processes between the layers. We assume that the sediment particles are well mixed across the height of each layer. So that their volume fraction can be effectively approximated by depth-averaged quantities. This model preferably has to be simple enough to be studied using even analytical methods. The main scientific question which is currently poorly understood is the influence of flow stratification on global flow patterns. For instance, the formation of current’s head (or the front region) and its steady-state velocity have to be carefully explained [47].

The present study is organized as follows. In Section 2 we state the physical and mathematical problem under consideration. In Section 3 we study the proposed model analytically and in Section 4 we validate its predictions by comparing them against experimental data. Some unsteady simulations and their relation to analytical (self-similar) solutions are presented in Section 4 as well. Finally, the main conclusions and perspectives of this study are outlined in Section 5.

Mathematical model

Consider an incompressible liquid which fills a two-dimensional fluid domain $\Omega$. A Cartesian coordinate system $O\,\boldsymbol{x}$, $\boldsymbol{x}\ =\ (x,\,y)$ is chosen in a classical way such that the axis $Oy$ points vertically upward and abscissa $O\,x$ is positive along the right horizontal direction. The fluid is bounded below by a solid non-erodible bottom $y\ =\ d\,(x)\,$. Above, the fluid can be assumed unbounded for the sake of simplicity, since our attention will be focused on processes taking place in the region close to the bottom.

The fluid is inhomogeneous and the flow can be conventionally divided in three parts. On the solid bottom there is a heavy fluid layer with density $\rho_{\,0}$ composed mainly of sedimentary deposits. Its thickness is $\zeta\,(x,\,t)\,$. Above, we have a muddle layer whose density will be denoted by $\bar{\rho}\,(\boldsymbol{x},\,t)$ composed of the sediments and the still water mix. Its thickness is $h\,(x,\,t)\,$. Finally, the whole domain above these two layers is filled with the still ambient water of the density $\rho_{\,a}\ >\ 0$ at rest ($\boldsymbol{u}_{\,a}\ \equiv\ \boldsymbol{0}$). This situation is schematically depicted in Figure 1. We include also into consideration the situation where a thin motionless layer of sediments with constant thickness $\zeta_{\,s}$ and density $\rho_{\,s}\ >\ 0$ covers the slope. It is usually the case in many practical situations and these sediment deposits may contribute to the flow head dynamics while propagating along the incline. One may imagine also that a certain mass of a sediment suspension or dense fluid is released into the flow domain at $x\ =\ 0$ with the mass flux $\rho_{\,0}\cdot\zeta\,(0,\,t)\cdot u\,(0,\,t)\,$.

After applying the standard long wave scaling (or, equivalently, the depth-averaging), one can derive the following system of equations (c.f.[47]):

[TABLE]

The sense of the variables is explained in Table 1 and also in Figure 1. For the dissipative term, we assume the following closure relation:

[TABLE]

In the sequel we are going to make the following choice of the entrainment rates $\chi^{\,\pm}$ (unless the contrary is stated explicitly):

[TABLE]

Moreover, we assume that the gravity force balances exactly the turbulent friction at the solid bottom, i.e.

[TABLE]

which implies that the kinematic Equation (2.1) is completely decoupled from the rest of the system (2.2) – (2.6). The information about $\zeta$ is transported along characteristics of this equation. So, if $\zeta$ is constant initially and this constant value is maintained at the channel inflow, $\zeta$ will remain so under the system dynamics. We will adopt this assumption as well in order to focus our attention on the mixing layer dynamics. Below we will consider only the middle mixing layer. Under these conditions, the equilibrium model becomes:

[TABLE]

where the constant $\delta$ is defined as:

[TABLE]

The system (2.7) – (2.10) can be equivalently recast in the conservative form which has an advantage to be valid for discontinuous solutions as well [72, 27]:

[TABLE]

Here, $\kappa\ \in\ \mathds{R}^{\,+}$ is a positive constant measuring the rate of turbulent dissipation [75, 46].

An important property of this model is the absence of the bottom friction term. This physical effect was taken into account, while excluding the bottom layer of thickness $\zeta(x,\,t)\,$. This gives us the mathematical reason for the absence of a friction term in model (2.11) – (2.14). Similarly, we can bring also a physical argument to support this fact. The horizontal velocity takes the maximum value on the boundary between the bottom sediment and mixing layers. Consequently, the Reynolds stress $\tau_{\ast}$ vanishes here.

Remark 1**.**

In fact, we can derive an additional balance law by combining together Equations (2.7) and (2.10):

[TABLE]

The last equation is not independent from Equations (2.11) – (2.14). Consequently, it does not bring new information about the equilibrium sedimentation-free model (2.7) – (2.10). Nevertheless, we provide it here for the sake of the exposition completeness.

The proposed model has an advantage of being simple, almost physically self-consistent333Only two constants $\sigma$ and $\delta$ need to be prescribed to close the system. and having the hyperbolic structure. It was derived for the first time in [48, Chapter 5]. In order to obtain a well-posed problem the system (2.11) – (2.14) has to be completed by corresponding boundary and initial conditions. In the following sections, the just proposed equilibrium system (2.11) – (2.14) will be studied in more details by analytical and numerical means.

Analytical study of the model

In this Section we discuss the steady state solutions and investigate the qualitative behaviour of two important special classes of solutions — travelling waves and similarity solutions.

Mixing layer formation

Mixing layers are formed between two fluid layers having different densities. The initial stages of mixing processes play a very important rôle in natural and laboratory environments. The physical mechanism of a mixing layer formation is given by various interfacial instabilities (*e.g.*Rayleigh–Taylor or Kelvin–Helmholtz). The most important mechanism is, of course, the Kelvin–Helmholtz instability [33] because of the shear velocity presence [67]. In supercritical flows the mixing intensity is considerably intensified. In this way, mixing layers are the most pronounced in transcritical flows over underwater obstacles (*i.e.*bathymetric features), since a jet is formed on the obstacle downstream (lee) side making the flow supercritical. The process of mixing layer formation over an inclined bottom was studied in [47] in the framework of a three-layer model. In the present study we apply a similar approach to the two-layer444We remind that the upper layer of the still water is assumed to be motionless in the present study. System (2.1) – (2.6). Below we derive a steady solution, which provides the boundary conditions for the unsteady propagation of the flow head (see Figure 1).

Stationary solutions to System (2.1) – (2.6) satisfy the following system of differential equations:

[TABLE]

The last system of equations can be seen as a quasilinear system with respect to spatial derivatives of unknown variables \bigl{(}h_{\,x},\,u_{\,x},\,b_{\,x},\,\zeta_{\,x},\,w_{\,x},\,q_{\,x}\bigr{)}\,. The determinant $\Delta\,(x)$ of this system can be easily computed:

[TABLE]

By assuming that $\Delta\,(x)\ \neq\ 0\,$, we obtain:

[TABLE]

where the coefficients $a_{\,1,\,2}$ are defined as

[TABLE]

The last system of equations can be seen as an Initial Value Problem (IVP) in the spatial variable $x\,$, if all spatial derivatives in the right-hand side are replaced by their expressions (we do not make this operation to keep the shorthand notation). We illustrate below the behaviour of solutions to System (3.2) – (3.7) on the example of the mixing layer formation problem over an inclined bottom. The real data to build this solution were taken from [68].

3.1.1  Problem statement and solution

In this Section we formulate the IVP inspired by the experimental study [68]. Consider a density current over the flat plane inclined with angle $\varphi$ with respect to the horizontal direction. Without any loss of generality, we postulate that the mixing layer starts to form at some location $x\ =\ x_{\,0}\,$, where we set the initial conditions for the steady System (3.2) – (3.7). Moreover, we adopt the following closure:

[TABLE]

We assume that in this point we have a supercritical flow with velocity $w_{\,0}\,$, height $\zeta_{\,0}$ and density $\rho_{\,0}\,$. The upper layer of still water with density $\rho_{\,a}$ is motionless. Thus, we have

[TABLE]

and $b_{\,0}$ is defined as

[TABLE]

The initial height of the mixing layer is $h_{\,0}\ =\ 0$ by our assumption on point $x_{\,0}\,$. The initial asymptotic stages of the mixing layer formation as $x\ \to\ x_{\,0}$ can be determined from the condition that the right hand sides in (3.2) – (3.7) are bounded [48, 47]:

[TABLE]

Under the condition $\zeta\ >\ 0\,$, we shall have $b\,(x)\ \equiv\ b_{\,L}\ =\ {\textstyle{1\over 2}}\,b_{\,0}\,$. This can be seen after substituting the closure relations for $\chi^{\,\pm}$ from Equation (3.8) into Equation (3.4). One readily obtains that

[TABLE]

From the condition $b_{\,L}\ =\ {\textstyle{1\over 2}}\,b_{\,0}\,$, it follows that $b\,(x)\ \equiv\ {\textstyle{1\over 2}}\,b_{\,0}\,$.

By taking into account the last asymptotics for the solution, we can construct it using standard numerical means until the point $x_{\,1}\,$, where the mixing layer will touch the bottom, *i.e.*the point where $\zeta(x_{\,1})\ \equiv\ 0\,$. The solution to system (3.2) – (3.7) taking into account the aforementioned asymptotics is depicted in Figure 2. We used the following set of physical parameters in our computations:

[TABLE]

The flow geometry along with the initial conditions were taken from [68]. On the upper panel of Figure 2 we show the comparison with experimental data provided by Armi & Pawlak (2000) [68], who employed the Laser Induced Fluorescence (LIF) and Digital Particle Imaging Velocimetry (DPIV) together with a concentration of Rhodamine dye for flow visualization. Color code corresponds to the density gradient (blue is lower, red is higher). Their results show that instabilities evolve in an asymmetrical manner. In Figure 2 with a dashed line we show the numerical prediction given by function $y\ =\ \zeta\,(x)$ (lower line) and $y\ =\ \zeta\,(x)\ +\ h\,(x)$ (upper line). The reasonable agreement between the experimental data and our numerical solution at the initial parts of the mixing layer validate the approximate model (at least for steady solutions). We would like to mention that in the experiment (as well as in the nature), there is a slight backward flow above the mixing layer. The typical velocities are of the order $0.5\ \thicksim\ 1\;\mathsf{\frac{{\sf cm}}{{\sf s}}}\,$. In our model this effect is not taken into account. It could be done by including the third (upper) layer into consideration [47]. On the left and right sub-plots in upper Figure 2, we show the experimental distribution of the velocity (solid line) and of the density (dashed line) at the beginning and at the end of the considered fluid domain.

The theoretical lower bound of the mixing layer coincides fairly well with the experimental estimation of the high density gradient area as it can be seen in Figure 2 (upper panel). On the other hand, the theoretical upper bound goes sufficiently higher than the coloured region measured experimentally. This little discrepancy comes from the fact that the experimental data report the high density gradient, while our model predicts the upper bound of large vortices appearing during nonlinear stages of the Kelvin–Helmholtz (KH) instability development [33]. This situation is described in more details in [48, Chapter 7, § 3]. The experimental evidence for the KH mechanism is shown in [68, Figure 8]. In particular, it is shown there that the mixing layer growth into the still water is

[TABLE]

In the same time, the effective thickness of the mixing layer grows according the computed velocity profile as

[TABLE]

This value corresponds much better to numerous experimental findings [48, Fig. 7.4]. Consequently, we may conclude that in the considered case of the stratified fluid flow down the incline, the effective height555We remind here the definition of the effective height of the mixing layer: it is the fluid body bounded from above and below by virtual surfaces where $w\ =\ 0.95\,\times\,w_{\,0}$ and $w\ =\ 0.1\,\times\,w_{\,0}$ correspondingly. The asymmetry in this definition comes from the nonlinear form of the vertical velocity profile. of the mixing layer will be smaller than the height of the computed turbulent layer.

In Figure 2 (lower panel) we show with dotted lines the (theoretically predicted) distribution of velocities $u\,(x)$ and $w\,(x)$ in the mixing layer. The slope rupture in these curves is visible. It happens since at $x_{\,c}\ \approx\ 26\;\textsf{cm}$ the determinant (3.1) vanishes, i.e.

[TABLE]

Physically speaking it means that the flow remains everywhere supercritical. However, for small changes of the flow parameters, it is possible that hydraulic jumps will appear along with following local subcritical zones. We note also that internal hydraulic jumps on the downstream side of an underwater obstacle is an important feature of stratified flows (both in the ocean and in the atmosphere). The inclusion of the mixing layer formation before the hydraulic jump is of capital importance in such situations.

3.1.2  Sediment layer

When the mixing layer approaches the lower boundary with sediment deposits, the mass entrainment from the lowest layer decelerates and stops completely. In the framework of our model, this effect could be realized by ensuring the transition from $\chi^{\,-}\ =\ -\,\sigma\,q$ towards $\chi^{\,-}\ =\ 0$ (and, correspondingly, from $\chi^{\,+}\ =\ 2\,\sigma\,q$ towards $\chi^{\,+}\ =\ \sigma\,q$). Physically it means that the maximum of the flow velocity is achieved somewhere near the boundary between the mixing layer and bottom layer of constant density $\rho_{\,0}\,$. The transition to $\chi^{\,-}\ =\ 0$ takes place in neighbourhood of the right boundary of the flow represented in Figure 2.

3.1.3  Intermediate conclusions

The mixing layer structure over a downhill determines the flow structure further down. In particular, the total mass flux (relative to fluid density $\rho_{\,a}$) $\mathcal{M}_{\,0}\ =\ b_{\,0}\,\zeta_{\,0}\,w_{\,0}$ is divided into two parts. Namely, this splitting takes place at the point of transition of the mixing layer into the turbulent jet and the undiluted bottom layer, *i.e.*at $x\ =\ x_{\,1}\,$.

The turbulent jet receives the buoyancy flux $\mathcal{M}_{\,j}\ =\ b_{\,1}\,h_{\,1}\,u_{\,1}\,$, while the bottom layer takes $\mathcal{M}_{\,b}\ =\ b_{\,0}\,\zeta_{\,1}\,w_{\,1}\,$. We notice also that for the slope angle $\varphi\ =\ 10.8^{\,\circ}$ and experimental facilities considered in Figure 2, the dominant part of the mass flow goes into the turbulent jet, i.e.$\mathcal{M}_{\,j}\ \gg\ \mathcal{M}_{\,b}\,$.

If we assume that the flow velocity at $x\ =\ x_{\,1}$ in the bottom layer stabilizes666Physically it happens when the gravity force projection along the slope is balanced by the friction force with the rigid bottom. due to small thickness of the layer, in the bottom layer we have $w\ \equiv\ w_{\,b}\,$, $\zeta\ \equiv\ \zeta_{\,b}$ and the density in the layer does not change anymore. Hence, the bottom layer becomes completely passive in the flow dynamics. But in the same time it influences the flow head propagation velocity for moderate bottom slopes777For large inclines the physical mechanisms are slightly different.. This point will be demonstrated below.

Steady flows

From numerical and experimental points of view, a significant number of studies has been devoted to the study of the transient gravity current problem with special focus on the density front formation and evolution [9, 34, 53]. The steady flow configuration has received less attention [25]. However, when one has a dynamical system in hands, it is natural to begin its study by looking for equilibrium points [40]. Thus, the study of a mathematical model is not complete if we do not discuss this class of solutions. Moreover, stationary solutions may be realized and observed on certain time scales in laboratory experiments, if certain conditions are met and maintained during sufficiently long time. For instance, the buoyant flow entering the inclined channel has to be maintained at constant rate during the whole experiment. If permanent boundary conditions are maintained not only in the mixing layer, but also in the bottom layer as well, then, the lower sediment layer will play a passive rôle in the flow, i.e.

[TABLE]

In other words, the friction and mixing between these two layers do not take place. It is not difficult to see that fluid flows of the form

[TABLE]

satisfy Equations (2.7) – (2.10) provided that the mixing layer thickness depends linearly on the coordinate $x$ along the channel, i.e.

[TABLE]

Here, quantities with the sub-script j denote values in the steady jet. The variable $m$ will be sometimes referred as the mass, however, physically it represents the excess of the fluid column weight with respect to the still water level due to the presence of heavy sediment suspensions. To satisfy the system (2.7) – (2.10), the following identities have to be satisfied:

[TABLE]

where $\upalpha\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ -\,d_{\,x}\ \geqslant\ 0\,$. An algebraic consequence of relations (3.9) – (3.11) above can be easily derived:

[TABLE]

The boundary conditions specify also the buoyancy influx $\mathcal{M}_{\,j}$ into the channel:

[TABLE]

All the relations above can be combined into a single equation for the quantity $\mathcal{A}_{\,j}\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ \dfrac{u_{\,j}}{U_{\,j}}$ using simple algebraic transformations (see also A.1):

[TABLE]

There exists a unique solution to this equation for any positive right hand side. Moreover, one can show that $\mathcal{A}_{\,j}\ >\ 1\,$. It corresponds to the supercritical flow with $u_{\,j}^{\,2}\ >\ m_{\,j}\,$. Then, once we determined $\mathcal{A}_{\,j}\,$, from Equations (3.9) – (3.11) we can determine the remaining quantities $u_{\,j}\,$, $m_{\,j}\,$, $q_{\,j}$ and $\varsigma\,$.

To make a conclusion, in turbidity gravity flows in an inclined channel, where the horizontal (depth-integrated) velocity maximum is achieved on the boundary between the sediment and turbulent mixing layers, the model (2.7) – (2.10) predicts the stationary supercritical flow with pure fluid entrainment from the upper layer.

Travelling waves

One of the main questions in the modelling of turbidity flows is to determine the density front velocity. The proposed base model (2.7) – (2.10) is sufficiently simple to address this question analytically.

3.3.1  General considerations

In this Section we describe the class of travelling wave solutions to System (2.11) – (2.14) in their generality. Since we are mainly interested in smooth solutions, we will use the non-conservative form (2.7) – (2.10) for the sake of convenience. Also, we assume that the bottom slope $\upalpha$ is constant, which is necessary for the existence of solutions with permanent shape.

The travelling wave ansatz takes the following form:

[TABLE]

where $\xi\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ x\ -\ c\cdot t$ and $c$ is a positive888If the wave travels to the rightwards direction, as we assume without any loss of generality. constant, which has the physical sense of the travelling wave speed (to be determined later). After substituting this ansatz into Equations (2.7) – (2.10), we obtain a system of Ordinary Differential Equations (ODEs):

[TABLE]

where the prime ′ denotes the differentiation operation with respect to $\xi\,$. This system of ODEs describes a transcritical gravity flow in a coordinate system which moves with velocity $c\,$. More precisely, the flow is of supercritical type in the avalanche core and it switches to the subcritical regime when we cross the wave front. By analogy with the detonation theory [23], the front velocity $c$ (under certain conditions) is given by the Chapman–Jouguet principle which states that a transcritical front propagates with the minimal admissible velocity [10].

After some computations the ODEs System (3.12) – (3.15) can be reduced to a single differential equation with implicit dependence on $\xi\,$:

[TABLE]

where the expression for $m\,(\xi)$ in terms of $u\,(\xi)$ is obtained by integrating Equation (3.13):

[TABLE]

The integration “constant” $\mathcal{M}\,(c)$ can be determined from the Rankine–Hugoniot conditions written at the wave front [72, 27].

3.3.2  Stability of equilibria

In this Section we study the existence and stability of equilibria points to Equation 3.16. One of the difficulties is that the right hand side depends on a free parameter $\delta\,$, whose variation has to be taken into account. The wave celerity $c\ >\ u$ has to be specified as well.

As the first step, we rewrite Equation (3.16) as a dynamical system in the plane \bigl{(}q,\,u\bigr{)}\, with $\xi$ being the evolution variable:

[TABLE]

The equilibria points might be of two kinds:

1. (1)

on the intersection of curves \Gamma_{\,\mathscr{A}}\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ \bigl{\{}\mathscr{A}\,(q,\,u)\ =\ 0\bigr{\}} and \Gamma_{\,\mathscr{B}}\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ \bigl{\{}\mathscr{B}\,(q,\,u)\ =\ 0\bigr{\}}\,, 2. (2)

on the intersection of \bigl{\{}\mathscr{D}\,(u)\ =\ 0\bigr{\}} and \bigl{\{}\mathscr{B}\,(q,\,u)\ =\ 0\bigr{\}}\,.

The linear stability of equilibria is determined by eigenvalues of the following Jacobian matrix:

[TABLE]

The elements of the Jacobian matrix can be computed by direct differentiation:

[TABLE]

The last expression can be simplified at equilibria locations by taking into account the fact that $\mathscr{A}\ =\ \mathscr{B}\ \equiv\ 0$ at the first kind and $\mathscr{D}\ =\ \mathscr{B}\ \equiv\ 0$ at equilibria states of the second kind.

Eigenvalues $\lambda_{1,\,2}\,(\delta)$ of the Jacobian matrix $\mathscr{J}$ as functions of the parameter $\delta$ for both types of equilibria points are shown in Figure 3. The physical parameters used in this numerical computation are $\sigma\ =\ 0.15\,$, $\upalpha\ =\ \tan\varphi\ =\ 1.0\,$, $c\ =\ 2.787$ and $\mathcal{M}\ =\ 0.1\,c\,$. From this illustration it is clear now that, at least for these values of parameters, the equilibria of the first kind are unstable spiral points, while the equilibria of the second kind are (linearly) stable spirals [51].

3.3.3  A particular class of travelling waves

Similarly to the construction of steady solutions presented in the previous Section 3.2, we are looking for travelling wave solutions to system (2.7) – (2.10) of the form:

[TABLE]

where $\xi\ \equiv\ x\ -\ c\cdot t$ is a combined independent variable introduced earlier. In other words, we consider the frame of reference where the travelling wave is steady. Notice that the model (2.7) – (2.10) does not possess the Galilean invariance property because of the velocity variables present in the right hand sides. This effect comes from the assumption that the ambient fluid remains in rest, which privileges this particular frame of reference. The constant $c\ >\ 0$ is the unknown wave celerity to be determined during the solution procedure.

Consider travelling waves of the following particular form:

[TABLE]

with $\varsigma\ >\ 0$ and $\xi\ <\ 0\,$. The travelling wave ansatz presented above satisfies the following relations:

[TABLE]

We use also an additional hypothesis that the flow in the coordinate frame moving with the travelling wave is critical, i.e.

[TABLE]

The last condition can be equivalently recast as

[TABLE]

where $\mathrm{Fr}_{\,f}$ is the Froude number with respect to the wave front [21]. This condition ensures the existence of the self-sustained regime of the wave propagation independently of small perturbations which might occur behind the wave front. It is completely analogous to the so-called Chapman–Jouguet condition for the propagation of a self-sustained detonation wave in gas dynamics [82].

The wave celerity $c$ might be eliminated from Equations (3.19) – (3.22) by introducing new variables:

[TABLE]

As a result, we come to the following closed system of equations:

[TABLE]

By assuming that $\sigma\cdot q^{\,\star}\ \neq\ 0\,$, from last equations we can derive a simple relation:

[TABLE]

and by using the scaled version (3.26) of relation (3.22), we obtain that

[TABLE]

or formally,

[TABLE]

provided that $u^{\,\star}\ >\ {\textstyle{1\over 2}}$ so that the value $q^{\,\star}\ \in\ \mathds{R}^{\,+}\,$. After dividing (3.24) by (3.23) we have:

[TABLE]

In other words, we have the following equation for $u^{\,\star}\,$:

[TABLE]

Under the same condition on $u^{\,\star}\,$, we can transform Equation (3.27) into the following algebraic equation:

[TABLE]

It can be shown (see A.2) that there exists a unique positive root to this equation, which belongs to the interval $u^{\,\star}\ \in\ (\,{\textstyle{1\over 2}},\,1\,)\,$. All other remaining quantities such as $\varsigma\,$, $m^{\,\star}$ and $q^{\,\star}$ are determined after finding $u^{\,\star}\,$. The velocity $c$ can be found from the mass conservation equation written at the wave front:

[TABLE]

If $m_{\,s}\ \equiv\ 0\,$, the last equation takes a particularly simple form:

[TABLE]

where we introduced some notations:

[TABLE]

We remind that the positive constant $c_{\,w}$ controls the magnitude of bottom friction effects. We thus arrive to the following cubic polynomial equation for the quantity $C\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ \dfrac{c}{U_{\,b}}\,$:

[TABLE]

It is not difficult to see that the polynomial $\mathcal{P}(C)$ is a monotonically increasing function of $C$ provided that $u^{\,\star}\ \leqslant\ 1$ and, thus, there exists a unique positive root (since $\mathcal{P}\,(0)\ =\ -1\ <\ 0$). The dependence of the computed in this way coefficient $C\ =\ C_{\,e}\,(\mathrm{Fr}_{\,b})\ =\ C_{\,e}\,(\varphi)$ on the channel slope angle $\varphi$ is shown in Figure 4 with the dashed line (for a fixed value of $c_{\,w}\ =\ 0.004$ and with $\sigma\ =\ 0.15\,$, $\delta\ =\ 4\,$). Note that in the experimental study [9] the notation $C\ \equiv\ \dfrac{c}{U_{\,0}}$ was used, while we consider another definition $C\ \equiv\ \dfrac{c}{U_{\,b}}\,$, where $U_{\,0}\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ \sqrt[3]{\mathcal{M}_{\,0}}\,$. We can see that this prediction does not compare very well with the experimental data [26, 81, 76, 9] for large slope angles $\varphi\,$. This drawback will be corrected below. The agreement for small angles of inclination is achieved since $\mathcal{M}_{\,j}\ \ll\ \mathcal{M}_{\,b}$ and, consequently, $U_{\,0}\ -\ U_{\,b}\ \ll\ U_{\,b}\,$. Therefore, the total mass flux $\mathcal{M}_{\,0}$ can be replaced by the bottom flux $\mathcal{M}_{\boldsymbol{b}}$ to find the density current head velocity. It is also worth to mention that the big scatter in experimental data depicted in Figure 4 can be also partly explained by the differences between the total $\mathcal{M}_{\,0}$ and bottom $\mathcal{M}_{\,b}$ mass fluxes, which depend on precise inflow conditions.

3.3.4  Further considerations

The stability of the constructed travelling wave solution with the linear growth of the profile (due to the perpetual entrainment of sediments and still water into the mixing layer) has to be studied separately. The constructed solution belongs to the important class of stratified flows with fluid particles (in some parts of this flow) moving faster than the wave front. It is known that for surface waves it leads to inevitable wave breaking [52]. The authors are not aware of any mathematical stability studies of such stratified flows.

In the modelling of density currents, where there are heavy particles (of density $\rho_{\,0}$) entrained into the flow, the “boundary” conditions imposed on the wave front have the capital importance. In the previous Section these conditions were determined by the flow structure and we saw that it leads to the poor prediction of the wave celerity comparing with experimental data [26, 81, 76, 9]. Here we propose another set of “boundary” conditions:

[TABLE]

which are composed of the mass conservation equation, the flow criticality condition relative to the wave front and also the momentum conservation equation. It is not difficult to see that for $u\ <\ c\,$, it follows that $u\ =\ {\textstyle{1\over 3}}\;c\,$ and for the dimensionless combination $C\ =\ \frac{c}{U_{\,b}}$ we obtain the following cubic equation:

[TABLE]

All intermediate computations are explained in B. For large values of the Froude number, the solution $C\ \approx\ 1.5\,$. This estimation is in good agreement with the empirical conjecture on the density current flow head velocity based on the experiments reported in [9]:

[TABLE]

In Figure 4 the solid line shows the dependence of the flow head celerity $C\ =\ C_{\,b}\,(\varphi)$ based on Equation (3.31). The overall good agreement of the lower solid curve with the experimental data [26, 81, 76, 9] certifies the model quality. It is worth to notice that the proportionality coefficient in experimental studies was determined based on velocity $\mathcal{M}_{\,0}^{\,\frac{1}{3}}$ and not on velocity $U_{\,b}\,$. As a result, for more accurate comparisons with considered models, one has to determine accurately, which part of the mass flow $\mathcal{M}_{\,0}$ entered into the boundary layer, since U_{\,b}\ =\ \mathcal{M}_{\,b}^{\,\frac{1}{3}}\ =\ \Bigl{(}\mu\,\mathcal{M}_{\,0}\Bigr{)}^{\,\frac{1}{3}}\,, with $\mu\ \leqslant\ 1\,$.

Similarity solutions

Additionally to special solutions considered in two previous Sections, System (2.7) – (2.10) admits also the following class self-similar solutions:

[TABLE]

where $\xi\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ \dfrac{x}{t^{\,\uptheta+1}}$ and $\uptheta\ \in\ \mathds{R}\,$. For flows with the sustained mass flow $\mathcal{M}_{\,b}\ \equiv\ \mathrm{const}$ or with the sustained sediment mass $m_{\,s}\ \equiv\ \mathrm{const}\,$, from conditions on the front it follows that $\uptheta\ \equiv\ 0\,$. For the sustained mass flow in the turbulent mixing layer $\mathcal{M}_{\,j}\ \equiv\ \mathrm{const}\,$, the same value for $\uptheta$ follows from boundary conditions, i.e.$\uptheta\ \equiv\ 0\,$. Henceforth, for all just mentioned cases, we look for solutions of the form:

[TABLE]

with $\xi\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ \dfrac{x}{t}\,$. Such solutions play a very important rôle in understanding the system (2.7) – (2.10), since dynamic solutions tend to self-similar ones at large times provided that constant mass fluxes are maintained.

Another important class of self-similar solutions (3.32) is realized when $\mathcal{M}_{\,0}\ \equiv\ 0\,$, $m_{\,s}\ \equiv\ 0\,$. This corresponds to the evolution of a finite mass of a heavier liquid, which propagates under the layer of a lighter fluid. In this case, the mass conservation law

[TABLE]

yields self-similar solutions of the form (3.32):

[TABLE]

where $\xi\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ \dfrac{x}{t^{\,2/3}}\,$. This self-similar solution indicates [71] that the wave front position $x_{\,f}$ behaves asymptotically in time as:

[TABLE]

From the last estimation, the following asymptotic behaviour of the front velocity $U_{\,f}$ can be readily deduced:

[TABLE]

The last decay law was validated experimentally in [4, 53]. In the present work this asymptotic behaviour will be used to validate the numerical simulations. This approach to the description of turbidity currents is referred as the “thermal theory” [4].

In turbidity flows, in order to construct self-similar solutions for density currents along a slope (as well as for travelling waves), it is of capital importance to prescribe the adequate boundary conditions. During the construction of travelling waves above, we see that the wave celerity depends on the relations imposed at the wave front. When the solution is given by ansatz (3.32), the flow head is uniquely determined. However, the applicability of such similarity solutions to flows in the absence of sediments (i.e.$m_{\,s}\ \equiv\ 0\,$, $\mathcal{M}_{\,b}\ \equiv\ 0$) has to be studied separately.

The situation changes when we consider the problem of sediments entrainment into the flow ($m_{\,s}\ >\ 0\,$, $\mathcal{M}_{\,b}\ \equiv\ 0$). In the framework of the conservative system (2.11) – (2.14) this problem can be interpreted as the mixed Initial–Boundary Value Problem (IBVP). Namely, at the initial moment of time we know the distribution of heavy sediments in the sediment layer $h_{\,s}\,(x)\,$, $m_{\,s}\,(x)\,$. The sediments might be at rest ($u_{\,s}\,(x)\ \equiv\ 0$) or moving with prescribed velocity $u_{\,s}\,(x)$ and turbulent kinetic energy $q_{\,s}\,(x)\,$. The perturbations entraining the sediments into the flow are entering the fluid domain from the left boundary (without any loss of generality). In Figure 1 we depicted the sketch of the fluid domain for an illustration. In this situation there is no need to separate the flow head and to write additional relations on the wave front, since the main flow characteristics are obtained in the process of solving the IBVP. However, the main question remains unanswered: which self-similar regime (3.32) will appear as the long time limit of the unsteady solution?

System (2.11) – (2.14) is of hyperbolic (and hydrodynamic) type and it corresponds to the flow of a barotropic gas with chemical reactions in compressible fluid dynamics. We have already mentioned above the analogy between density currents and the detonation theory. During the normal detonation, the fluid flow satisfying the Chapman–Jouguet conditions downstream the flow head is not always realized (especially in the presence of accompanying chemical reactions). Under certain conditions the wave of detonation might propagate with the velocity exceeding that of perturbations behind the wave front [23]. The main conclusion that we can draw from this analogy is that the realizability of self-similar solutions has to be studied separately. Such an analytical study might turn out to be very complex. Below we will show by numerical means that the propagation of the flow head with the speed higher than expected is possible. We would like to mention also that similarity solutions for gravity currents were constructed also in [71] and [39, Sections 3 & 4].

Model validation and unsteady simulations

Strictly speaking, we already validated steady solutions to System (3.2) – (3.7) by making comparisons with the experiments from [68]. We show that this system is able to predict qualitatively and quantitatively the development of the mixing layer over a slope.

Problem formulation

In this Section we continue the validations by considering unsteady solutions hereinafter. Moreover, we shall consider the applicability of an even simpler one-layer model:

[TABLE]

We apply it to simulate the sediments entrainment process by density currents over moderate (finite) slopes. Numerical solutions will be compared with experimental data reported in [70] as well as with exact special solutions of travelling wave type. In their experiments the dense fluid consisted either of saltwater or of the sawdust particle suspension. The shape of these particles was irregular in accordance with suspended snow flakes. The spatial growth of the cloud was determined from the side view images recorded with a video camera. A $5\ {\sf cm}$ square grid was drawn on the side glass to facilitate the front position determination. Moreover, we shall show below that under certain initial conditions the numerical solution will tend asymptotically to self-similar solutions described earlier.

We made our choice for experimental data of Rastello & Hopfinger (2004) [70] for the following reasons:

• •

This experiment corresponds fairly well to the scope and purpose of our numerical model

• •

The model (4.1) – (4.4) is suitable for the simulation of gravity currents even over moderate and large slopes used in the experimental study [70] (see Table 2 for the values of the slope parameter). When using other models, one has to check that the influence of the sediment bottom layer is taken into account to represent correctly the front dynamics

• •

The experiment was conducted for sufficiently long time. In this way we are able to check our model validity for the acceleration and deceleration stages of the flow. Finally, we were even able to check the asymptotic behaviour of the flow, which is very typical for buoyant flows over a slope.

We already mentioned above that the equilibrium system (4.1) – (4.4) possesses the hyperbolic structure similar to the equations of gas dynamics with two “sonic” and two contact characteristics [72, 27]. There is a wide choice for the numerical discretization of such systems. However, taking into account that we deal with a system of conservation (balance) laws, it seems to be natural to opt for finite volume schemes [3]. Henceforth, in order to solve numerically the equilibrium system (4.1) – (4.4) we use the classical and robust finite volume discretization and the widely used Godunov scheme [28, 30, 29]. The explicit Euler scheme is used in time discretization.

Sketch of the numerical experiment is given in Figure 5. Our goal is to reproduce in silico some of the laboratory experiments reported in [70]. The values of all physical parameters are given in Table 2. The initial and boundary conditions are rather standard as well. Initially, at $t\ =\ 0\,$, the distribution of evolutionary quantities is given on the computational domain $[\,0,\,\ell\,]\,$:

[TABLE]

On the right boundary we set wall boundary conditions for simplicity999In any case, the simulation stops before the mass reaches the right boundary. So, the influence of the right boundary condition on presented numerical results is completely negligible.. On the left extremity of the computational domain the boundary conditions depend on the flow regime in the vicinity of the left boundary. If the flow there is supercritical, we impose Cauchy’s data. Otherwise, we impose a wall boundary condition as well. The case we study below is depicted in Figure 5 and it will be described in more details below.

Experimental set-up

In order to reproduce in silico the experiments of Rastello & Hopfinger (2004) [70] we use the following configuration of the numerical tank. The channel length $\ell$ is equal to $200\;{\sf cm}\,$. As we already mentioned, on boundaries we impose wall boundary conditions, in other words the channel is closed as in experiments. A heavier fluid of buoyancy b_{\,\ell}\ \in\ \bigl{\{}1,\,19\bigr{\}}\;\dfrac{{\sf cm}}{{\sf s}^{\,2}} fills an initially closed container of the length of $20\;{\sf cm}$ and of variable height $h_{\,\ell}\,$, which changed from one experiment to another. This configuration corresponds to the classical lock-exchange experiment. The heavy fluid “mass101010We employ the term mass in the sense of the relative weight of the dense fluid.” $m_{\,\ell}\ =\ b_{\,\ell}\cdot h_{\,\ell}\,$. At the distance of $50\;{\sf cm}$ from this recipient the slope is covered by initially motionless sediment layer of the height $\zeta_{\,s}\ =\ 0.2\;{\sf cm}$ and $2\;{\sf cm}\,$. The “mass” of sediments is $m_{\,s}\ =\ b_{\,s}\cdot\zeta_{\,s}\,$. The length of sediments layer is $100\;{\sf cm}\,$. During the propagation of the heavy fluid head all these sediments were entrained into the flow. In some experiments the sediment layer was absent, i.e.$m_{\,s}\ \equiv\ 0\,$. In this case the flow simply propagates over the rigid inclined bottom. In order to avoid the degeneration of certain equations, in numerical experiments we cover the whole slope with a micro-layer of sediments with mass $m_{\,s}^{\,\circ}\ >\ 0\,$, such that

[TABLE]

The values of all physical parameters are given in Table 2. The influence of the model parameters on predicted values was found to be rather weak. Consequently, in all numerical simulations reported below we used the following values of these parameters:

[TABLE]

The results of the critical comparisons with laboratory data are discussed below.

Numerical results

In Figure 6 we report the results of numerical simulations showing the spatial profiles of four quantities $h\,$, $u\,$, $m$ and $q\,$. These simulations were performed without the sediments layer, i.e.$m_{\,s}\ \equiv\ 0\,$. The panels correspond to:

**Left panel: **

Slope angle $\varphi\ =\ 32^{\circ}$, the snapshot is taken at $T\ =\ 10\;{\sf s}$

**Right panel: **

Slope angle $\varphi\ =\ 45^{\circ}$, the snapshot is taken at $T\ =\ 20\;{\sf s}\,$.

All quantities reported in Figure 6 are given in dimensional variables. The position of the density current front is clearly visible in both cases despite the fact that the IBVP was solved in the domain $[\,0,\,\ell\,]\times[\,0,\,T\,]$ using simple shock-capturing methods (*i.e.*no special treatment was necessary to detect the front). This property of the employed numerical scheme will be used to analyze below the asymptotic behaviour of unsteady solutions.

In Figure 7 we show the dependence of the flow head velocity $U_{\,f}$ on the distance $x_{\,f}$ traveled by the front. The parameters of these numerical/laboratory Experiments 1 & 2 are given in Table 2. Please, note that the heavy fluid density in Experiment 1 is much higher. The last observation explains why the head velocity in Experiments 1 is higher than in Experiments 2, even if the slope is bigger in Experiments $2\,$. Various symbols ($\bigcirc$ and $\bigtriangledown$) correspond to laboratory measurements of the head velocity and are taken from [70, Fig. 11]. Solid and dash-dotted lines (with $\bigcirc$ and $\bigtriangledown$) indicate the experiments with ($m_{\,s}\ >\ 0$) and without ($m_{\,s}\ =\ 0$) sediments correspondingly. The presence of sediments before the head front yields the initial acceleration of the flow. This acceleration phase is followed by deceleration, since sediments are distributed over a limited distance in our experiments. A very good agreement with numerical predictions can be observed in Figure 7.

In Figure 8 we show the long time behaviour of the front velocity from the Experiment $2$ reported in Figure 7 with parameters from Table 2. In order to perform this simulation we increased the computational domain111111This operation can be made easily in numerical experiments contrary to laboratory experiments, where the channel size is rather fixed. from $200$ to $500\;{\sf cm}\,$. Sediments were absent along the slope in this computation (i.e.$m_{\,s}\ \equiv\ 0$). The excellent agreement with theoretically established asymptotics $U_{\,f}\ \thicksim\ \dfrac{1}{\sqrt{x_{\,f}}}$ can be observed in this Figure 8. It validates both the numerics and the underlying theoretical argument.

Discussion

Above we have proposed and tested a new model for density currents. We outline below the main conclusions and perspectives of this study.

Conclusions

In the present study we considered the problem of density current modelling propagating down the slope in the presence of sediment deposits on the rigid bottom. In order to describe this flow in mathematical terms, we divide it into three zones (from the bottom vertically upwards): the layer of sediments (i), the mixing layer (ii) and the still water layer (iii). In order to simplify the mathematical description, we assumed that the upper layer (iii) of the still water was motionless. The depth-averaged description was adopted in each of two remaining layers (i) and (ii). As a result, we arrived to a shallow water two-layer system including turbulent modelling, which can be recast in a conservative form (2.1) – (2.6) of quasi-linear balance laws [72]. The initially proposed system contains only six evolution equations to describe the complex density current. However, it can be further simplified if we assume that the flow in the bottom layer occupied by dense liquid or suspension with uplifted sediments is in the equilibrium state, *i.e.*the gravity force is exactly balanced by the friction forces with the rigid bottom. In this way, we can remove two equations corresponding to the bottom layer (i). The sediment equilibrium model is given by System (4.1) – (4.4). This concludes the modelling part of our study.

Then, the proposed equilibrium model is studied using analytical means. Namely, we were interested in some special, but very important classes of solutions — steady states, travelling waves and similarity flows. The linear stability of travelling waves was studied. The initial conditions which yield asymptotically (in time) some special solutions (such as self-similar and travelling waves) were discussed. The velocity of travelling waves is of capital importance for the understanding of turbidity flows propagation. Namely, the travelling wave speed gives an estimation of the flow head propagation along the slope. This prediction of our travelling wave analysis was checked against the experimental data from numerous previous studies [26, 81, 76, 9].

The mixing layer formation was discussed in the framework of steady solutions and our theoretical prediction was compared against experimental data from [67]. Moreover, unsteady solutions predicted by our model were compared with density flow experiments made in the LEGI (UMR 5519) Laboratory by Rastello & Hopfinger (2004) [70]. The model predictions are in good agreement with their measurements of the front velocity. Finally, the asymptotic behaviour for the flow head velocity as a function of the front position was measured in our numerical computations and an excellent agreement was found with existing theoretical estimations. This concludes the validation of the model and also of several particular solutions derived in this study.

It is important to underline that the proposed model describes three different flow regimes:

1. (1)

Exchange flow, 2. (2)

Thermals, 3. (3)

Underwater avalanches.

In the first case, the permanent (or transient) inflow of a heavy fluid mass accelerates downslope. In natural environments this situation can be realized when a heavy fluid overflows a bottom obstacle or a hill, which creates a deposit of the heavy fluid mass behind the obstacle. Our mathematical model showed that the initial accelerating phase in the development of the supercritical flow is of capital importance since it yields the formation of the so-called mixing layer. In particular, the density current head velocity depends directly on the portion of the mass entering the bottom layer. It is this sea-bed underflow, which supplies the head directly almost without mixing. This fact may explain also the relatively important scattering of observed front velocities in laboratory experiments (see *e.g.*Figure 4).

In order to describe the evolution of thermals over a downslope, one can use the simplified model (2.11) – (2.14), which is a direct generalization of usual nonlinear shallow water equations taking into account the ambient quiescent fluid entrainment into the flow. It is important to notice also that the entrainment speed in our model is determined automatically together with other flow parameters. Figure 8 showed also that this model can be used to describe thermal flows during the acceleration and deceleration phases simultaneously.

Density currents over a sloping bottom featuring sediments entrainment into the flow have numerous applications in geophysics such as snow and underwater avalanches, seasonal mass exchanges in deep lakes, *etc.*These flows are situated in the center of different mathematical modelling efforts mentioned in Introduction. Unfortunately, contrary to the previous two situations (lock-exchange and thermal flows), laboratory experiments and real-world measurements still remain seldom. It is in this field of applications that the rôle of simple mathematical models increases substantially. Such models should allow to predict theoretically at least the front velocity and, if possible, several other parameters of the self-sustained current. In the particular situation where we neglect the sedimentation velocity, the one layer model (2.11) – (2.14) is suitable for the description of underwater avalanches. Moreover, our computation showed that the flow achieves quickly a self-similar regime, which will be studied in future publications. However, in contrast to exchange flows, the condition (3.30) of the flow criticality behind the wave front (analogous to the Chapman–Jouguet condition in detonation) is not fulfilled in the presence of sediments entrainment. It is another flow regime, which is realized in practice: the flow turns out to be substantially supercritical in the frame of reference moving with the front. This situation is usually referred to as “overdriven detonation”, when the main flow parameters, such as detonation pressure and propagating velocity, exceed the corresponding Chapman–Jouguet values [50]. A further study is necessary to understand the front velocity selection mechanism for such flows. In this way we come naturally to the description of some other perspectives of our study.

Perspectives

In the present study, we have employed only basic legacy first order Godunov finite volume schemes [28, 30, 29]. Their principal advantages are the ease of implementations and the robustness of numerical results. However, the accuracy might be also important in some applications, where a quantitative prediction is critically important. Consequently, some high order finite volume well-balanced schemes have to be developed to solve numerically the density current models proposed in our study. This technology is relatively well mastered nowadays [16, 14, 15, 31].

The system of equations presented earlier in this article is a simplified equilibrium model based on a number of physical assumptions. Of course, more complete models of higher physical fidelity should be developed as well in the future along the lines outlined in the monograph [48]. For instance, in this manuscript we considered density currents with uniform cross-sections, *i.e.*2D flows. In future works 3D effects have to be included. Moreover, in all examples considered earlier, the uniform bed slopes was used. In part, it comes from the experimental set-up used in previous investigations. However, in upcoming works the interaction between the fluid flow with the bed morphology has to be investigated to advance our understanding of these processes.

Acknowledgments

This research was supported by the project Aval (AAP Montagne 2015, University Savoie Mont Blanc) and Russian Science Foundation (project 15-11-20013). V. Liapidevskii acknowledges the CNRS/ASR Cooperation Program under the project №23975 and the hospitality of the University Savoie Mont Blanc during his numerous visits. The authors are grateful to Professors L. Armi & G. Pawlak for their help in preparing this manuscript.

Appendix A An exact solution

In the Appendix we present the construction of a particular steady solution. Using the information presented above in Sections 3.2 and 3.3, we construct in this Section the profile of a particular exact solution to Equations (2.7) – (2.10). We assume that in the bottom layer $\zeta\ \equiv\ \zeta_{\,b}$ and $w\ \equiv\ w_{\,b}\,$. Moreover, we suppose that the total mass inflow $\mathcal{M}_{\,0}\ =\ b_{\,0}\,\zeta_{\,0}\,w_{\,0}$ is prescribed on the left boundary. This incident mass flux $\mathcal{M}_{\,0}$ splits into two streams:

[TABLE]

We assume also that the ratio $\mu\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ \dfrac{\mathcal{M}_{\,j}}{\mathcal{M}_{\,0}}\ \in\ (0,\,1\,]$ is fixed. Thus, $\mathcal{M}_{\,b}$ can be also easily deduced.

The flow at any fixed moment of time $t\ =\ t_{\,m}\ >\ 0$ is composed of two parts:

1. (1)

Stationary flow: $u\,(x,\,t)\ \equiv\ u_{\,j}\,$, $m\,(x,\,t)\ \equiv\ m_{\,j}\,$, $q\,(x,\,t)\ \equiv\ q_{\,j}\,$ and $h\,(x,\,t)\ =\ h_{\,j}\ +\ \varsigma_{\,j}\cdot x\,$, where $x\ \in\ [\,0,\,L_{\,j}\,]\,$, $L_{\,j}\ =\ u_{\,j}\cdot t_{\,m}\,$. 2. (2)

Travelling wave: This part depends on the variable $\xi\ =\ x\ -\ c\cdot t$ and $u\,(\xi)\ \equiv\ u_{\,f}\,$, $m\,(\xi)\ \equiv\ m_{\,f}\,$, $q\,(\xi)\ \equiv\ q_{\,f}\,$ and $h\,(\xi)\ =\ h_{\,f}\ -\ \varsigma_{\,f}\cdot(x\ -\ L_{\,t})\,$, where $\varsigma_{\,f}\ >\ 0\,$, $L_{\,t}\ -\ L_{\,f}\ \leqslant\ x\ \leqslant\ L_{\,t}$ and $L_{\,f}\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ u_{\,f}\cdot t_{\,m}\,$, $L_{\,t}\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ c\cdot t_{\,m}\,$.

The problem consists in finding all solution parameters. The solution algorithm is summarized below.

Steady flow

Let $a\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ \dfrac{u_{\,j}}{U_{\,j}}$ be the ratio of two speeds with $U_{\,j}\ =\ \sqrt[3]{\mathcal{M}_{\,j}}$ and $z\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ a^{3}\,$. Then, one has to find the unique root $z\ >\ 1$ to the following equation:

[TABLE]

which is to be compared with Equation (3.27). It is equivalent to finding the root to the following cubic polynomial:

[TABLE]

Let $z^{\star}$ be the required root. Then, we find $a\ =\ \sqrt[3]{z^{\star}}$ and $u_{\,j}\ =\ a\,U_{\,j}\,$. We know that $\mathcal{M}_{\,j}\ =\ U_{\,j}^{\,3}\ =\ m_{\,j}\,u_{\,j}\,$. Hence, $m_{\,j}\ =\ \dfrac{U_{\,j}^{\,2}}{a}\,$. The remaining variables are:

[TABLE]

Finally, the steady part of the turbulent jet layer depth is equal to

[TABLE]

Travelling waves

For the travelling wave part we have $\mathcal{M}_{\,b}\ =\ U_{\,b}^{\,3}$ and $\mathrm{Fr}_{\,b}\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ \dfrac{w_{\,b}}{\sqrt{b_{\,0}\,\zeta_{\,b}}}\,$. Let $u^{\,\star}\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ \frac{u_{\,f}}{c}$ be the dimensionless flow velocity. It can be found by solving the following algebraic equation, which possesses a real root on the interval u^{\,\star}\ \in\ \bigl{(}{\textstyle{1\over 2}},\,1\bigr{)}\,:

[TABLE]

where the constant $\upbeta$ was defined earlier in (3.27). The polynomial $\mathcal{P}_{\,2}\,(u^{\,\star})$ can be expanded:

[TABLE]

Let us assume that we have found the required root u^{\,\star}\ \in\ \bigl{(}{\textstyle{1\over 2}},\,1\bigr{)}\,. Then, we solve another polynomial equation to determine $C\ =\ \frac{c}{U_{\,b}}\,$:

[TABLE]

The last equation admits a unique positive root $C^{\,\star}\ =\ C^{\,\star}\,(\mathrm{Fr}_{\,b})\ =\ C_{\,e}\,(\varphi)$ as well. The last dependence is shown in Figure 4 with the dashed line. Then, we consider the flow in the bottom layer. We have the following relations:

[TABLE]

[TABLE]

We can notice that $w_{\,b}$ can be also written as $w_{\,b}\ =\ \mathrm{Fr}_{\,b}^{\,2/3}\,U_{\,b}\,$. Now we come back to the travelling wave:

[TABLE]

and from the relation

[TABLE]

we can find

[TABLE]

provided that $(c\ -\ u_{\,f})\cdot b_{\,0}\ \neq\ 0\,$. Finally, we find two last elements of the solution:

[TABLE]

The travelling wave profile is then given by

[TABLE]

The last profile is located in the segment $x\ \in\ [\,L_{\,t}\ -\ L_{\,f},\,L_{\,t}\,]\,$. The solution exists if $u_{\,j}\ \leqslant\ u_{\,f}\,$. We underline that in the region between the stationary and travelling wave parts the flow in general is unsteady. Thus, the indicated boundaries $L_{\,t}\ -\ L_{\,f}$ and $L_{\,t}$ are rather approximations to the reality. Two such exact solutions for two different values of the parameter $\mu\ =\ 0.05$ (left panel) and $\mu\ =\ 0.654$ (right panel) are depicted in Figure 9. This picture shows that at least two different configurations can be realized in practice:

**Two-wave (left panel (a)): **

when the main portion of the heavy fluid enters into the boundary layer ($\mu\ =\ 0.05$)

**One-wave (right panel (b)): **

when the flow in the buoyant jet reaches the head front ($\mu\ =\ 0.654$).

Remark 2**.**

In experiments in order to highlight the processes happening in the flow, the heavy fluid of density $\rho_{\,0}$ is in general coloured. Other experimental visualization techniques are also available, but the fluid painting is the most widely used one. By $\rho_{\,a}$ we denote the density of the light ambient fluid. In such experimental conditions the visible interface between two fluids is located, where the following condition is satisfied:

[TABLE]

where $\rho_{\,v}$ is the visible density and $\varpi$ is a constant whose approximate value belongs to the segment $[\,0.01,\,0.1\,]\ \ni\ \varpi\,$. For the sake of illustration, we depicted this visible interfaces in Figure 9 with solid lines, while the exact analytical solutions were shown with dashed lines. In the preparation of the visible interface we assumed that the fluid density varies linearly inside the turbulent layer from $\rho_{\,a}$ to $\rho_{\,0}\,$. In other words, the virtual visible interface $y\ =\ y_{\,v}$ is determined by the following relation:

[TABLE]

For the sake of convenience, the last condition can be equivalently reformulated in terms of the buoyancy variable:

[TABLE]

Appendix B Derivation of the speed-Froude relation

From “boundary” conditions (3.29), (3.30) we deduce the following system of two equations with respect to unknowns $u$ and $c\,$:

[TABLE]

By solving this system we obtain the following solution ($u\ \neq\ 1$):

[TABLE]

Furthermore, by definition we have

[TABLE]

Thus,

[TABLE]

Hence,

[TABLE]

Taking into account the last results, the “boundary” condition (3.28) becomes:

[TABLE]

Or equivalently we have:

[TABLE]

By rearranging the terms in the last equation we obtain:

[TABLE]

By introducing the dimensionless velocity $C\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ \frac{c}{U_{\,b}}$ we obtain the required Equation (3.31).

Appendix C Nomenclature

In the main text above we used the following notations (this list is not exhaustive):

\equiv\**: **

equal identically

\propto\**: **

proportional

\mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\**: **

the left hand side is defined

\mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{=:}}\,}\**: **

the right hand side is defined

\apprle\**: **

smaller than the approximate upper bound121212We underline the fact that this bound is given approximatively.

\apprge\**: **

greater than the approximate lower bound††footnotemark:

\varphi\**: **

angle of the bottom slope

\upalpha\**: **

local bottom slope, i.e. $\upalpha\ =\ \tan\varphi$

h_{\,j}\**: **

total depth of the layer $j$

\rho_{\,j}\**: **

fluid density in the layer $j$

\zeta\**: **

total depth of the bottom layer

b_{\,j}\**: **

buoyancy of the layer $j$

u_{\,j}\**: **

depth-averaged horizontal velocity of fluid particles in the layer $j$

v\**: **

vertical velocity component

w\**: **

horizontal velocity component in the bottom layer

\mathrm{w}\**: **

transversal velocity component (in 3D)

p\**: **

fluid pressure

m_{\,j}\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ b_{\,j}\,h_{\,j}\**: **

“mass” contained in the fluid column

\mathcal{M}_{\,j}\ \mathop{\stackrel{{\scriptstyle\,\mathrm{def}}}{{:=}}\,}\ m_{\,j}\,u_{\,j}\**: **

mass flux in the layer $j$

q\**: **

depth-averaged turbulent mean square velocity in the mixing layer

\xi\**: **

characteristics speed (*i.e.*eigenvalues of the Jacobian matrix of the hyperbolic fluxes)

\lambda_{\,1,\,2}\**: **

eigenvalues in the stability studies

c\**: **

dimensional velocity of the travelling wave

C\**: **

dimensionless wave velocity

\mathrm{Fr}\**: **

the dimensionless Froude number

\ell,\,L\**: **

length scales

t\**: **

time variable

x\**: **

“horizontal” coordinate along the bottom slope

y\**: **

“vertical” coordinate normal to the bottom

d\,(x)\**: **

the bathymetry (depth) function

g\**: **

gravity acceleration

\tilde{g}\**: **

reduced gravity acceleration

\chi\**: **

entrainment rate

u_{\,\ast}\**: **

friction velocity at the solid bottom