# Long nonlinear internal waves and mixing in a three-layer stratified flow in
the Boussinesq approximation

A. A. Chesnokov1<EMAIL_ADDRESS>S. L. Gavrilyuk2  V. Yu.
Liapidevskii1  1Lavrentyev Institute of Hydrodynamics, 15 Lavrentyev prospect,
630090 Novosibirsk, Russia 2Aix-Marseille Université, UMR CNRS 7343, IUSTI, 5
rue E. Fermi, 13453 Marseille CEDEX 13, France

###### Abstract

A one-dimensional long-wave model of an unsteady three-layer flow of a
stratified fluid under a lid is proposed, taking into account turbulent mixing
in the intermediate layer. In the Boussinesq approximation, the equations of
motion are reduced to an evolutionary system of balance laws, which is
hyperbolic for a small difference in velocities in the layers. Classes of
stationary solutions are studied and the concept of subcritical
(supercritical) three-layer stratified flow is introduced. Oscillating
solutions are constructed that describe the spatial evolution of the mixing
layer. The problem of transcritical flow over an obstacle is considered.
Solutions are obtained that describe qualitatively different flow regimes on
the leeward side of the obstacle. The proposed model is validated using
experimental data and field observations on the entrainment of ambient fluid.

###### keywords:

internal waves; stratified flows; turbulent mixing

##  1 Introduction

Stratified flows over topography are ubiquitous in geophysical and
environmental settings. For instance, the flow of a dense fluid along a
sloping bottom occurs in deep ocean outflows and severe downslope winds. A
detailed overview of the diversity of such flows and their relevance can be
found in the monographs of Baines (2005) and Simpson (1997).

A distinguishing feature of stratified flow over topography is the formation
of a bifurcation enclosing partly mixed fluid. Detailed observations of the
establishment of such topographic flows over the Knight Inlet sill (British
Columbia, Canada) are described in Farmer & Armi (1999); Armi & Farmer (2002);
Cummins et al. (2006). In particular, in the last reference, the upstream
formation of a strong internal hydraulic jump is reported and studied. Similar
phenomena are reproduced in the laboratory experiments of Pawlak & Armi (1998,
2000) where the special character of shear instability in steep downslope
flows is illustrated. The enhanced entrainment efficiency due to the
Kelvin–Helmholtz instability leads to the changing density of the trapped
fluid, thus providing the link between small-scale processes and the larger-
scale response. A Navier–Stokes type model of stratified flows is used in Lamb
(2004) to describe the formation of a strong supercritical flow beneath a
large breaking lee wave in tidal flow over the Knight Inlet sill. Further
study of topographic effects in stratified flows is carried out Winters & Armi
(2014). A set of example flows is introduced that systematically leads to the
observation that approximately blocked and topographically controlled flows
are produced whenever a uniform inflow encounters an obstacle that is
sufficiently tall. The dynamic stability and connection between topographic
control and wave excitation aloft of stratified flow configurations
characteristic of hydraulically controlled downslope flow over topography are
studied in recent works of Jagannathan et al. (2017, 2020).

In theoretical modelling of turbulent shear flows, such as a surface or near-
bottom jet, or a mixing layer, an important question is how turbulence is
generated and how it maintains itself in a stably gravity-stratified shear
flow. Chu & Baddour (1984) considered the characteristic features of the
stratification influence on mixing processes in the near-field and far-field.
Sher & Woods (2015, 2017) performed experiments of a turbulent gravity current
and presented measurements of the entrainment of ambient fluid into gravity
currents produced by a steady flux of buoyancy. It is known that the process
of entrainment produces a deepening mixing layer at the interface, which
increases the gradient Richardson number of this layer and eventually may
suppress further entrainment. Taking this fact into account, Horsley & Woods
(2018) constructed an analytical solution to the simplified depth-averaged
model and discussed its properties. Yuan & Horner-Devine (2017) presented
experimental investigation of large-scale vortices in a freely spreading
gravity currents propagating into laterally confined and unconfined
environments.

Internal turbulent hydraulic jump is another possible mechanisms of mixing.
Ogden & Helfrich (2016) investigated internal hydraulic jumps in flows with
upstream shear using two-layer shock-joining theories. The models can be
modified to allow entrainment and to indirectly account for continuous
velocity profiles, producing solutions where the basic theories failed. Ogden
& Helfrich (2020) studied these jumps to illuminate the changing physics as
shear increases. Baines (2016) presented a two-layer model of internal
hydraulic jumps which incorporates mixing between the layers within the jump.
Thorpe & Li (2014); Thorpe et al. (2018) derived another model of internal
hydraulic jumps which adopts continuous profiles of velocity and density both
upstream and downstream of the transition. They also applied this model to
describe observations made at several locations in the ocean.

The aim of the present paper is to propose a mathematical model for the
formation and evolution of a mixing layer resulting from the interaction
between two co-flowing layers of homogeneous fluids over topography.
Experimental data and field observations on gravity currents presented in the
above-mentioned works indicate the presence of a fairly clear boundaries
between the regions of potential flow of homogeneous outer layers with
different densities and an intermediate non-homogeneous layer of turbulent
mixing. The outer layers are potential and can be approximately described by a
shallow water type model. In the intermediate layer, the flow is sheared and
is described by the equations of weakly sheared flows (Teshukov (2007)). The
interaction between the outer layers and interlayer is taken into account
through a natural mixing process, where the mixing velocity is proportional to
the intensity of large eddies in the interlayer. This approach has already
been used to simulate a mixing layer in a homogeneous fluid (Liapidevskii &
Chesnokov (2014); Chesnokov & Liapidevskii (2020)) and to describe the
evolution of spilling breakers in shallow water (Gavrilyuk et al. (2016)). Its
extension to the case of inhomogeneous stratified flows is proposed by
Gavrilyuk et al. (2019). A similar approach based on the layered description
of the stratified flow is used by Liapidevskii (2004); Liapidevskii et al.
(2018) to simulate the mixing layer on the leeward side of the obstacle and
sediment laden gravity currents. However, in the last references, the process
of entrainment of liquid into the mixing layer is taken into account in a
slightly different way.

In the present work, we assume that the upper boundary is fixed and the
density difference in the layers is small, so we can use the Boussinesq
approximation. This makes it possible to obtain a depth-averaged model of a
three-layer system consisting of two outer homogeneous layers and an
intermediate mixing layer in the form of a system of six balance laws. The
system is hyperbolic for small relative velocities in the layers. This model
allows for a simple numerical implementation and is used to describe the basic
mixing features in stratified flows over a flat and uneven bottom.

The three-layer flow scheme makes it possible to eliminate contradictions that
cannot be resolved in the framework of a two-layer scheme. Among them a
problem of choosing relations between the downstream and upstream states of an
internal hydraulic jump, as well as the need to take into account the non-
uniformity of velocity profile and mass transfer between the layers by
additional empirical relations for hydraulic jumps (Ogden & Helfrich (2016)).
The model derived in this paper is a closed system of integral conservation
laws describing a three-layer shallow water flow over topography. The model
allows one to describe turbulent mixing between homogeneous layers as a
nonlinear stage of the Kelvin–Helmholtz instability development at the
interface between layers both in stationary and non-stationary flows. The non-
uniform velocity profile (long-wave horizontal vorticity) in the mixing layer
is taken into account. The main attention in the work is focused on the study
of the capabilities of the presented model for the mathematical description of
the effects of mixing and topography on the structure of supercritical and
subcritical flows of stratified fluids. Applications of this model for the
quantitative description of laboratory experiments and field observations are
presented to show its ability to describe complex stratified flows.

In the following section, a one-dimensional long-wave model of a three-layer
stratified flow under a lid is derived. The model takes into account turbulent
mixing between adjacent layers. In Section 3, the Boussinesq approximation is
performed and the characteristic velocities are determined. Stationary
solutions are studied in Section 4, where the definition of subcritical
(supercritical) three-layer flows is proposed. Then, continuous and
discontinuous oscillating solutions are constructed that describe the spatial
evolution of the mixing layer. In Section 5, the model is validated by
comparison with known experimental data and field observations. Finally,
conclusions are drawn in Section 6.

##  2 Long-wave model of a three-layer flow with mixing

Consider two-dimensional non-stationary flows of an inviscid incompressible
non-homogeneous fluid confined between the moving bottom $z=Z(t,x)$ and rigid
upper lid $z=H_{0}={\rm const}$. The corresponding dimensionless Euler
equations are

$\begin{array}[]{l}\displaystyle
u_{x}+w_{z}=0,\quad\rho_{t}+(u\rho)_{x}+(w\rho)_{z}=0,\\\\[8.53581pt]
\displaystyle(\rho u)_{t}+(\rho u^{2})_{x}+(\rho
uw)_{z}+p_{x}=0,\\\\[8.53581pt] \displaystyle\varepsilon^{2}\big{(}(\rho
w)_{t}+(\rho uw)_{x}+(\rho w^{2})_{z}\big{)}+p_{z}=-g\rho.\end{array}$ (1)

The kinematic boundary conditions at $z=Z$ and $z=H_{0}$ are

$Z_{t}+uZ_{x}-w\big{|}_{z=Z}=0,\quad w\big{|}_{z=H_{0}}=0.$ (2)

In what follows, we assume that the flow has a three-layer structure (Fig. 1).
The internal boundaries $z=z_{1}=Z+h_{1}(t,x)$ and $z=z_{2}=H_{0}-h_{2}(t,x)$
separate the outer layers, where the flow is almost potential and homogeneous
with constant densities $\rho_{1}$ and $\rho_{2}$, and the intermediate
turbulent non-homogeneous layer with density $\rho(t,x,z)$. The vertical
density distribution is assumed to be continuous so that

$\rho\big{|}_{z=z_{1}+0}=\rho_{1},\quad\rho\big{|}_{z=z_{2}-0}=\rho_{2}.$ (3)

The depth-averaged density of the intermediate layer is defined as

$\bar{\rho}(t,x)=\frac{1}{\eta}\int_{z_{1}}^{z_{2}}\rho(t,x,z^{\prime})dz^{\prime}.$

We will show later that the vertical density distribution will be stable for
any time, i.e. $\rho_{2}<\bar{\rho}(t,x)<\rho_{1}$, if it was initially
stable. At the internal boundaries the kinematic conditions are satisfied:

$z_{1t}+uz_{1x}-w\big{|}_{z=z_{1}}=-M_{1},\quad
z_{2t}+uz_{2x}-w\big{|}_{z=z_{2}}=M_{2}.$ (4)

The right-hand sides $M_{1}$ and $M_{2}$ responsible for the mixing between
layers will be precised later. In the above equations $u=V^{-1}\hat{u}$,
$w=(d_{0}V)^{-1}l_{0}\hat{w}$, $p=(\rho_{0}V^{2})^{-1}\hat{p}$,
$x=l_{0}^{-1}\hat{x}$, $z=d_{0}^{-1}\hat{z}$, $t=Vl_{0}^{-1}\hat{t}$, and
$g=d_{0}V^{-2}\hat{g}$ are dimensionless components of the velocity, pressure,
Cartesian coordinates, time, and gravitational acceleration, respectively;
$\hat{u}$, $\hat{w}$, $\hat{p}$, $\hat{x}$, $\hat{z}$, $\hat{t}$, and
$\hat{g}$ are the corresponding dimensional variables. The parameters $V$,
$\rho_{0}$, $d_{0}$, $l_{0}$ denote the characteristic velocity, density, and
the characteristic vertical and horizontal scales, respectively. The
dimensionless parameter $\varepsilon=d_{0}/h_{0}$ is the ratio of vertical and
horizontal scales. System (1) admits the conservation of energy

$E_{t}+\big{(}(E+p)u\big{)}_{x}+\big{(}(E+p)w\big{)}_{z}=0,$ (5)

where $E=(u^{2}/2+\varepsilon^{2}w^{2}/2+gz)\rho$. We suppose that the waves
are long, i.e. $\varepsilon\ll 1$ and the terms of order $O(\varepsilon^{2})$
are neglected in equations (1), (5). In this case, the last equation in system
(1) is reduced to the hydrostatic law of pressure distribution over the depth

$p=p^{*}(t,x)+g\int\limits_{z}^{H_{0}}\rho(t,x,z^{\prime})\,dz^{\prime},$ (6)

where $p^{*}$ is the pressure at the upper lid.

Figure 1: Three-layer stratified flow over topography.

To describe the three-layer flow, we define the following depth-averaged
variables:

$\begin{array}[]{l}\displaystyle
u_{1}=\frac{1}{h_{1}}\int\limits_{Z}^{z_{1}}u\,dz,\quad
u_{2}=\frac{1}{h_{2}}\int\limits_{z_{2}}^{H_{0}}u\,dz,\quad\bar{u}=\frac{1}{\eta}\int\limits_{z_{1}}^{z_{2}}u\,dz,\\\\[8.53581pt]
\displaystyle\quad\quad\quad
q^{2}=\frac{1}{\eta}\int\limits_{z_{1}}^{z_{2}}(u-\bar{u})^{2}\,dz,\quad\bar{\rho}=\frac{1}{\eta}\int\limits_{z_{1}}^{z_{2}}\rho\,dz.\end{array}$
(7)

Here $\eta=z_{2}-z_{1}$ is the thickness of the intermediate layer. The
variable $q$ measures the distortion of the velocity profile in the mixing
layer. Note that the geometric constraint yields $h_{1}+\eta+h_{2}=H_{0}-Z$.
Further we derive a one-dimensional closed system of depth-averaged equations.

###  2.1 Depth-averaged equations for the outer layers

As mentioned above, in the lower and upper layers the fluid is homogeneous
($\rho_{i}={\rm const}$, $i=1,2$). From (6) it follows that the pressure in
these layers is

$p=\left\\{\begin{array}[]{ll}(z_{1}-z)g\rho_{1}+g\bar{\rho}\eta+g\rho_{2}h_{2}+p^{*},&\
z\in(Z,z_{1})\\\\[5.69054pt] (H_{0}-z)g\rho_{2}+p^{*}&\
z\in(z_{2},H_{0})\end{array}\right.$ (8)

The derivation of averaged over depth equations for outer layers is similar.
For definiteness, we will derive in details the equations for the upper layer.
Integrating the incompressibility, horizontal momentum and energy equations
with respect to $z$ over the interval $(z_{2},H_{0})$ and using the boundary
conditions, we get the following exact integral relations:

$\begin{array}[]{l}\displaystyle
h_{2t}+(u_{2}h_{2})_{x}=-M_{2},\\\\[8.53581pt]
\displaystyle(u_{2}h_{2})_{t}+\bigg{(}\int\limits_{z_{2}}^{H_{0}}\Big{(}u^{2}+\frac{p}{\rho_{2}}\Big{)}\,dz\bigg{)}_{x}=-\Big{(}M_{2}u+\frac{p}{\rho_{2}}z_{2x}\Big{)}\Big{|}_{z=z_{2}+0}\,,\\\\[8.53581pt]
\displaystyle\bigg{(}\int\limits_{z_{2}}^{H_{0}}E\,dz\bigg{)}_{t}+\bigg{(}\int\limits_{z_{2}}^{H_{0}}(E+p)u\,dz\bigg{)}_{x}=-\big{(}M_{2}(E+p)+h_{2t}p\big{)}\big{|}_{z=z_{2}+0}\,.\end{array}$
(9)

The calculation of the integrals in system (9) is based on the estimates given
below. We will say that the flow is weakly sheared if
$u_{z}=O(\varepsilon^{\alpha})$, $\alpha>0$. For a homogeneous fluid the flow
vorticity $\omega=u_{z}-\varepsilon^{2}w_{x}$ conserves along the
trajectories: $\omega_{t}+u\omega_{x}+w\omega_{z}=0$. The long-wave vorticity
$\omega=u_{z}$ also satisfies this equation with error $O(\varepsilon^{2})$.
Therefore, if $u_{z}=O(\varepsilon^{\beta})$ for $t=0$, then
$u_{z}=O(\varepsilon^{\beta})$ for any $t>0$, with $\beta=\min(\alpha,2)$. One
can prove that, if the flow is weakly sheared, then

$\frac{1}{h_{1}}\int\limits_{Z}^{z_{1}}u^{2}\,dz=u_{1}^{2}+O(\varepsilon^{2\beta}),\quad\frac{1}{h_{2}}\int\limits_{z_{2}}^{H_{0}}u^{3}\,dz=u_{2}^{3}+O(\varepsilon^{3\beta}).$
(10)

(for the proof, see Barros et al. (2007); Gavrilyuk et al. (2019)).

We say that the fluid flow in the outer layers is almost potential, if
$\alpha>1$. Since terms of order $O(\varepsilon^{2})$ and higher are not taken
into account in the considered long-wave approximation, this means the terms
$O(\varepsilon^{2\beta})$ and $O(\varepsilon^{3\beta})$ in formulae (10) are
ignored. Neglecting higher order terms system (9) takes the form

$\begin{array}[]{l}\displaystyle
h_{2t}+(u_{2}h_{2})_{x}=-M_{2},\\\\[8.53581pt]
\displaystyle(u_{2}h_{2})_{t}+\Big{(}u_{2}^{2}h_{2}+\frac{gh_{2}^{2}}{2}+\frac{p^{*}h_{2}}{\rho_{2}}\Big{)}_{x}=-M_{2}u\big{|}_{z=z_{2}+0}+\Big{(}\frac{p^{*}}{\rho_{2}}+gh_{2}\Big{)}h_{2x}\,,\\\\[8.53581pt]
\displaystyle\bigg{(}\frac{u_{2}^{2}h_{2}}{2}+gH_{0}h_{2}-\frac{gh_{2}^{2}}{2}\bigg{)}_{t}+\bigg{(}\frac{u_{2}^{3}h_{2}}{2}+\Big{(}\frac{p^{*}}{\rho_{2}}+gH_{0}\Big{)}u_{2}h_{2}\bigg{)}_{x}\\\\[8.53581pt]
\displaystyle\quad\quad=-\bigg{(}\frac{u^{2}}{2}\bigg{|}_{z=z_{2}+0}+gH_{0}+\frac{p^{*}}{\rho_{2}}\bigg{)}M_{2}-\Big{(}\frac{p^{*}}{\rho_{2}}+gh_{2}\Big{)}h_{2t}\,.\end{array}$
(11)

Then we transform the energy equation of the system:

$\begin{array}[]{l}\displaystyle\bigg{(}\frac{u_{2}^{2}}{2}+gH_{0}+\frac{p^{*}}{\rho_{2}}\bigg{)}\big{(}h_{2t}+(u_{2}h_{2})_{x}\big{)}+u_{2}h_{2}\Big{(}u_{2t}+u_{2}u_{2x}+\frac{1}{\rho_{2}}p^{*}_{x}\Big{)}\\\\[8.53581pt]
\displaystyle\quad\quad=-\bigg{(}\frac{u^{2}}{2}\Big{|}_{z=z_{2}+0}+gH_{0}+\frac{p^{*}}{\rho_{2}}\bigg{)}M_{2}\end{array}$

Using the first and second equations in (11), we get

$\frac{M_{2}}{2}(u_{2}-u)^{2}\big{|}_{z=z_{2}+0}=0.$ (12)

Similarly, integrating equations (1), (5) and using boundary conditions (2),
(4) and the pressure representation (8), we obtain the depth-averaged
equations for the lower layer:

$\begin{array}[]{l}\displaystyle
h_{1t}+(u_{1}h_{1})_{x}=-M_{1},\\\\[8.53581pt]
\displaystyle(u_{1}h_{1})_{t}+\Big{(}u_{1}^{2}h_{1}+\frac{gh_{1}^{2}}{2}+(p^{*}+g\rho_{2}h_{2}+g\bar{\rho}\eta)\frac{h_{1}}{\rho_{1}}\Big{)}_{x}\\\\[8.53581pt]
\displaystyle\quad=-M_{1}u\big{|}_{z=z_{1}-0}+\frac{1}{\rho_{1}}(p^{*}+g\rho_{2}h_{2}+g\bar{\rho}\eta)h_{1x}-gh_{1}Z_{x}\,,\\\\[8.53581pt]
\displaystyle\bigg{(}\frac{u_{1}^{2}h_{1}}{2}+gZh_{1}+\frac{gh_{1}^{2}}{2}\bigg{)}_{t}+\bigg{(}\frac{u_{1}^{3}h_{1}}{2}+\Big{(}\frac{p^{*}}{\rho_{1}}+gh_{2}\frac{\rho_{2}}{\rho_{1}}+g\eta\frac{\bar{\rho}}{\rho_{1}}+gz_{1}\Big{)}u_{1}h_{1}\bigg{)}_{x}\\\\[8.53581pt]
\displaystyle\quad=-\frac{M_{1}u^{2}}{2}\Big{|}_{z=z_{1}-0}-gz_{1}M_{1}-(M_{1}+h_{1t})\bigg{(}\frac{p^{*}}{\rho_{1}}+gh_{2}\frac{\rho_{2}}{\rho_{1}}+g\eta\frac{\bar{\rho}}{\rho_{1}}\bigg{)}+gh_{1}Z_{t}.\end{array}$
(13)

Taking into account the first two equations of (13), the energy equation
yields:

$\frac{M_{1}}{2}(u_{1}-u)^{2}\big{|}_{z=z_{1}-0}=0.$ (14)

Suppose that the variables $M_{1}$ and $M_{2}$ are not identically equal to
zero. Then, by virtue of (14) and (12), the compatibility condition between
the energy, momentum and mass equations for flows in the outer layers gives us
only one possibility

$u|_{z=z_{1}-0}=u_{1},\quad u|_{z=z_{2}+0}=u_{2}.$ (15)

Below we assume that conditions (15) are satisfied.

###  2.2 Depth-averaged equations for the intermediate non-homogeneous shear
layer

As we already said above, the velocity and density profiles at the internal
interfaces are supposed to be continuous, i.e.

$(u,\rho)\big{|}_{z=z_{1}+0}=(u_{1},\rho_{1}),\quad(u,\rho)\big{|}_{z=z_{2}-0}=(u_{2},\rho_{2}).$
(16)

Integrating the first three equations in system (1) and energy equation (5)
with respect to $z$ over the intermediate layer thickness and using boundary
conditions (4), (3), (15), and (16), we get:

$\begin{array}[]{l}\displaystyle\eta_{t}+(\bar{u}\eta)_{x}=M_{1}+M_{2},\quad(\bar{\rho}\eta)_{t}+\bigg{(}\int\limits_{z_{1}}^{z_{2}}u\rho\,dz\bigg{)}_{x}=\rho_{1}M_{1}+\rho_{2}M_{2},\\\\[8.53581pt]
\displaystyle\bigg{(}\int\limits_{z_{1}}^{z_{2}}u\rho\,dz\bigg{)}_{t}+\bigg{(}\int\limits_{z_{1}}^{z_{2}}(u^{2}\rho+p)\,dz\bigg{)}_{x}\\\\[8.53581pt]
\displaystyle\quad\quad\quad\quad\quad\quad\quad\quad=M_{1}\rho_{1}u_{1}+M_{2}\rho_{2}u_{2}+z_{2x}p\big{|}_{z=z_{2}}-z_{1x}p\big{|}_{z=z_{1}}\,,\\\\[8.53581pt]
\displaystyle\bigg{(}\int\limits_{z_{1}}^{z_{2}}E\,dz\bigg{)}_{t}+\bigg{(}\int\limits_{z_{1}}^{z_{2}}(E+p)u\,dz\bigg{)}_{x}=M_{2}(E+p)\big{|}_{z=z_{2}}\\\\[8.53581pt]
\displaystyle\quad\quad\quad\quad\quad\quad\quad\quad+M_{1}(E+p)\big{|}_{z=z_{1}}-z_{2t}p\big{|}_{z=z_{2}}+z_{1t}p\big{|}_{z=z_{1}}\,,\end{array}$
(17)

where $E=\rho(u^{2}/2+gz)$ is the long-wave energy. As before, we need to
estimate the integrals equations (17).

Suppose that the stratification is weak and, following Teshukov (2007), we
represent the dimensionless density in the form

$\rho=\rho_{c}+\delta\tilde{\rho}(t,x,z),$ (18)

where $\rho_{c}=(\rho_{1}+\rho_{2})/2$ and $\delta$ is a small parameter such
that $\varepsilon^{2}\ll\delta\ll 1$. It is easy to see that, by virtue of the
first two equations of system (1), the variable $\tilde{\rho}_{z}$ satisfies
the equation

$\frac{d\tilde{\rho}_{z}}{dt}=u_{x}\tilde{\rho}_{z}-u_{z}\tilde{\rho}_{x}.$
(19)

Here $d/dt=\partial_{t}+u\partial_{x}+w\partial_{z}$ is the material
derivative. Neglecting the terms $O(\varepsilon^{2})$ in system (1) and taking
into account representation (18), we can derive the following equation for the
long-wave vorticity $u_{z}$:

$\frac{du_{z}}{dt}=\frac{\delta p_{x}}{\rho^{2}}\tilde{\rho}_{z}+\frac{\delta
g}{\rho}\tilde{\rho}_{x}.$ (20)

Equations (19) and (20) imply: if initially the variables $u_{z}$ and
$\tilde{\rho}_{z}$ are small: $u_{z}=O(\gamma)$ and
$\tilde{\rho}_{z}=O(\gamma)$ with $\varepsilon^{2}\ll\gamma\ll 1$, then for
any time

$u_{z}=O(\gamma+\delta),\quad\tilde{\rho}_{z}=O(\gamma+\delta).$ (21)

The proof is a consequence of the linearity of equations for $u_{z}$ and
$\tilde{\rho}_{z}$ (see Teshukov (2007) for details).

Formulae (21) yield the estimates

$|u-\bar{u}|=O(\gamma+\delta),\quad|\rho-\bar{\rho}|=O((\gamma+\delta)\delta).$
(22)

Indeed, for any $z$ belonging to the interval $(z_{1},z_{2})$ we have

$\begin{array}[]{l}\displaystyle|u-\bar{u}|=\bigg{|}\int\limits_{z_{1}}^{z}u_{z}\,dz-\frac{1}{\eta}\int\limits_{z_{1}}^{z_{2}}\bigg{(}\int\limits_{z_{1}}^{z^{\prime}}u_{z}\,dz\bigg{)}\,dz^{\prime}\bigg{|}\\\\[8.53581pt]
\displaystyle\quad\quad\quad\leq\max|u_{z}|\bigg{(}\int\limits_{z_{1}}^{z}\,dz+\frac{1}{\eta}\int\limits_{z_{1}}^{z_{2}}(z^{\prime}-z_{1})\,dz^{\prime}\bigg{)}\leq\frac{3\eta}{2}\max|u_{z}|.\end{array}$

Taking into account representation (18), we similarly obtain the second
estimate in (22).

Using the identities $\rho=\bar{\rho}+(\rho-\bar{\rho})$,
$u=\bar{u}+(u-\bar{u})$ and formulae (22), we obtain the following asymptotic
estimates of integrals in (17):

$\begin{array}[]{l}\displaystyle\int\limits_{z_{1}}^{z_{2}}\rho
u\,dz=\bar{u}\bar{\rho}\eta+O((\gamma+\delta)^{2}\delta),\quad\int\limits_{z_{1}}^{z_{2}}\rho
u^{2}\,dz=(\bar{u}^{2}+q^{2})\bar{\rho}\eta+O((\gamma+\delta)^{2}\delta),\\\\[8.53581pt]
\displaystyle\int\limits_{z_{1}}^{z_{2}}\rho
u^{3}\,dz=(\bar{u}^{2}+3q^{2})\bar{u}\bar{\rho}\eta+O((\gamma+\delta)^{2}\delta),\quad\int\limits_{z_{1}}^{z_{2}}z\rho\,dz=\Big{(}Z+h_{1}+\frac{\eta}{2}\Big{)}\bar{\rho}\eta+O((\gamma+\delta)\delta),\end{array}$

with the classical definition of $q^{2}$:

$q^{2}=\frac{1}{\eta}\int\limits_{z_{1}}^{z_{2}}(u-\bar{u})^{2}dz.$ (23)

In view of (6) and (22) the pressure in this layer is

$p=p^{*}+\rho_{2}gh_{2}+g\bar{\rho}(z_{2}-z)+O((\gamma+\delta)\delta)$

and, consequently,

$\int\limits_{z_{1}}^{z_{2}}(p+g\rho
z)u\,dz=(p^{*}+g\rho_{2}h_{2}+g\,\bar{\rho}(Z+h_{1}+\eta))\bar{u}\eta+O((\gamma+\delta)\delta).$

Note that these estimates are similar to those obtained in Gavrilyuk et al.
(2019). Neglecting the terms of order $O((\gamma+\delta)\delta)$ in the
previous integrals, we present system (17) in the form:

$\begin{array}[]{l}\displaystyle\eta_{t}+(\bar{u}\eta)_{x}=M_{1}+M_{2},\quad(\bar{\rho}\eta)_{t}+(\bar{u}\bar{\rho}\eta)_{x}=M_{1}\rho_{1}+M_{2}\rho_{2},\\\\[8.53581pt]
\displaystyle(\bar{u}\bar{\rho}\eta)_{t}+\Big{(}(\bar{u}^{2}+q^{2})\bar{\rho}\eta+(p^{*}+g\rho_{2}h_{2})\eta+\frac{g\bar{\rho}\eta^{2}}{2}\Big{)}_{x}\\\\[8.53581pt]
\displaystyle\quad\quad\quad=M_{1}\rho_{1}u_{1}+M_{2}\rho_{2}u_{2}+(p^{*}+g\rho_{2}h_{2})\eta_{x}-g\bar{\rho}\eta(Z+h_{1})_{x},\\\\[8.53581pt]
\displaystyle\bigg{(}\frac{(\bar{u}^{2}+q^{2})}{2}\bar{\rho}\eta+g(Z+h_{1}+\eta)\bar{\rho}\eta-\frac{g\bar{\rho}\eta^{2}}{2}\bigg{)}_{t}+\bigg{(}\frac{(\bar{u}^{2}+3q^{2})}{2}\bar{u}\bar{\rho}\eta\\\\[11.38109pt]
\displaystyle\quad+(p^{*}+g\rho_{2}h_{2})\bar{u}\eta+g(Z+h_{1}+\eta)\bar{u}\bar{\rho}\eta\bigg{)}_{x}=M_{1}\rho_{1}\Big{(}\frac{u_{1}^{2}}{2}+g(Z+h_{1})\Big{)}\\\\[8.53581pt]
\displaystyle\quad\quad+M_{2}\rho_{2}\Big{(}\frac{u_{2}^{2}}{2}+g(Z+h_{1}+\eta)\Big{)}+(p^{*}+g\rho_{2}h_{2})(M_{1}+M_{2}-\eta_{t})\\\\[8.53581pt]
\displaystyle\quad\quad\quad+g(M_{1}+Z_{t}+h_{1t})\bar{\rho}\eta-d.\end{array}$
(24)

To account for the energy dissipation in the last equation of (24) an extra
term $d$ has been added (formula for $d$ will be proposed later).

###  2.3 Differential consequence of the energy equation

Let us derive an equation for the variable $q$ measuring the distortion of the
velocity profile. For this, we first note that the first three equations (24)
imply:

$\begin{array}[]{l}\displaystyle\bar{\rho}_{t}+\bar{u}\bar{\rho}_{x}=\frac{1}{\eta}\Big{(}(\rho_{1}-\bar{\rho})M_{1}+(\rho_{2}-\bar{\rho})M_{2}\Big{)},\\\\[8.53581pt]
\displaystyle\bar{u}_{t}+\bar{u}\bar{u}_{x}+2qq_{x}+gh_{1x}+\Big{(}g+\frac{q^{2}}{\eta}\Big{)}\eta_{x}+\frac{g\rho_{2}}{\bar{\rho}}h_{2x}+\frac{1}{\bar{\rho}}\Big{(}q^{2}+\frac{g\eta}{2}\Big{)}\bar{\rho}_{x}+\frac{1}{\bar{\rho}}p^{*}_{x}\\\\[8.53581pt]
\displaystyle\quad\quad\quad=-gZ_{x}+\frac{1}{\bar{\rho}\eta}\Big{(}(u_{1}-\bar{u})\rho_{1}M_{1}+(u_{2}-\bar{u})\rho_{2}M_{2}\Big{)}.\end{array}$
(25)

The non-conservative form of the last equation in (24) is

$\begin{array}[]{l}\displaystyle\Big{(}\frac{\bar{u}^{2}+q^{2}}{2}+g(Z+h_{1}+\eta)\Big{)}\big{(}(\bar{\rho}\eta)_{t}+(\bar{u}\bar{\rho}\eta)_{x}\big{)}+q\bar{\rho}\eta\big{(}q_{t}+(\bar{u}q)_{x}\big{)}\\\\[8.53581pt]
\displaystyle\quad+\bar{u}\bar{\rho}\eta\Big{(}\bar{u}_{t}+\bar{u}\bar{u}_{x}+2qq_{x}+gh_{1x}+\Big{(}g+\frac{q^{2}}{\eta}\Big{)}\eta_{x}+\frac{g\rho_{2}}{\bar{\rho}}h_{2x}+\frac{1}{\bar{\rho}}\Big{(}q^{2}+\frac{g\eta}{2}\Big{)}\bar{\rho}_{x}+\frac{1}{\bar{\rho}}p^{*}_{x}\Big{)}\\\\[8.53581pt]
\displaystyle\quad\quad-\frac{g\eta^{2}}{2}\big{(}\bar{\rho}_{t}+\bar{u}\bar{\rho}_{x}\big{)}+g\bar{u}\bar{\rho}\eta
Z_{x}=M_{1}\rho_{1}\Big{(}\frac{u_{1}^{2}}{2}+g(Z+h_{1})\Big{)}\\\\[8.53581pt]
\displaystyle\quad\quad\quad+M_{2}\rho_{2}\Big{(}\frac{u_{2}^{2}}{2}+g(Z+h_{1}+\eta)\Big{)}+M_{1}g\bar{\rho}\eta-d.\end{array}$

Taking into account relations (25) and the second equation in (24), we obtain

$\begin{array}[]{l}\displaystyle
q\bar{\rho}\eta\big{(}q_{t}+(\bar{u}q)_{x}\big{)}=\frac{M_{1}\rho_{1}}{2}\big{(}(u_{1}-\bar{u})^{2}-q^{2}-g\eta\big{)}\\\\[8.53581pt]
\displaystyle\quad\quad\quad+\frac{M_{2}\rho_{2}}{2}\big{(}(u_{2}-\bar{u})^{2}-q^{2}+g\eta\big{)}+\frac{g\bar{\rho}\eta}{2}(M_{1}-M_{2})-d.\end{array}$
(26)

From this consequence and the mass equation for $\eta$ one can derive the
transport equation for $q/\eta$ which can be interpreted as the evolution
equation for the flow vorticity in the intermediate layer. A priori, the
vorticity can change its sign during the flow evolution.

###  2.4 Final three-layer system

The final system can be written as

$\begin{array}[]{l}\displaystyle
h_{1t}+(u_{1}h_{1})_{x}=-M_{1},\quad\eta_{t}+(\bar{u}\eta)_{x}=M_{1}+M_{2},\quad
h_{2t}+(u_{2}h_{2})_{x}=-M_{2},\\\\[8.53581pt] \displaystyle
u_{1t}+\bigg{(}\frac{u_{1}^{2}}{2}+gh_{1}+g\eta\frac{\bar{\rho}}{\rho_{1}}+gh_{2}\frac{\rho_{2}}{\rho_{1}}+\frac{p^{*}}{\rho_{1}}\bigg{)}_{x}=-gZ_{x},\\\\[11.38109pt]
\displaystyle
u_{2t}+\bigg{(}\frac{u_{2}^{2}}{2}+\frac{p^{*}}{\rho_{2}}\bigg{)}_{x}=0,\quad
Q_{t}+\bigg{(}\rho_{1}u_{1}^{2}h_{1}+(\bar{u}^{2}+q^{2})\bar{\rho}\eta+\rho_{2}u_{2}^{2}h_{2}+Hp^{*}\\\\[8.53581pt]
\displaystyle\quad+\frac{gh_{1}^{2}\rho_{1}}{2}+\frac{g\eta^{2}\bar{\rho}}{2}+\frac{gh_{2}^{2}\rho_{2}}{2}+g\bar{\rho}h_{1}\eta+g\rho_{2}h_{2}(h_{1}+\eta)\bigg{)}_{x}\\\\[8.53581pt]
\displaystyle\quad\quad=-(gh_{1}\rho_{1}+g\eta\bar{\rho}+gh_{2}\rho_{2}+p^{*})Z_{x},(\bar{\rho}\eta)_{t}+(\bar{u}\bar{\rho}\eta)_{x}=M_{1}\rho_{1}+M_{2}\rho_{2},\\\\[8.53581pt]
\displaystyle
q_{t}+(\bar{u}q)_{x}=\frac{M_{1}\rho_{1}}{2q\eta\bar{\rho}}\big{(}(u_{1}-\bar{u})^{2}-q^{2}-g\eta\big{)}\\\\[8.53581pt]
\displaystyle\quad\quad+\frac{M_{2}\rho_{2}}{2q\eta\bar{\rho}}\big{(}(u_{2}-\bar{u})^{2}-q^{2}+g\eta\big{)}+\frac{g}{2q}(M_{1}-M_{2})-\frac{d}{q\bar{\rho}\eta},\end{array}$
(27)

where $Q=\rho_{1}u_{1}h_{1}+\bar{\rho}\bar{u}\eta+\rho_{2}u_{2}h_{2}$ is the
total mass discharge, and $H=H_{0}-Z$ is the total fluid thickness. To
represent the depth-averaged equations in the form of balance laws, we
combined the momentum equations in the outer layers and intermediate layer
into the total momentum equation for $Q$.

Instead of the cumbersome total energy equation, we use its differential
consequence for the $q$ variable. This replacement of the balance law does not
affect the construction of solutions in the class of continuous flows, but in
the event of discontinuities, the solutions may differ. However, in the case
of small jump amplitude this difference is negligible. An example of such an
approach is shown by Lipatov et al. (2021) for the case of compressible flows.
This difference manifests itself in a fairly small region, outside of which
the solutions are close or even almost coincide. The appropriateness of using
the equation for the velocity distortion $q$ is also underlined by Gavrilyuk
et al. (2016, 2019) where free surface two-layer flows taking into account
mixing were considered.

To close the model it is necessary to determine the mass entrainment terms
$M_{1}$, $M_{2}$, and the energy dissipation term $d$. Since stratification is
considered weak, we suppose that the entrainment of fluid from the outer
layers into the intermediate one is symmetric. Following Gavrilyuk et al.
(2016, 2019); Liapidevskii et al. (2018), we take the entrainment velocities
$M_{1}$ and $M_{2}$ in the form

$M_{1}=M_{2}=\sigma q,\quad\sigma={\rm const}>0.$ (28)

The interpretation of $\sigma$ parameter comes from the theory of plane
turbulent shear. More, exactly, $2\sigma$ is the ratio of the shear stress to
the turbulent energy $k_{T}$, $k_{T}\approx q^{2}/2$. This ratio is
approximately constant and equal to $0.3$ (see Pope (2000), Table 5.4, p.
157). Thus, in the following, we always take $\sigma=0.15$. A detailed
justification of the entrainment terms can be found in Gavrilyuk et al.
(2016).

To account for the energy transfer from large scale eddies to small scale
eddies we added the extra term $d$ in the last equation of (24):

$d=\frac{\sigma\kappa}{2}\bar{\rho}|q|^{3},\quad\kappa={\rm const}>0.$ (29)

Such a dissipation term is classical in the theory of turbulent flows of
homogeneous fluids. Indeed, the rate of turbulent energy dissipation is
usually written as (see Pope (2000), p. 244):

$\epsilon=\frac{k_{T}^{3/2}}{L_{11}}\left(\frac{L_{11}}{L}\right),$ (30)

where $L_{11}$ is the length scale of the energy-containing eddies, and
$\varkappa=L_{11}/L$ is a dimensionless parameter which tends asymptotically
for high Reynolds numbers to a constant value (see Pope (2000), p. 244–245).
In our case, the horizontal velocity in the intermediate shear layer can be
approximated by a linear profile :
$u=\bar{u}+\omega\left(z-(z_{1}+z_{2})/2\right)$. Here $\omega={\rm const}$ is
the component of the horizontal vorticity. Since the velocity profile is
symmetric with respect to $\bar{u}$, one can conclude that the size of the
energy-containing eddies is $\eta/2$. Hence,

$\epsilon=\frac{k_{T}^{3/2}}{L_{11}}\left(\frac{L_{11}}{L}\right)=\frac{|q|^{3}}{2\sqrt{2}L_{11}}\left(\frac{L_{11}}{L}\right)=\frac{|q|^{3}}{\eta\sqrt{2}}\left(\frac{L_{11}}{L}\right)=\frac{|q|^{3}}{2\eta}\left(\frac{\sqrt{2}L_{11}}{L}\right)$
(31)

Thus

$\frac{L_{11}}{L}=\frac{\kappa\sigma}{\sqrt{2}}.$ (32)

The experimentally observed data for the ratio $L_{11}/L$ is (Pope (2000), p.
245):

$0.8>\frac{L_{11}}{L}=\frac{\kappa\sigma}{\sqrt{2}}>0.43.$ (33)

For $\sigma\approx 0.15$ it implies the following estimation for $\kappa$ :

$\simeq 8>\kappa>\simeq 4.$ (34)

Finally, $\kappa$ is the only phenomenological parameter in our model. As we
have observed in our numerical experiments, the results do not depend too much
on the specific choice of $\kappa\in(4,8)$.

Thus, equations (27) with the additional relation $h_{1}+\eta+h_{2}=H_{0}-Z$
form a closed system for nine unknown functions $h_{1}$, $\eta$, $h_{2}$,
$u_{1}$, $\bar{u}$, $u_{2}$, $\bar{\rho}$, $q$ and $p^{*}$.

##  3 Governing equations in the Boussinesq approximation

Let us return back to dimensional variables. The depth-averaged equations will
not change if we drop the ‘hats’ on dimensional variables. Further we apply
the Boussinesq approximation assuming that the ratio
$(\rho_{1}-\rho_{2})/\rho_{2}\ll 1$ is negligible, but not the buoyancy terms

$b=g\frac{\rho_{1}-\rho_{2}}{\rho_{2}},\quad\bar{b}=g\frac{\bar{\rho}-\rho_{2}}{\rho_{2}}\,.$

Since the fourth, sixth, and seventh equations in (27) admit the
representation

$\begin{array}[]{l}\displaystyle
u_{1t}+\bigg{(}\frac{u_{1}^{2}}{2}+\frac{p^{*}}{\rho_{1}}+gh_{1}\frac{\rho_{1}-\rho_{2}}{\rho_{1}}+g\eta\frac{\bar{\rho}-\rho_{2}}{\rho_{1}}\bigg{)}_{x}=-gZ_{x}\frac{\rho_{1}-\rho_{2}}{\rho_{1}},\\\\[8.53581pt]
\displaystyle
Q_{t}+\bigg{(}\rho_{1}u_{1}^{2}h_{1}+(\bar{u}^{2}+q^{2})\bar{\rho}\eta+\rho_{2}u_{2}^{2}h_{2}+\frac{gh_{1}^{2}}{2}(\rho_{1}-\rho_{2})+gh_{1}\eta(\bar{\rho}-\rho_{2})\\\\[8.53581pt]
\displaystyle\quad\quad\quad+\frac{g\eta^{2}}{2}(\bar{\rho}-\rho_{2})+Hp^{*}\bigg{)}_{x}=-(gh_{1}(\rho_{1}-\rho_{2})+g\eta(\bar{\rho}-\rho_{2})+p^{*})Z_{x},\\\\[8.53581pt]
\displaystyle\big{(}(\bar{\rho}-\rho_{2})\eta\big{)}_{t}+\big{(}(\bar{\rho}-\rho_{2})\bar{u}\eta\big{)}_{x}=(\rho_{1}-\rho_{2})M_{1}\end{array}$

and the ratios $1/\rho_{1}$ and $1/\bar{\rho}$ can be approximated as

$\frac{1}{\rho_{1}}=\frac{1}{\rho_{2}+(\rho_{1}-\rho_{2})}\approx\frac{1}{\rho_{2}}-\frac{\rho_{1}-\rho_{2}}{\rho_{2}^{2}}\,,\quad\frac{1}{\rho_{1}}\approx\frac{1}{\rho_{2}}-\frac{\bar{\rho}-\rho_{2}}{\rho_{2}^{2}},$

the governing equations take the form

$\begin{array}[]{l}\displaystyle h_{1t}+(u_{1}h_{1})_{x}=-\sigma
q,\quad\eta_{t}+(\bar{u}\eta)_{x}=2\sigma q,\quad
h_{2t}+(u_{2}h_{2})_{x}=-\sigma q,\\\\[8.53581pt] \displaystyle
u_{1t}+\bigg{(}\frac{u_{1}^{2}}{2}+\frac{p^{*}}{\rho_{2}}+bh_{1}+\bar{b}\eta\bigg{)}_{x}=-bZ_{x},\quad
u_{2t}+\bigg{(}\frac{u_{2}^{2}}{2}+\frac{p^{*}}{\rho_{2}}\bigg{)}_{x}=0,\\\\[11.38109pt]
\displaystyle\bar{Q}_{t}+\bigg{(}u_{1}^{2}h_{1}+(\bar{u}^{2}+q^{2})\eta+u_{2}^{2}h_{2}+\frac{bh_{1}^{2}}{2}+\bar{b}h_{1}\eta+\frac{\bar{b}\eta^{2}}{2}+\frac{Hp^{*}}{\rho_{2}}\bigg{)}_{x}\\\\[11.38109pt]
\displaystyle\quad\quad=-\bigg{(}bh_{1}+\bar{b}\eta+\frac{p^{*}}{\rho_{2}}\bigg{)}Z_{x},\quad(\bar{b}\eta)_{t}+(\bar{u}\bar{b}\eta)_{x}=\sigma
qb,\quad q_{t}+(\bar{u}q)_{x}=\varphi.\end{array}$ (35)

Here we have already taken into account formulae (28), (29) and used the
notation

$\varphi=\frac{\sigma}{2\eta}\Big{(}(u_{1}-\bar{u})^{2}+(u_{2}-\bar{u})^{2}-(2+\kappa\,{\rm
sign}\,q)q^{2}-b\eta\Big{)}$

and $\bar{Q}=u_{1}h_{1}+\bar{u}\eta+u_{2}h_{2}$. In view of the first three
equations (35) we have $\bar{Q}_{x}=-H_{t}=Z_{t}$. Therefore, the variable
$\bar{Q}$ is known (up to an arbitrary function of time).

One can derive from (24), (13) and (11) the conservation equation of the total
energy of a three-layer flow in the Boussinesq approximation:

$\begin{array}[]{l}\displaystyle\bigg{(}\frac{u_{1}^{2}h_{1}}{2}+\frac{(\bar{u}^{2}+q^{2})\eta}{2}+\frac{u_{2}^{2}h_{2}}{2}+\frac{bh_{1}^{2}}{2}+\bar{b}\eta\Big{(}h_{1}+\frac{\eta}{2}\Big{)}\bigg{)}_{t}+\bigg{(}\frac{u_{1}^{3}h_{1}}{2}\\\\[8.53581pt]
\displaystyle\quad+\frac{(\bar{u}^{2}+3q^{2})\bar{u}\eta}{2}+\frac{u_{2}^{3}h_{2}}{2}+(bh_{1}+\bar{b}\eta)u_{1}h_{1}+\bar{u}\bar{b}\eta(h_{1}+\eta)+p^{*}\bar{Q}\bigg{)}_{x}\\\\[8.53581pt]
\displaystyle\quad\quad=p^{*}Z_{t}-(u_{1}h_{1}b+v\eta\bar{b})Z_{x}-\frac{\sigma\kappa}{2}|q|^{3}\,.\end{array}$
(36)

As mentioned above, we have replaced this cumbersome equation by its more
convenient differential consequence for the variable $q$. A priori, this
procedure is equivalent for smooth solutions, but not for discontinuous
solutions. In the next section, we will show that such a change in the
conservation law does not significantly affect the structure of discontinuous
solutions.

It is convenient to derive the following consequences of system (35) for
variables $\bar{b}$, $s=q/\eta$ and $\bar{u}$:

$\begin{array}[]{l}\displaystyle\bar{b}_{t}+\bar{u}\bar{b}_{x}=\frac{\sigma
q}{\eta}(b-2\bar{b}),\quad
s_{t}+\bar{u}s_{x}=\frac{1}{\eta}\Big{(}\varphi-\frac{2\sigma
q^{2}}{\eta}\Big{)},\\\\[8.53581pt]
\displaystyle\bar{u}_{t}+\bar{u}\bar{u}_{x}+2qq_{x}+\bar{b}h_{1x}+\Big{(}\bar{b}+\frac{q^{2}}{\eta}\Big{)}\eta_{x}+\frac{\eta}{2}\bar{b}_{x}+\frac{1}{\rho_{2}}p^{*}_{x}\\\\[8.53581pt]
\displaystyle\quad\quad\quad=\frac{\sigma
q}{\eta}(u_{1}+u_{2}-2\bar{u})-\bar{b}Z_{x}.\end{array}$ (37)

Note that these equations can be obtained from (25) and (26) by passing to the
Boussinesq approximation and using formulae (28) and (29). We write now the
first equation of (37):

$\bar{b}_{t}+\bar{u}\bar{b}_{x}=\frac{\sigma q}{\eta}(b-2\bar{b}).$ (38)

The buoyancy $\bar{b}$ has the following property: if, initially,
$\bar{b}>b/2$, then this property is valid any time. Moreover, $\bar{b}$ is
always positive. These properties follow directly from (38). In particular,
this implies that for any time $\rho_{2}<\bar{\rho}<\rho_{1}$, i.e. the
stratification stays always stable.

Let us eliminate the pressure $p^{*}$ on the upper lid. For this, we introduce
new variables

$r=u_{1}-u_{2},\quad R=\bar{Q}-u_{2}H=(u_{1}-u_{2})h_{1}+(\bar{u}-u_{2})\eta$

and note that the variable $u_{2}H$ satisfies equation

$(u_{2}H)_{t}+\bigg{(}\frac{u_{2}^{2}H}{2}+\frac{Hp^{*}}{\rho_{2}}\bigg{)}_{x}=-\Big{(}Z_{t}+\frac{u_{2}}{2}Z_{x}\Big{)}u_{2}-\frac{p^{*}}{\rho_{2}}Z_{x}.$

Subtracting this balance law from the sixth equation (35), as well as the
fifth from the fourth, we get an evolutionary system.

Thus, to define six unknowns $(h_{1},h_{2},r,R,\bar{b},q)$ we obtain the
closed system of balance laws

$\begin{array}[]{l}\displaystyle h_{1t}+(u_{1}h_{1})_{x}=-\sigma
q,\quad\eta_{t}+(\bar{u}\eta)_{x}=2\sigma q,\\\\[8.53581pt]
\displaystyle(\bar{b}\eta\big{)}_{t}+(\bar{u}\bar{b}\eta)_{x}=\sigma qb,\quad
r_{t}+\bigg{(}\frac{u_{1}^{2}-u_{2}^{2}}{2}+bh_{1}+\bar{b}\eta\bigg{)}_{x}=-bZ_{x},\\\\[8.53581pt]
\displaystyle
R_{t}+\bigg{(}u_{1}^{2}h_{1}+(\bar{u}^{2}+q^{2})\eta+\Big{(}h_{2}-\frac{H}{2}\Big{)}u_{2}^{2}+\frac{bh_{1}^{2}}{2}+\bar{b}h_{1}\eta+\frac{\bar{b}\eta^{2}}{2}\bigg{)}_{x}\\\\[8.53581pt]
\displaystyle\quad\quad\quad=u_{2}Z_{t}+\Big{(}\frac{u_{2}^{2}}{2}-bh_{1}-\bar{b}\eta\Big{)}Z_{x},\quad
q_{t}+(\bar{u}q)_{x}=\varphi,\end{array}$ (39)

where the variables $h_{2}$, $u_{1}$, $u_{2}$ and $\bar{u}$ can be expressed
as

$h_{2}=H-h_{1}-\eta,\quad u_{2}=\frac{\bar{Q}-R}{H},\quad
u_{1}=r+u_{2},\quad\bar{u}=\frac{R-rh_{1}}{\eta}+u_{2}.$

System (39) describes non-stationary three-layer hydrostatic flows with mixing
in the Boussinesq approximation.

###  3.1 Characteristics of equations (39)

Let us rewrite system (39) in the form

$\mathbf{U}_{t}+\mathbf{A}\mathbf{U}_{x}=\mathbf{F},$ (40)

where $\mathbf{U}=(\bar{b},s,h_{1},\eta,r,R)^{\rm T}$ is the vector of
unknowns, $\mathbf{F}$ is the right-hand which doesn’t contain derivatives of
$\mathbf{U}$, and

$\mathbf{A}=\begin{pmatrix}\bar{u}&0&0&0&0&0\\\ 0&\bar{u}&0&0&0&0\\\
0&0&u_{1}&0&h_{1}&-h_{1}/H\\\ 0&0&-r&u_{2}&-h_{1}&1-\eta/H\\\
\eta&0&b&\bar{b}&u_{1}&-r/H\\\ a_{1}&a_{2}&a_{3}&a_{4}&a_{5}&a_{6}\\\
\end{pmatrix}$

is the $6\times 6$ matrix. The last row of the matrix is

$\begin{array}[]{l}\displaystyle
a_{1}=\Big{(}h_{1}+\frac{\eta}{2}\Big{)}\eta,\quad a_{2}=2\eta^{2}q,\quad
a_{3}=(u_{1}-\bar{u})^{2}-(u_{2}-\bar{u})^{2}+bh_{1}+\bar{b}\eta,\\\\[5.69054pt]
\displaystyle a_{4}=3q^{2}-(u_{2}-\bar{u})^{2}+(h_{1}+\eta)\bar{b},\quad
a_{5}=2(u_{1}-\bar{u})h_{1},\quad
a_{6}=u_{2}+2\Big{(}\bar{u}-\frac{\bar{Q}}{H}\Big{)}.\end{array}$

The eigenvalues of $\mathbf{A}(\mathbf{U})$ are the characteristic velocities
of system (40). The sixth order polynomial equation ${\rm
det}(\mathbf{A}-\lambda\mathbf{I})=0$ has the root $\lambda=\bar{u}$ of
multiplicity two. To determine the remaining four roots, one has to solve the
equation

$\begin{array}[]{l}\displaystyle\chi(\lambda)=\big{(}(u_{1}-\lambda)^{2}-(b-\bar{b})h_{1}\big{)}\big{(}(u_{2}-\lambda)^{2}-\bar{b}h_{2}\big{)}\eta\\\\[5.69054pt]
\displaystyle\quad+\big{(}(\bar{u}-\lambda)^{2}-3q^{2}\big{)}\Big{(}\big{(}(u_{1}-\lambda)^{2}-(b-\bar{b})h_{1}\big{)}h_{2}+\big{(}(u_{2}-\lambda)^{2}-\bar{b}h_{2}\big{)}h_{1}\Big{)}=0.\end{array}$
(41)

Consider the characteristic polynomial $\chi(\lambda)$ in the case of equal
velocities $u_{1}=u_{2}=\overline{u}$. Due to the Galilean invariance of the
governing equations in the case $\mathbf{F}=\mathbf{0}$, one can always take
the depth averaged velocities vanishing. However, $q$ does not vanish. Then
equation (41) multiplied by $(h_{1}\eta h_{2})^{-1}$ is

$\begin{array}[]{l}\displaystyle\bigg{(}\frac{1}{h_{1}h_{2}}+\frac{1}{h_{1}\eta}+\frac{1}{\eta
h_{2}}\bigg{)}\lambda^{4}-\bigg{(}\frac{\bar{b}}{h_{1}}+\frac{b-\bar{b}}{h_{2}}+\frac{b}{\eta}+\frac{3q^{2}}{\eta}\Big{(}\frac{1}{h_{1}}+\frac{1}{h_{2}}\Big{)}\bigg{)}\lambda^{2}\\\\[11.38109pt]
\displaystyle\quad\quad\quad\quad\quad+\frac{3q^{2}b}{\eta}+\bar{b}(b-\bar{b})=0.\end{array}$
(42)

Let $\mu_{1}=\lambda_{1}^{2}$ and $\mu_{2}=\lambda_{2}^{2}$ be the roots of
the corresponding bi-quadratic equation for $\mu=\lambda^{2}$. Since the
stratification is stable, one has $b>\overline{b}>0$. Then one obviously has

$\mu_{1}+\mu_{2}>0,\quad\mu_{1}\mu_{2}>0.$

Let us first show that for $q=0$ all roots $\mu_{i}$ are positive. Indeed, for
this it is sufficient to show that the corresponding discriminant $\Delta$ is
positive :

$\begin{array}[]{l}\displaystyle\Delta=\left(\frac{\bar{b}^{2}}{h_{1}^{2}}+\frac{(b-\bar{b})^{2}}{h_{2}^{2}}+\frac{b^{2}}{\eta^{2}}+\frac{2{\bar{b}}(b-\bar{b})}{h_{1}h_{2}}+\frac{2b(b-\bar{b})}{\eta
h_{2}}+\frac{2b\bar{b}}{h_{1}\eta}\right)\\\\[11.38109pt]
\displaystyle\quad\quad\quad-4\bar{b}(b-\bar{b})\left(\frac{1}{h_{1}h_{2}}+\frac{1}{h_{1}\eta}+\frac{1}{\eta
h_{2}}\right).\end{array}$

Let $\mathbf{z}=(1/h_{1},1/h_{2},1/\eta)^{\rm T}$. Then one can write

$\Delta=\mathbf{z}^{\rm
T}\cdot\mathbf{B}\cdot\mathbf{z},\quad{\mathbf{B}}=\left(\begin{array}[]{ccc}{\bar{b}}^{2}&-{\bar{b}}(b-{\bar{b}})&-{\bar{b}}(b-2{\bar{b}})\\\
-{\bar{b}}(b-{\bar{b}})&(b-{\bar{b}})^{2}&(b-{\bar{b}})(b-2{\bar{b})}\\\
-{\bar{b}}(b-2{\bar{b}})&(b-{\bar{b}})(b-2{\bar{b})}&b^{2}\\\
\end{array}\right).$

Both the second and third principal minors of $\mathbf{B}$ are identically
zero, so the matrix is only non-negative definite: $\mathbf{B}\geq 0$. The
case $q\neq 0$ improves the situation: the corresponding matrix is always
degenerate, but first and second principal minor are positive. Overall, the
roots $\mu_{i}$ are positive, but they can coincide. This means that in the
absence of a velocity shear, all eigenvalues are real.

When the governing equations are hyperbolic the concepts of supercritical and
subcritical flow can be introduced. We say that the flow is supercritical if
all roots $\lambda=\lambda_{i}$ of the characteristic equation (41) are
positive, and subcritical if there is at least one negative root.

##  4 Stationary solutions

Stationary solutions of system (35) (or (39)) are determined by the equations

$\begin{array}[]{l}\displaystyle(u_{1}h_{1})^{\prime}=-\sigma
q,\quad(\bar{u}\eta)^{\prime}=2\sigma
q,\quad(\bar{u}\eta\bar{b})^{\prime}=\sigma qb,\\\\[8.53581pt]
\displaystyle(\bar{u}q)^{\prime}=\varphi,\quad\Big{(}\frac{u_{1}^{2}-u_{2}^{2}}{2}+bh_{1}+\bar{b}\eta\Big{)}^{\prime}=-bZ^{\prime},\\\\[8.53581pt]
\displaystyle\bar{u}\bar{u}^{\prime}+2qq^{\prime}+\bar{b}h_{1}^{\prime}+\Big{(}\bar{b}+\frac{q^{2}}{\eta}\Big{)}\eta^{\prime}+\frac{\eta}{2}\bar{b}^{\prime}-u_{2}u_{2}^{\prime}=\frac{\sigma
q}{\eta}(u_{1}-2\bar{u}+u_{2})-\bar{b}Z^{\prime}.\end{array}$ (43)

Here ‘prime’ means the derivative with respect to $x$. Since
$u_{2}=(\bar{Q}-u_{1}h_{1}-\bar{u}\eta)/h_{2}$ and $h_{2}=H_{0}-Z-h_{1}-\eta$,
equations (43) form a closed system for unknowns $h_{1}$, $\eta$, $u_{1}$,
$\bar{u}$, $q$ and $\bar{b}$. Let us rewrite this system in normal form

$\begin{array}[]{l}\displaystyle u_{1}^{\prime}=-\frac{\sigma
q+u_{1}h_{1}^{\prime}}{h_{1}},\quad\bar{u}^{\prime}=\frac{2\sigma
q-\bar{u}\eta^{\prime}}{\eta},\quad\quad\bar{b}^{\prime}=\frac{\sigma
q}{\bar{u}\eta}(b-2\bar{b}),\\\\[11.38109pt] \displaystyle
q^{\prime}=\frac{\varphi-q\bar{u}^{\prime}}{\bar{u}},\quad\eta^{\prime}=\frac{C_{1}-A_{1}h_{1}^{\prime}}{A_{2}},\quad
h_{1}^{\prime}=\frac{B_{2}C_{1}-A_{2}C_{2}}{A_{1}B_{2}-A_{2}^{2}}\,,\end{array}$
(44)

where

$\begin{array}[]{l}\displaystyle
A_{1}=\frac{u_{1}^{2}}{h_{1}}+\frac{u_{2}^{2}}{h_{2}}-b,\quad
A_{2}=\frac{u_{2}^{2}}{h_{2}}-\bar{b},\quad C_{1}=\eta\bar{b}^{\prime}-\sigma
q\Big{(}\frac{u_{1}}{h_{1}}-\frac{u_{2}}{h_{2}}\Big{)}-\Big{(}\frac{u_{2}^{2}}{h_{2}}-b\Big{)}Z^{\prime},\\\\[8.53581pt]
\displaystyle B_{2}=\frac{\bar{u}^{2}-3q^{2}}{\eta}+A_{2},\quad
C_{2}=\frac{\eta\bar{b}^{\prime}}{2}-\frac{\sigma
q}{\eta}\Big{(}u_{1}-4\bar{u}+\frac{h_{2}-\eta}{h_{2}}u_{2}\Big{)}+\frac{2q}{\bar{u}}\Big{(}\varphi-\frac{2\sigma
q^{2}}{\eta}\Big{)}-A_{2}Z^{\prime}.\end{array}$

We note that

$\chi(0)=(A_{1}B_{2}-A_{2}^{2})h_{1}\eta h_{2}$

and, consequently, if the denominator of the right-hand side in the last
equation (44) vanishes, then $\lambda=0$ is the root of characteristic
equation (41). We say that a stationary flow is supercritical if $\chi(0)>0$,
and subcritical if $\chi(0)<0$.

To solve ODEs (44) numerically, we use the standard ode45 procedure of the
MATLAB package that implements the fourth-order Runge–Kutta method.

###  4.1 Mixing layer formation

To construct a stationary solution to the mixing layer problem, it is
necessary to determine the values of $\bar{u}$, $q$ and $\bar{b}$ as $\eta\to
0$. Without loss of generality, we assume that $\eta=0$ at $x=0$, and the
values of the functions at this point are marked with the subscript ‘zero’.
Velocities $u_{10}>0$, $u_{20}>0$ ($u_{10}\neq u_{20}$) and thicknesses
$h_{10}>0$, $h_{20}>0$ for the outer layers are assumed to be specified. We
also assume that there are finite limits of all functions and their
derivatives. By virtue of the second, third, fourth, and sixth equations of
system (43) for $\eta\to 0$, we obtain

$\begin{array}[]{l}\displaystyle\eta^{\prime}\to\frac{2\sigma
q_{0}}{\bar{u}_{0}},\quad\bar{b}_{0}=\frac{b}{2},\quad(u_{10}-\bar{u}_{0})^{2}+(u_{20}-\bar{u}_{0})^{2}=(2+\kappa)q_{0}^{2},\\\\[8.53581pt]
\displaystyle\eta^{\prime}\to\frac{\sigma}{q_{0}}(u_{10}+u_{20}-2\bar{u}_{0}).\end{array}$

This means that

$q_{0}=\sqrt{(u_{10}+u_{20}-2\bar{u}_{0})\bar{u}_{0}/2}$ (45)

($q_{0}>0$ since during the mixing layer formation, it expands) and
$\bar{u}_{0}$ is the root of the quadric equation

$\bar{u}_{0}^{2}-\frac{(6+\kappa)(u_{10}+u_{20})}{2(4+\kappa)}\bar{u}_{0}+\frac{u_{10}^{2}+u_{20}^{2}}{4+\kappa}=0.$
(46)

Variation of the dissipation parameter $\kappa$ allows us to obtain the unique
root of this equation $\bar{u}_{0}>(u_{10}+u_{20})/2$ lying between the values
of $u_{10}$ and $u_{20}$. This can always be achieved for not too large
absolute value of the relative velocity $|u_{20}-u_{10}|$.

These restrictions on the variables $\bar{u}_{0}$ and $q_{0}$ apply only to
the initial stage of mixing, where the thickness of the interlayer is close to
zero. If non-stationary equations (39) are used to calculate the evolution of
the mixing layer, then instead of conditions (45), (46) we can set, for
example $\bar{u}_{0}=(u_{1}+u_{2})/2$, $q_{0}=0$. It will be shown below that
when the solution of equations (39) with the indicated conditions for
$\bar{u}_{0}$ and $q_{0}$ reaches the stationary regime, it coincides with the
corresponding solution of stationary equations (44) everywhere, except for a
small neighbourhood of the cross-section $x=0$. To construct a numerical
solution of stationary or non-stationary equations, it is necessary to set a
small but positive value of the interlayer thickness $\eta$ at $x=0$.

###  4.2 The structure of stationary solutions

The conservative form of equations (39) allows one to construct solutions with
a hydraulic jump that transforms a supercritical flow into a subcritical one.
Across the jump the Rankine–Hugoniot relations are satisfied :

$[Q_{i}]=0,\quad[Q_{m}]=0,\quad[\bar{u}\eta\bar{b}]=0,\quad[J_{s}]=0,\quad[J]=0,\quad[W]=0,$
(47)

where

$\begin{array}[]{l}\displaystyle Q_{i}=u_{i}h_{i},\quad
Q_{m}=\bar{u}\eta,\quad J_{s}=\bar{u}q,\quad
J=\frac{u_{1}^{2}-u_{2}^{2}}{2}+bh_{1}+\bar{b}\eta,\\\\[8.53581pt]
\displaystyle
W=u_{1}^{2}h_{1}+(\bar{u}^{2}+q^{2})\eta+u_{2}^{2}\Big{(}h_{2}-\frac{H}{2}\Big{)}+\frac{bh_{1}^{2}}{2}+\Big{(}h_{1}+\frac{\eta}{2}\Big{)}\eta\bar{b}\end{array}$
(48)

$i=1,2$ and the square brackets mean the jump of the corresponding quantities.
These relations follow from conservation laws (39). By virtue of (47),
variables defined by formulas (48), as well as the buoyancy $\bar{b}$, are
continuous through the jump.

Let the supercritical flow ahead of the jump at $x=x_{s}-0$ be known. Then
behind the jump (values of the functions at $x=x_{s}\pm 0$ are denoted with
superscript ‘$\pm$’) we have

$u_{1}^{+}=\frac{Q_{1}}{h_{1}^{+}},\quad\bar{u}^{+}=\frac{Q_{m}}{\eta^{+}},\quad
u_{2}^{+}=\frac{Q_{2}}{h_{2}^{+}},\quad
q^{+}=\frac{J_{s}}{Q_{m}}\eta^{+}\quad(h_{2}^{+}=H-h_{1}^{+}-\eta^{+})$ (49)

and the following system of algebraic equations for finding $h_{1}^{+}$,
$\eta^{+}$:

$\begin{array}[]{l}\displaystyle\frac{1}{2}\bigg{(}\frac{Q_{1}^{2}}{(h_{1}^{+})^{2}}-\frac{Q_{2}^{2}}{(h_{2}^{+})^{2}}\bigg{)}+bh_{1}^{+}+\bar{b}\eta^{+}=J,\\\\[11.38109pt]
\displaystyle\frac{Q_{1}^{2}}{h_{1}^{+}}+\frac{Q_{2}^{2}}{h_{2}^{+}}\bigg{(}1-\frac{H}{2h_{2}^{+}}\bigg{)}+\frac{Q_{m}^{2}}{\eta^{+}}+\frac{J_{s}^{2}}{Q_{m}^{2}}(\eta^{+})^{3}+\frac{b}{2}(h_{1}^{+})^{2}+\bigg{(}h_{1}^{+}+\frac{\eta^{+}}{2}\bigg{)}\bar{b}\eta^{+}=W.\end{array}$
(50)

The quantities $Q_{i}$, $Q_{m}$, $J_{s}$, $J$ and $W$ are continuous across
the jump, so we know them a priori. Thus, the construction of a discontinuous
flow is reduced to finding a solution to algebraic equations (50)
corresponding to a subcritical flow.

Instead of the conservation law for $q$ one can also use the energy equation
(36). For this, we can eliminate the pressure $p^{*}$ from equation (36) by
subtracting the fifth equation in (35) multiplied by $\bar{Q}$. In this case,
instead of $J_{s}$, the variable

$M=\frac{u_{1}^{3}h_{1}}{2}+\frac{(\bar{u}^{2}+3q^{2})\bar{u}\eta}{2}+\frac{u_{2}^{2}(u_{2}h_{2}-\bar{Q})}{2}+(bh_{1}+\bar{b}\eta)u_{1}h_{1}+\bar{u}\bar{b}\eta(h_{1}+\eta)$

conserves across the jump. By combining $W$ and $M$ which are also conserved
across the jump, we can eliminate $q$ and obtain two algebraic equations for
$h_{1}^{+}$ and $\eta^{+}$. In this case, the second equation in (50) should
be replaced by

$\begin{array}[]{l}\displaystyle\frac{1}{2}\bigg{(}\frac{3Q_{m}}{\eta^{+}}-\frac{Q_{1}}{h_{1}^{+}}\bigg{)}\frac{Q_{1}^{2}}{h_{1}^{+}}+\frac{Q_{m}^{3}}{(\eta^{+})^{2}}+\bigg{(}\frac{3Q_{m}}{\eta^{+}}\Big{(}\frac{H}{2}-h_{1}^{+}-\eta^{+}\Big{)}-Q_{2}+\bar{Q}\bigg{)}\frac{Q_{2}^{2}}{2(h_{2}^{+})^{2}}\\\\[8.53581pt]
\displaystyle\quad\quad+\frac{3bQ_{m}(h_{1}^{+})^{2}}{4\eta^{+}}-(bh_{1}^{+}+\bar{b}\eta^{+})Q_{1}+\frac{Q_{m}\bar{b}}{2}\Big{(}h_{1}^{+}-\frac{\eta^{+}}{2}\Big{)}=\frac{3Q_{m}}{2\eta^{+}}W-M.\end{array}$
(51)

If $h_{1}^{+}$ and $\eta^{+}$ are defined, then the velocities behind the jump
are found from the first three relations (49). Then we find $q^{+}$ using the
continuity of $W$ across the jump, since $J_{s}$ is not conserved in this
case.

Let us consider the formation of a stationary mixing layer in a supercritical
flow over flat bottom ($Z=0$). All variables considered below are
dimensionless. We take $\kappa=6$, $H_{0}=10$ and $b=1$. At $x=0$ we take
$u_{10}=1.8$, $u_{20}=0.6$, $h_{10}=1$, $\eta_{0}=0.005$ (to carry out
computations, it is necessary to set a small but positive value of the
intermediate layer thickness). According to formulas (46) and (45) at $x=0$ we
obtain $\bar{u}_{0}=1.118$, $q_{0}=0.303$ (the second root of equation (46)
$\bar{u}_{0}=0.322$ lies outside the interval $(u_{10},u_{20})$ and therefore
this root is not physically admissible). In virtue of the third equation (44)
the buoyancy $\bar{b}_{0}=b/2=0.5$ remains constant for all $x>0$. It is easy
to verify that the data prescribed at $x=0$ correspond to supercritical flow,
since all roots $\lambda=\lambda_{i}$ of characteristic equation (41) are
positive. The corresponding continuous solution of equations (44) (interfaces
$z=h_{1}$ and $z=h_{1}+\eta$) are shown in Fig. 2 a by solid lines. As can be
seen from the graph, the monotonic expansion of the mixing layer takes place
up to $x\approx 14$. Then its quasi-periodic contraction and expansion occurs.
This behaviour is determined by a change in the sign of the variable $q$
(‘shear velocity’), which is responsible for the entrainment of liquid into
the intermediate mixing layer. The oscillatory nature of the solution is
illustrated in Fig. 2 b, which shows the considered solution in the plane
$(q,\bar{u}\eta)$ (solid curve). We note that for stationary solutions the
extrema of the flow rate $\bar{u}\eta$ in the mixing layer (points of maximum
and minimum fluid entrainment) correspond to $q=0$. When $x\to\infty$, the
solution tends towards equilibrium. This flow is supercritical everywhere.

Figure 2: Stationary mixing layer: a — solid curves 1, 2 are interfaces
$z=h_{1}$ and $z=h_{1}+\eta$ for continuous solution of equations (44), dashed
and dash-dotted curves 3 and 4 for discontinuous solution with a jump at
$x=10$ (the dashed lines correspond to the conservation of $\bar{u}q$, while
dash-dotted lines correspond to the conservation of energy); b — solid curve 5
and dashed curve 6 correspond to the continuous and discontinuous solutions on
the $(q,\bar{u}\eta)$-plane calculated on the interval $x\in(0,300)$.

As noted above, it is possible to construct a discontinuous solution. Suppose
that in the previous supercritical flow at some point $x=x_{s}$ there is a
shock. Let us choose $x_{s}=10$. Using the known values of the solution before
at $x_{x}-0$ ($h_{1}^{-}=0.869$, $\eta^{-}=0.649$, $u_{1}^{-}=1.689$,
$\bar{u}^{-}=1.028$, $q^{-}=0.135$) we solve equations (50) and find the
subcritical state behind the shock ($h_{1}^{+}=0.991$, $\eta^{+}=1.121$,
$u_{1}^{+}=1.490$, $\bar{u}^{+}=0.594$, $q^{+}=0.233$). The choice of this
solution comes from the analysis of the Rankine–Hugoniot relations (50) and
will be discussed later. Then we solve ODEs (44) with these data at $x=x_{s}$.
The obtained discontinuous solution is shown in Fig. 2 by dashed curves. This
flow is subcritical in the region $x>x_{s}$ and also has an oscillatory
character.

Nonlinear algebraic system (50) in the considered case has four admissible
solutions (such that $h_{1}^{+}>0$, $\eta^{+}>0$, $h_{1}^{+}+\eta^{+}<H$ and
$(h_{1}^{+},\eta^{+})\neq(h_{1}^{-},\eta^{-})$). One of them
($h_{1}^{+}=0.561$, $\eta^{+}=7.912$) transforms the supercritical flow into
the supercritical one. For another solution ($h_{1}^{+}=8.512$,
$\eta^{+}=0.037$) the characteristic equation (41) has complex roots. Both
solutions are likely to be unstable and are not realized. Finally, the last
admissible solutions of equations (50) ($h_{1}^{+}=0.991$, $\eta^{+}=1.121$,
and $h_{1}^{+}=3.051$, $\eta^{+}=4.946$) correspond to the subcritical flow
behind the shock. We chose the solution having the shock of smallest
amplitude. To justify this, we will use now non-stationary equations (39) and
will construct a numerical discontinuous solution arising in the case of flow
over an obstacle. By choosing the height of the obstacle, it is possible to
achieve a quasi-stationary regime in which the solution $h_{1}^{+}=0.991$ and
$\eta^{+}=1.121$ is realized.

Figure 3: Quasi-stationary discontinuous flow: a — interfaces $z=Z+h_{1}$ and
$z=H_{0}-h_{2}$ (curves 1, 2) and bottom topography $z=Z(x)$ (curve 3); b —
velocities $u_{1}$, $\bar{u}$ and $u_{2}$ in the layers (curves 4, 5 and 6)
and variable $q$ (curve 7). Solid curves — solution of equations (39) at
$t=400$, dashed curves — corresponding stationary solution.

It is interesting to note that the use of the energy equation across the jump,
i.e. of the system consisting of the first equation of (50) and (51) under the
same conditions ahead of the jump gives only two admissible solutions
$h_{1}^{+}=8.526$, $\eta^{+}=0.025$ and $h_{1}^{+}=0.950$, $\eta^{+}=0.991$.
The first of these solutions gives the complex value of the variable $q$ and
therefore is not physical. The second solution corresponds to
$u_{1}^{+}=1.544$, $\bar{u}^{+}=0.672$ and $q=0.302$. All roots of the
characteristic polynomial (41) are real. One of the roots is negative, the
other roots are positive. This corresponds to the subcritical flow regime
behind the jump. The corresponding downstream solution (for $x>x_{s}$) is
shown in Fig. 2 by a dash-dotted curve. It can be clearly seen that replacing
the cumbersome energy equation by its differential consequence does not lead
to a significant difference in the solution behaviour.

To conclude this Section, let us consider the formation of a quasi-stationary
flow with a jump over an obstacle. As before, we take $H_{0}=10$, $b=1$ and
set the same data $u_{10}=1.8$, $u_{20}=0.6$, $h_{10}=1$, $h_{20}=8.995$
$\eta_{0}=0.005$, $\bar{u}_{0}=1.118$, $q_{0}=0.303$ and $\bar{b}_{0}=0.5$ on
the left boundary of the computational domain $x=0$. The initial conditions
for unknowns are the same excepting $h_{1}$ for which one takes
$h_{1}(x,0)=h_{10}-Z(x)$, $x>0$. We set ‘soft’ conditions
$\mathbf{U}_{N}=\mathbf{U}_{N-1}$ on the right boundary of the computational
domain $x=25$. Here $\mathbf{U}_{j}$ is the value of unknown vector-function
$\mathbf{U}$ in the nodal point $x_{j}$. The bottom is chosen as

$Z(x)=z_{0}\exp(-0.25(x-x_{0})^{2})$

with $x_{0}=20$ and $z_{0}=0.135$. The bottom is required to create a left-
facing shock. The height of the obstacle $z_{0}$ is selected in such a way
that the shock propagates slowly or stops at some point $x=x_{s}$. This
corresponds to the formation of a stationary discontinuous solution. To solve
balance laws (39) numerically, we implemented the second-order central scheme
proposed by Nessyahu & Tadmor (1990). The calculation results on a uniform
grid with the number of nodes $N=1000$ on the interval $x\in(0,10)$ at $t=400$
are shown in Fig. 3 by solid curves. The position of the shock front is at
$x=10$. This shock solution corresponds to our stationary solution having the
smallest amplitude.

The corresponding discontinuous stationary solution of equations (44) is shown
in Fig. 3 by dashed lines. Visually, there is no difference between the
stationary solution and that obtained from non-stationary computations. We
note that the discontinuous solutions shown in Fig. 2 and 3) coincide, further
difference is related to the bottom topography. The flow shown in Fig. 3
continuously passes from subcritical regime to supercritical one at the point
$x_{c}\approx 19.7$. Thus, at the end of the computational domain the flow
again becomes supercritical.

##  5 Model validation

In this section we compare numerical simulation and experimental or field data
found in the literature. In subsections 5.1 and 5.2 we use the SGS system of
units (cm, g, s) while in 5.3 the International System of units is used (m,
kg, s).

###  5.1 Transcritical flow over an obstacle

We consider the formation of the mixing layer in a down-flow current of a
dense fluid in the framework of the proposed three-layer model. The
formulation of the problem is close to Pawlak & Armi (2000), where mixing and
entrainment during the evolution of stratified flows were experimentally
studied.

Figure 4: Flow over an obstacle: a — subcritical regime; b — supercritical
flow on the leeward side of the obstacle (curve 3). Solid curves 1 and 2 show
the interfaces $z=Z+h_{1}$ and $z=H_{0}-h_{2}$ at $t=t_{*}$ obtained by non-
stationary model (39), line 4 presents the sign of the characteristic function
$\chi(0)$. Dashed curves — corresponding stationary solution of equations
(44).

In this test we take $\kappa=6$ and choose $H_{0}=50$, $b=1.2$. Computations
are carried out in the domain $x\in[0,L]$, $L=140$. The bottom topography is
shown in Fig. 4 (curve 3). The obstacle is located on the interval
$x\in(10,117.5)$ and its maximum height $z_{\max}=14.5$ is reached at the
point $x=50$. The leeward side of the obstacle is rectilinear with a slope of
$k=-0.24$ (approximately $13.5^{\circ}$). Near the outlet section of the
channel ($x>130$), the bottom has a slope of $k=-0.5$. On the left boundary
$x=0$ we set $h_{1}=28$, $\eta=0.5$, $u_{1}=U$, $\bar{u}=U/2$, $u_{2}=0$,
$\bar{b}=b/2$ and $q=0$. First, we construct stationary solutions for
$U\approx 1.104$. At such a velocity, the flow becomes critical over the
obstacle since $\chi(0)$ vanishes at $x\approx 53.6$ (see the definition (41)
of the characteristic polynomial $\chi(\lambda)$). Fig. 4 a (dashed curves) is
obtained for $U=1.1039$ by using the standard ode45 procedure of the MATLAB
package that implements the fourth-order Runge–Kutta method. As can be seen
from the graph, intensive mixing occurs on the leeward side of the obstacle,
and fluid from the outer layers is entrained into the mixing layer. The fluid
velocities in the layers and the variable $q$, which is responsible for the
mixing process, are shown in Fig. 5 (solid curves). For $x<x_{c}$ ($x_{c}$ is
a critical point near the top of the obstacle, where $\chi(0)$ vanishes), the
flow velocity in the lower and intermediate layers increases. For $x>x_{c}$,
the heavy fluid continues to accelerate along the leeward slope, while the
fluid of intermediate density in the mixing layer, on the contrary, is
decelerated. The velocity of the light fluid in the upper layer changes
insignificantly, decreasing monotonically over the entire interval. In this
case, the solution is subcritical everywhere, with the exception of a small
area above the obstacle.

A slight increase of the lower layer velocity leads to a significant change of
the flow for $x>x_{c}$. Fig. 4 b (dashed curves) is obtained for $U=1.1052$.
With this velocity, the flow over the leeward side of the obstacle is
supercritical and the mixing layer thickness is noticeably smaller than in the
case of subcritical flow. Both the heavier fluid in the lower layer and the
lighter fluid in the intermediate mixing layer are accelerated over the
leeward slope of the obstacle (Fig. 5, dash-dotted curves). Before reaching
the flat bottom, the transition from supercritical to subcritical flow occurs
through a hydraulic jump at $x\approx 109.35$ (this point is chosen on the
basis of a non-stationary computation). In the subcritical flow region, the
mixing layer thickness increases, and the fluid velocity in the intermediate
layer slows down. Further, in the vicinity of the outlet section, the flow
again becomes supercritical. We note that by virtue of the third equation (44)
and the condition $\bar{b}=b/2$ at $x=0$, one obtains $\bar{b}=b/2$ in the
entire flow region.

Figure 5: Velocities in a transcritical stationary flow. Solid curves
correspond to the flow regime shown in Fig. 4 a, dash-dotted lines — Fig. 4 b.
Curves 1, 2 and 3 are the velocities in the lower, middle and upper layers;
curve 4 shows the variable $q$.

The constructed stationary solutions can also be obtained as a result of the
numerical solution of non-stationary equations (39). The computations are
carried out on a uniform grid with a number of nodes $N=1000$. As the initial
data at $t=0$, we take the above-mentioned values of the functions at $x=0$,
with the exception of $h_{1}$ and $u_{1}$, which are defined as
$h_{1}=h_{10}-Z(x)$, $u_{1}=(Q-\bar{u}_{0}\eta_{0})/h_{1}$. Here
$Q=(h_{10}+\eta_{0}/2)U$ and index ‘0’ corresponds to the values of functions
at $x=0$. The ‘soft’ conditions $\mathbf{U}_{N}=\mathbf{U}_{N-1}$ are imposed
on the right boundary $x=L$ of the computational domain. Numerical solutions
of equations (39) obtained for $U=1.10$ and $U=1.11$ are shown in Fig. 4 a and
Fig. 4 b, respectively. Both graphs are shown at $t=500$ s, which corresponds
to reaching the stationary regime. As can be seen from the figure, there is
good agreement between stationary and non-stationary computations, both for
the subcritical and supercritical regimes of flow over the leeward side of the
obstacle. Note that in non-stationary computations, the threshold value of the
velocity $U$ can vary slightly depending on the chosen numerical scheme and
spatial resolution.

###  5.2 Evolution of a mixing layer in a supercritical flow: comparison with
experiment

In this section we consider the evolution of the mixing layer flowing down the
inclined bottom and compare solutions of stationary equations (44) with the
experimental data by Pawlak & Armi (2000). In these experiments the channel
height was $H_{0}=30$ and a constant bottom slope was $10.8^{\circ}$. At
$x=x_{0}=10.5$ the depths of the layers and corresponding average velocities
are as follows : $h_{1}=6.25$, $\eta_{0}=0.75$, $h_{2}=8$ and $u_{1}=4$,
$\bar{u}=u_{1}/2$, $u_{2}=0$. We also choose $b=1.4$, $\bar{b}=b/2$. This is
consistent with the flow parameters in Pawlak & Armi (2000) (Fig. 9 therein).
Finally, $q$ for $x=x_{0}$ can take arbitrary value from the interval $[0,1]$.
This has no noticeable effect on the position of the interfaces $z=Z+h$ and
$z=Z+h+\eta$. To be specific, we take $q=0$. It is easy to verify that such a
choice of flow parameters at $x=x_{0}$ corresponds to the supercritical
regime, since $\chi(0)>0$. The dissipation parameter is taken as $\kappa=8$.

Figure 6: The interfaces $z=Z+h$ and $z=Z+h+\eta$ according to equations (44)
are shown by thick solid curves. Color picture presents experimental data of
the density gradient field images with velocity and density profiles (thin
transverse curves) Pawlak & Armi (2000) (Fig. 9 therein). The picture is
rotated so that the slope becomes horizontal.

In Fig. 6 bold solid lines show the boundaries $z=h_{1}$ and $z=h_{1}+\eta$ of
the mixing layer in the interval $x\in(10.5,72)$ determined by model (44). The
constant slope bottom has been rotated to become horizontal, and the abscissa
on Fig. 6 is the downslope distance. Coloured picture represents the density
gradient field image (blue is for a lower density gradient, red is for a
higher one) obtained by Pawlak & Armi (2000). As we can see, there is a fairly
good agreement between the experimental data and our numerical solution. We
have to mention that in the experiment (as well as in the nature), there is a
slight backward flow above the mixing layer. This fact is reflected in our
model. In general, the nature of the flow is similar to the example considered
above in the case of a supercritical flow over the leeward side of the
obstacle (see Fig. 4 b and 5).

Note that a similar comparison of the mixing layer boundaries with the
experimental data of Pawlak & Armi (2000) was carried out by Liapidevskii et
al. (2018) (see Fig. 2 in this work) based on a different (albeit similar)
depth-averaged model. The constitutive equations proposed by Liapidevskii et
al. (2018) better describe the mixing layer for $x\in(10,25)$, but give
excessively overestimated sizes of the mixing region for $x>30$. Our proposed
model gives better results for a sufficiently developed mixing layer.

###  5.3 Evolution of a mixing layer in a subcritical flow: comparison with
field observations

Field observations of stratified flows over topography in Knight Inlet
(British Columbia, Canada) are given by Farmer & Armi (1999); Cummins et al.
(2006). The stratified fluid flow over a sill illustrates the upstream
formation of a strong internal bore. In particular, a flow regime was recorded
over the leeward side of the obstacle similar to that shown in Fig. 4 a. An
image of acoustic backscatter obtained with a 200-kHz echosounder is presented
by Cummins et al. (2006) (see in Fig. 3 therein), where blue and red colours
correspond to 30 dB and 55 dB, respectively. This shows a plunging flow on the
lee side of the sill with large-amplitude instabilities on the interface of
the downslope flow. The Kelvin–Helmholtz instability develops as a result of
the large shear found along the interfaces. The formation of the flow
structure shown in the figure has been thoroughly documented by Farmer & Armi
(1999).

Figure 7: The interfaces $z=Z+h$ and $z=Z+h+\eta$ according to model (39) are
shown by bold solid curves. Coloured picture presents field data of acoustic
backscatter obtained by Cummins et al. (2006) (Fig. 3 therein).

We do not claim here to accurately reproduce the field data. Nevertheless,
choosing the flow parameters close to those considered by Cummins et al.
(2006), our model describes the characteristic features of the mixing layer
evolution. We perform computations using three-layer equations (39) in the
domain $x\in[-50,150]$ (in meters). The upper fluid level is constant and
equal to zero. On the interval $x\in[-50,0]$, the total depth $H_{0}-Z(x)$ is
$62$, and then, according to the parabolic law, increases to $80$ for
$x\in[0,72]$ and then takes this constant value. Buoyancy $b$ is $0.06$. On
the left boundary $x=-50$ of the calculated interval, the layers thicknesses,
corresponding velocities and buoyancy are set as follows $h_{1}=41$,
$\eta=3.5$ and $h_{2}=17.5$; $u_{1}=1$, $\bar{u}=0.59$ and $u_{2}=0.18$,
$\bar{b}=b/2$, $q=u_{2}$. The same values of the functions are chosen as the
initial conditions, excepting those for $h_{1}$ and $u_{1}$. The last are
determined taking into account the bottom topography and the constancy of the
total flow rate. As before, ‘soft’ boundary conditions are imposed at $x=150$.
The dissipation constant is taken as $\kappa=8$. Over time, the numerical
solution reaches a quasi-stationary regime. In fig. 7, bold solid lines show
the interface between the layer obtained according to model (39) at $t=4000$.
The coloured image on this figure shows the acoustic backscatter field data
Cummins et al. (2006) (Fig. 3 therein). As can be seen from this figure, the
interfaces predicted by our model are in a good agreement with the field data
of the acoustic backscatter Cummins et al. (2006), where the interfaces
correspond to the strongest acoustic response. We note that the obtained
numerical solution is subcritical $\chi(0)<0$ in all computational domain.
Obviously, this solution is similar to that considered above in the case of a
subcritical flow on the leeward side of the obstacle (see Fig. 4 a and Fig.
5).

##  6 Conclusion

A long wave approximation of the non-homogeneous Euler equations was used for
the construction of a depth-averaged model of three layer flows. The outer
layers are homogeneous and almost potential and can be described by Saint-
Venant-like equations (11) and (13). The fluid flow in the intermediate layer
is sheared and inhomogeneous, and is described by equations (24). The
interfaces separating these layers are considered as fronts through which the
turbulent mixing of fluids occurs. The general model (27) reveals the main
mechanisms of the mixing layer development. The model is reduced to the
classical equations of three-layer shallow water with a sheared intermediate
layer if the mass transfer between the layers is absent.

In the Boussinesq approximation, the proposed model is simplified and takes
the form (39), which allows a simple numerical implementation. It was proved
that for small relative velocities in the layers, the model is hyperbolic. The
concepts of supercritical, subcritical and transcritical flows are defined in
terms of the signs of the eigenvalues of the corresponding characteristic
polynomial. It makes it possible to formulate in classical terms the
conditions for the control of a subcritical flow by a downstream located
obstacle. The system of conservation laws representing our model uniquely
determines the jump relations. It is established that the use of an additional
conservation law for the shear velocity instead of the energy equation
simplifies the study of discontinuous solutions and does not lead to a
significant change in their structure. The model is able to control the mixing
process by changing the position of the stationary hydraulic jump upstream of
the obstacle. For supercritical and subcritical co–current flows new
oscillating structure of the mixing layer is found.

Considerable attention is paid to the validation of the proposed model (39) by
comparing its solutions with experimental data and field observations found in
the literature. We study transcritical three-layer flows over an obstacle.
This problem is close to that studied in Pawlak & Armi (2000), where mixing
and entrainment during the evolution of stratified flows were experimentally
explored. It is shown that the proposed model describes quite accurately the
development of the mixing layer on the leeward side of the obstacle in the
supercritical flow regime (Fig. 6). Field observations of the interaction of a
stratified flow with topography presented by Farmer & Armi (1999); Cummins et
al. (2006) were also used to validate the model. Comparison of the numerical
results with the observations over the Knight Inlet sill illustrating the
formation of a mixing layer, is shown in Fig. 7. These results demonstrate the
ability of our model to describe the characteristic features of gravity
currents.

## Acknowledgements

A.C. and V.L. were supported by the Russian Science Foundation (grant No.
21-71-20039). S.G. has been partially funded by the Excellence Initiative of
Aix-Marseille UniversityA*Midex, a French Investissements d’Avenir programme
AMX-19-IET-010.

## References

  * Armi & Farmer (2002) Army, L. & Farmer, D. 2002 Stratified flow over topography: bifurcation fronts and transition to the uncontrolled state. Proc. R. Soc. Lond. A 458, 513–538.
  * Baines (2005) Baines, P.G. 2005 Topographic Effects in Stratified Flows. Cambridge University Press.
  * Baines (2016) Baines, P.G. 2016 Internal hydraulic jumps in two-layer systems. J. Fluid Mech. 787, 1–15.
  * Barros et al. (2007) Barros, R., Gavrilyuk, S. & Teshukov, V. 2007 Dispersive nonlinear waves in two-layer flows with free surface. I. Model derivation and general properties. Stud. Appl. Math. 119, 191–211.
  * Chesnokov & Liapidevskii (2020) Chesnokov, A.A. & Liapidevskii, V.Yu. 2020 Mixing layer and turbulent jet flow in a Hele–Shaw cell. Int. J. Non-Linear. Mech. 125, 103534.
  * Chu & Baddour (1984) Chu, V.H. & Baddour, R. 1984 Turbulent gravity-stratified shear flows. J . Fluid Mech. 138, 353–378.
  * Cummins et al. (2006) Cummins, P.F., Armi, L. & Vagle, S. 2006 Upstream internal hydraulic jumps. J. Phys. Oceanogr. 36, 753–769.
  * Farmer & Armi (1999) Farmer, D. & Armi, L. 1999 Stratified flow over topography: The role of small-scale entrainment and mixing in flow establishment. Proc. R. Soc. Lond. A 455, 3221–3258.
  * Gavrilyuk et al. (2016) Gavrilyuk, S.L., Liapidevskii, V.Yu. & Chesnokov, A.A. 2016 Spilling breakers in shallow water: applications to Favre waves and to the shoaling and breaking of solitary waves. J. Fluid Mech. 808, 441–468.
  * Gavrilyuk et al. (2019) Gavrilyuk, S.L., Liapidevskii, V.Yu. & Chesnokov, A.A. 2019 Interaction of a subsurface bubble layer with long internal waves. Eur. J. Mech. B Fluids 73, 157–169.
  * Horsley & Woods (2018) Horsley, M.C. & Woods, A.W. 2018 A note on analytic solutions for entraining stratified gravity currents. J. Fluid Mech. 836, 260–276.
  * Jagannathan et al. (2017) Jagannathan, A., Winters, K.B. and Armi, L. 2017 Stability of stratified downslope flows with an overlying stagnant isolating layer. J. Fluid Mech. 810, 392–411.
  * Jagannathan et al. (2020) Jagannathan, A., Winters, K.B. & Armi, L. 2020 The effect of a strong density step on blocked stratified flow over topography. J. Fluid Mech. 889, A23.
  * Lamb (2004) Lamb, K.G. 2004 On boundary–-layer separation and internal wave generation at the Knight Inlet sill, Proc. R. Soc. Lond. A.  460, 2305–2337.
  * Liapidevskii (2004) Liapidevskii 2004 Mixing layer on the lee side of an obstacle. J. Appl. Mech. Tech. Phys. 45, 199–203.
  * Liapidevskii & Chesnokov (2014) Liapidevskii, V.Yu. & Chesnokov, A.A. 2014 Mixing layer under a free surface. J. Appl. Mech. Tech. Phys. 55, 299–310.
  * Liapidevskii et al. (2018) Liapidevskii, V.Yu., Dutykh, D. & Gisclon, M. 2018 On the modelling of shallow turbidity flows. Adv. Water Resour. 113, 310–327.
  * Lipatov et al. (2021) Lipatov, I.I., Liapidevskii, V.Yu. & Chesnokov, A.A. Forced oscillations of a pseudoshock in transonic gas flow in a diffuser, Fluid Dyn. 56 (2021) 860–869.
  * Nessyahu & Tadmor (1990) Nessyahu, H. & Tadmor, E. 1990 Non-oscillatory central differencing schemes for hyperbolic conservation laws. J. Comp. Phys. 87, 408–463.
  * Ogden & Helfrich (2016) Ogden, K.A. & Helfrich, K.R. 2016 Internal hydraulic jumps in two-layer flows with upstream shear. J. Fluid Mech. 789, 64–92.
  * Ogden & Helfrich (2020) Ogden, K.A. & Helfrich, K.R. 2020 Internal hydraulic jumps in two-layer flows with increasing upstream shear. Phys. Rev. Fluids 5, 074803.
  * Pawlak & Armi (1998) Pawlak, G. & Armi, L. 1998 Vortex dynamics in a spatially accelerating shear layer. J. Fluid Mech. 376, 1–35.
  * Pawlak & Armi (2000) Pawlak, G. & Armi, L. 2000 Mixing and entrainment in developing stratified currents. J. Fluid Mech. 424, 45–73.
  * Pope (2000) Pope, S. B. 2000 Turbulent Flows. Cambridge University Press.
  * Sher & Woods (2015) Sher, D. & Woods, A.W. 2015 Gravity currents: entrainment, stratification and self-similarity. J. Fluid Mech. 784, 130–162.
  * Sher & Woods (2017) Sher, D. & Woods, A.W. 2017 Mixing in continuous gravity currents. J. Fluid Mech. 818, R4.
  * Simpson (1997) Simpson, J. 1997 Gravity Currents. Cambridge University Press,.
  * Teshukov (2007) Teshukov, V.M. 2007 Gas-dynamic analogy in the theory of stratified liquid flows with a free boundary. Fluid Dynamics 42, 807–817.
  * Thorpe & Li (2014) Thorpe, S.A. & Li, L. Turbulent hydraulic jumps in a stratified shear flow. Part 2. J. Fluid Mech. 758, 94–120.
  * Thorpe et al. (2018) Thorpe, S.A., Malarkey, J., Voet, G., Alford, M.H. Girton, J.B. & Carter, G.S. 2018 Application of a model of internal hydraulic jumps. J. Fluid Mech. 834, 125–148.
  * Winters & Armi (2014) Winters, K.B. & Armi, L. 2014 Topographic control of stratified flows: upstream jets, blocking and isolating layers. J. Fluid Mech. 753, 80–103.
  * Yuan & Horner-Devine (2017) Yuan, Y. & Horner-Devine, A.R. Experimental investigation of large-scale vortices in a freely spreading gravity current. Phys. Fluids 29, 106603


# On Some Summation Formulas

Taekyun Kim  Department of Mathematics, Kwangwoon University, Seoul 139-701,
Republic of Korea<EMAIL_ADDRESS>,  DAE SAN KIM  Department of Mathematics,
Sogang University, Seoul 121-742, Republic of Korea<EMAIL_ADDRESS>,
Hyunseok Lee  Department of Mathematics, Kwangwoon University, Seoul 139-701,
Republic of Korea<EMAIL_ADDRESS>and  Jongkyum Kwon  Department of
Mathematics Education, Gyeongsang National University, Jinju 52828, Republic
of Korea<EMAIL_ADDRESS>
###### Abstract.

In this note, we derive a finite summation formula and an infinite summation
formula involving Harmonic numbers of order up to some order by means of
several definite integrals.

###### Key words and phrases:

summation formula; Harmonic numbers of order $r$

###### 2010 Mathematics Subject Classification:

11B83; 26A42

##  1\. Introduction

For $s\in\mathbb{C}$ with $\mathrm{Re}(s)>0$, the gamma function is defined by

(1) $\Gamma(s)=\int_{0}^{\infty}e^{-t}t^{s-1}dt,\quad(\mathrm{see}\ [5]).$

From (1), we note that $\Gamma(s+1)=s\Gamma(s)$ and $\Gamma(n+1)=n!$,
$(n\in\mathbb{N})$, (see [1,2,3]).

As is well known, the Beta function is defined for
$\mathrm{Re}\,\alpha>0,\,\mathrm{Re}\,\beta>0$ by

(2)
$B(\alpha,\beta)=\int_{0}^{1}t^{\alpha-1}(1-t)^{\beta-1}dt,\quad(\mathrm{see}\
[1,5]).$

From (2) we note that

(3)
$B(\alpha,\beta)=\frac{\Gamma(\alpha)\Gamma(\beta)}{\Gamma(\alpha+\beta)},\quad(\mathrm{see}\
[5]).$

The Harmonic numbers are defined by

(4) $H_{n}=1+\frac{1}{2}+\frac{1}{3}+\cdots+\frac{1}{n},\quad(n\geq 1),$

and more generally, for any $r\in\mathbb{N}$, the Harmonic numbers of order
$r$ are given by

$H_{n}^{(r)}=1+\frac{1}{2^{r}}+\frac{1}{3^{r}}+\cdots+\frac{1}{n^{r}},\quad(n\geq
1).$

In this note, we derive a finite summation formula in (5) and an infinite
summation formula (6), where $f_{n,r}(0)$ is determined by the recurrence
relation in (7) and given by Harmonic numbers of order $j$, for
$j=1,2,\dots,r$. Our results are illustrated for $r=3$ and $r=4$ in the
Example below. It is amusing that the infinite sum in (6) involving Harmonic
numbers of order $\leq r$ boils down to $(-1)^{r}(r+1)!$. They are derived
from several definite integrals in an elementary way. After considering the
two summation formulas in (5) and (6) for $r=1$ and $r=2$, we will show the
results in Theorem 1.

######  Theorem 1.

Let $n,r$ be positive integers. Then we have

(5)
$\displaystyle(-1)^{r}r!\sum_{k=0}^{n}\binom{n}{k}(-1)^{k}\frac{1}{(2k+1)^{r+1}}=\frac{2^{2n}}{(2n+1)\binom{2n}{n}}f_{n,r}(0),$
(6)
$\displaystyle\sum_{n=1}^{\infty}\frac{2^{2n-1}}{n(2n+1)\binom{2n}{n}}f_{n,r}(0)=(-1)^{r}(r+1)!,$

where $f_{n,s}(x)$ are determined by the recurrence relation

(7) $\displaystyle
f_{n,s+1}(x)=\beta_{n}(x)f_{n,s}(x)+\frac{d}{dx}f_{n,s}(x),\,\,(s\geq
1),\,\,f_{n,1}(x)=\beta_{n}(x)=-\sum_{k=0}^{n}(x+2k+1)^{-1},$

so that $f_{n,s}(x)$ is a polynomial in $x$ involving
$\beta_{n}(x),\beta_{n}^{(1)}(x),\cdots,\beta_{n}^{(s-1)}(x)$, with

(8)
$\displaystyle\beta_{n}^{(j)}(0)=(-1)^{j-1}j!\big{(}H_{2n+1}^{(j+1)}-\frac{1}{2^{j+1}}H_{n}^{(j+1)}\big{)},\quad(j\geq
0).$

Example:
(a)

$\displaystyle-3!\sum_{k=0}^{n}\binom{n}{k}(-1)^{k}\frac{1}{(2k+1)^{4}}=\frac{2^{2n}}{(2n+1)\binom{2n}{n}}$
$\displaystyle\times\left\\{-(H_{2n+1}-\frac{1}{2}H_{n})^{3}-3(H_{2n+1}-\frac{1}{2}H_{n})(H_{2n+1}^{(2)}-\frac{1}{4}H_{n}^{(2)})-2(H_{2n+1}^{(3)}-\frac{1}{8}H_{n}^{(3)})\right\\},$
$\displaystyle\sum_{n=1}^{\infty}\frac{2^{2n-1}}{n(2n+1)\binom{2n}{n}}\left\\{-(H_{2n+1}-\frac{1}{2}H_{n})^{3}-3(H_{2n+1}-\frac{1}{2}H_{n})(H_{2n+1}^{(2)}-\frac{1}{4}H_{n}^{(2)})-2(H_{2n+1}^{(3)}-\frac{1}{8}H_{n}^{(3)})\right\\}$
$\displaystyle=-4!.$

(b)

$\displaystyle
4!\sum_{k=0}^{n}\binom{n}{k}(-1)^{k}\frac{1}{(2k+1)^{5}}=\frac{2^{2n}}{(2n+1)\binom{2n}{n}}\left\\{(H_{2n+1}-\frac{1}{2}H_{n})^{4}+6(H_{2n+1}-\frac{1}{2}H_{n})^{2}(H_{2n+1}^{(2)}-\frac{1}{4}H_{n}^{(2)})\right.$
$\displaystyle\quad\quad\quad\quad\quad\quad\quad\quad\left.3(H_{2n+1}^{(2)}-\frac{1}{4}H_{n}^{(2)})^{2}+8(H_{2n+1}-\frac{1}{2}H_{n})(H_{2n+1}^{(3)}-\frac{1}{8}H_{n}^{(3)})+6(H_{2n+1}^{(4)}-\frac{1}{16}H_{n}^{(4)})\right\\},$
$\displaystyle
5!=\sum_{n=1}^{\infty}\frac{2^{2n-1}}{n(2n+1)\binom{2n}{n}}\left\\{(H_{2n+1}-\frac{1}{2}H_{n})^{4}+6(H_{2n+1}-\frac{1}{2}H_{n})^{2}(H_{2n+1}^{(2)}-\frac{1}{4}H_{n}^{(2)})\right.$
$\displaystyle\quad\quad\quad\quad\quad\quad\quad\quad\left.3(H_{2n+1}^{(2)}-\frac{1}{4}H_{n}^{(2)})^{2}+8(H_{2n+1}-\frac{1}{2}H_{n})(H_{2n+1}^{(3)}-\frac{1}{8}H_{n}^{(3)})+6(H_{2n+1}^{(4)}-\frac{1}{16}H_{n}^{(4)})\right\\}.$

##  2\. Derivation of Summation Formulas

First, we observe that

(9) $\displaystyle\int_{0}^{1}(1-x^{2})^{n}dx$
$\displaystyle=\sum_{k=0}^{n}\binom{n}{k}(-1)^{k}\int_{0}^{1}x^{2k}dx$
$\displaystyle=\sum_{k=0}^{n}\binom{n}{k}(-1)^{k}\frac{1}{2k+1}.$

On the other hand,

(10) $\displaystyle\int_{0}^{1}(1-x^{2})^{n}dx$
$\displaystyle=\frac{1}{2}\int_{0}^{1}(1-y)^{n}y^{-\frac{1}{2}}dy=\frac{1}{2}\int_{0}^{1}(1-y)^{n+1-1}\cdot
y^{\frac{1}{2}-1}dy$
$\displaystyle=\frac{1}{2}B\bigg{(}n+1,\frac{1}{2}\bigg{)}=\frac{1}{2}\frac{\Gamma(n+1)\Gamma\big{(}\frac{1}{2}\big{)}}{\big{(}n+\frac{1}{2}\big{)}\big{(}n-\frac{1}{2}\big{)}\big{(}n-\frac{3}{2}\big{)}\cdots\frac{1}{2}\Gamma\big{(}\frac{1}{2}\big{)}}$
$\displaystyle=\frac{2^{n}n!}{(2n+1)(2n-1)(2n-3)\cdots
1}=\frac{2^{2n}n!n!}{(2n+1)(2n)!}=\frac{2^{2n}}{(2n+1)\binom{2n}{n}}.$

Thus, from (9) and (10) we note that

(11)
$\binom{2n}{n}\sum_{k=0}^{n}\binom{n}{k}\frac{(-1)^{k}}{2k+1}=\frac{2^{2n}}{2n+1}.$

From (11), we have

(12)
$\sum_{n=1}^{\infty}\frac{2^{2n}}{n(2n+1)\binom{2n}{n}}=\sum_{n=1}^{\infty}\frac{1}{n}\int_{0}^{1}(1-x^{2})^{n}dx=-2\int_{0}^{1}\log
x\ dx=2.$

Let

(13) $\displaystyle F(x)$
$\displaystyle=\int_{0}^{1}(1-t^{2})^{n}t^{x}dt=\sum_{k=0}^{n}\binom{n}{k}(-1)^{k}\int_{0}^{1}t^{2k+x}dt=\sum_{k=0}^{n}\binom{n}{k}(-1)^{k}\frac{1}{2k+x+1}.$

Note that

(14) $\displaystyle
F^{\prime}(x)=\frac{d}{dx}\int_{0}^{1}(1-t^{2})^{n}t^{x}dt=\int_{0}^{1}(1-t^{2})^{n}t^{x}\log
tdt=-\sum_{k=0}^{n}\binom{n}{k}(-1)^{k}\frac{1}{(2k+1+x)^{2}}.$

On the other hand

(15) $\displaystyle\int_{0}^{1}(1-t^{2})^{n}t^{x}dt$
$\displaystyle=\frac{1}{2}\int_{0}^{1}(1-y)^{n}y^{\frac{x+1}{2}-1}dy=\frac{1}{2}B\bigg{(}n+1,\frac{x+1}{2}\bigg{)}$
$\displaystyle=\frac{1}{2}\frac{\Gamma(n+1)\Gamma\big{(}\frac{x+1}{2}\big{)}}{\Gamma\big{(}n+1+\frac{x+1}{2}\big{)}}=\frac{1}{2}\frac{n!\Gamma\big{(}\frac{x+1}{2}\big{)}}{\big{(}n+\frac{x+1}{2}\big{)}\big{(}n+\frac{x-1}{2}\big{)}\big{(}n+\frac{x-3}{2}\big{)}\cdots\big{(}\frac{x+1}{2}\big{)}\Gamma\big{(}\frac{x+1}{2}\big{)}}$
$\displaystyle=\frac{n!2^{n}}{(2n+x+1)(2n+x-1)(2n+x-3)\cdots(x+1)}=\frac{n!2^{n}}{\prod_{k=0}^{n}(x+2k+1)}.$

Thus, we have

(16) $\displaystyle F^{\prime}(x)$
$\displaystyle=\frac{d}{dx}F(x)=\frac{d}{dx}\int_{0}^{1}(1-t^{2})^{n}t^{x}dx=\frac{d}{dx}\bigg{(}\frac{n!2^{n}}{\prod_{k=0}^{n}(x+2k+1)}\bigg{)}$
$\displaystyle=-\bigg{(}\frac{n!2^{n}}{\prod_{k=0}^{n}(x+2k+1)}\bigg{)}\bigg{(}\sum_{k=0}^{n}\frac{1}{x+2k+1}\bigg{)}.$

Thus, by (16), we get

(17)
$F^{\prime}(0)=-\bigg{(}\frac{n!n!2^{2n}}{(2n+1)(2n)!}\bigg{)}\bigg{(}\sum_{k=0}^{n}\frac{1}{2k+1}\bigg{)}=-\frac{2^{2n}}{(2n+1)\binom{2n}{n}}\bigg{(}H_{2n+1}-\frac{1}{2}H_{n}\bigg{)}.$

From (14) and (17), we note that

(18)
$\displaystyle\binom{2n}{n}\sum_{k=0}^{n}\binom{n}{k}(-1)^{k}\frac{1}{(2k+1)^{2}}=\frac{2^{2n}}{2n+1}\bigg{(}H_{2n+1}-\frac{1}{2}H_{n}\bigg{)}.$

By (18), we have

(19)
$\displaystyle\sum_{n=1}^{\infty}\frac{2^{2n}}{n(2n+1)\binom{2n}{n}}\big{(}H_{2n+1}-\frac{1}{2}H_{n}\big{)}$
$\displaystyle=\sum_{n=1}^{\infty}\frac{1}{n}\sum_{k=0}^{n}\binom{n}{k}\frac{(-1)^{k}}{(2k+1)^{2}}$
$\displaystyle=\sum_{n=1}^{\infty}\frac{1}{n}\int_{0}^{1}\int_{0}^{1}(1-x^{2}y^{2})^{n}dxdy$
$\displaystyle=-2\int_{0}^{1}\int_{0}^{1}(\log x+\log y)dxdy=4.$

Thus, we have

(20)
$\displaystyle\sum_{n=1}^{\infty}\frac{2^{2n}}{n(2n+1)\binom{2n}{n}}\bigg{(}H_{2n+1}-\frac{1}{2}H_{n}\bigg{)}=4.$

From (13), we note that

(21) $\displaystyle F^{\prime\prime}(x)=\frac{d^{2}}{dx^{2}}F(x)$
$\displaystyle=\frac{d^{2}}{dx^{2}}\int_{0}^{1}(1-t^{2})^{n}t^{x}dt$
$\displaystyle=\int_{0}^{1}(1-t^{2})^{n}\big{(}\log t\big{)}^{2}t^{x}dt$
$\displaystyle=2!\sum_{k=0}^{n}\binom{n}{k}(-1)^{k}\frac{1}{(2k+1+x)^{3}}.$

Hence, by (21), we get

(22) $F^{\prime\prime}(x)=\int_{0}^{1}(1-t^{2})^{n}\big{(}\log
t\big{)}^{2}t^{x}dt=2!\sum_{k=0}^{n}\binom{n}{k}(-1)^{k}\frac{1}{(2k+1+x)^{3}}.$

From (15) and (16), we have

(23) $\displaystyle F^{\prime\prime}(x)$
$\displaystyle=\frac{d^{2}}{dx^{2}}\int_{0}^{1}(1-t^{2})^{n}t^{x}dt$
$\displaystyle=\frac{n!2^{n}}{(2n+1+x)(2n+x-1)\cdots(x+1)}\bigg{\\{}\bigg{(}\sum_{k=0}^{n}\frac{1}{2k+x+1}\bigg{)}^{2}+\sum_{k=0}^{n}\frac{1}{(2k+x+1)^{2}}\bigg{\\}}.$

Thus, we note that

(24) $\displaystyle F^{\prime\prime}(0)$
$\displaystyle=\frac{n!2^{n}}{(2n+1)(2n-1)(2n-3)\cdots
1}\bigg{\\{}\bigg{(}\sum_{k=0}^{n}\frac{1}{2k+1}\bigg{)}^{2}+\sum_{k=0}^{n}\frac{1}{(2k+1)^{2}}\bigg{\\}}$
$\displaystyle=\frac{n!n!2^{2n}}{(2n+1)(2n)!}\bigg{\\{}\big{(}H_{2n+1}-\frac{1}{2}H_{n}\big{)}^{2}+\bigg{(}H_{2n+1}^{(2)}-\frac{1}{4}H_{n}^{(2)}\bigg{)}\bigg{\\}}.$

Therefore, by (21) and (24), we get

(25) $\displaystyle
2!\binom{2n}{n}\sum_{k=0}^{n}\binom{n}{k}\frac{(-1)^{k}}{(2k+1)^{3}}=\frac{2^{2n}}{(2n+1)}\bigg{\\{}\bigg{(}H_{2n+1}-\frac{1}{2}H_{n}\bigg{)}^{2}+\bigg{(}H_{2n+1}^{(2)}-\frac{1}{4}H_{n}^{(2)}\bigg{)}\bigg{\\}}.$

Note that

(26)
$\displaystyle\sum_{n=1}^{\infty}\frac{1}{n}\bigg{\\{}\sum_{k=0}^{n}\binom{n}{k}\frac{(-1)^{k}}{(2k+1)^{3}}\bigg{\\}}$
$\displaystyle=\sum_{n=1}^{\infty}\frac{1}{n}\int_{0}^{1}\int_{0}^{1}\int_{0}^{1}(1-x^{2}y^{2}z^{2})^{n}dxdydz$
$\displaystyle=-2\int_{0}^{1}\int_{0}^{1}\int_{0}^{1}(\log x+\log y+\log
z)dxdydz=6.$

From (25) and (26), we have

(27)
$\displaystyle\sum_{n=1}^{\infty}\frac{2^{2n-1}}{n(2n+1)\binom{2n}{n}}\bigg{\\{}\bigg{(}H_{2n+1}-\frac{1}{2}H_{n}\bigg{)}^{2}+\bigg{(}H_{2n+1}^{(2)}-\frac{1}{4}H_{n}^{(2)}\bigg{)}^{2}\bigg{\\}}=6.$

Now, we begin to prove Theorem 1. Let $F(x)=\int_{0}^{1}(1-t^{2})^{n}t^{x}dt$
be as in (13).
Let $r$ be a positive integer. Then repeated integrating by parts gives us

$\displaystyle F^{(r)}(x)$
$\displaystyle=\bigg{(}\frac{d}{dx}\bigg{)}^{r}\int_{0}^{1}(1-t^{2})^{n}t^{x}dt=\int_{0}^{1}(1-t^{2})^{n}(\log
t)^{r}t^{x}dt$
$\displaystyle=\sum_{k=0}^{n}\binom{n}{k}(-1)^{k}\int_{0}^{1}(\log
t)^{r}t^{x+2k}dt$ (28)
$\displaystyle=\sum_{k=0}^{n}\binom{n}{k}(-1)^{k}\frac{-r}{x+2k+1}\int_{0}^{1}(\log
t)^{r-1}t^{x+2k}dt$
$\displaystyle=\sum_{k=0}^{n}\binom{n}{k}(-1)^{k}\frac{-r}{x+2k+1}\cdots\frac{-1}{x+2k+1}\int_{0}^{1}t^{x+2k}dt$
$\displaystyle=(-1)^{r}r!\sum_{k=0}^{n}\binom{n}{k}(-1)^{k}\frac{1}{(x+2k+1)^{r+1}}.$

As we saw in (15) and (16), $F(x)$ is alternatively expressed by

(29) $\displaystyle
F(x)=\int_{0}^{1}(1-t^{2})^{n}t^{x}dt=\frac{n!2^{n}}{\prod_{k=0}^{n}(x+2k+1)},\,\,\,F(0)=\frac{2^{2n}}{(2n+1)\binom{2n}{n}},$

and its derivative is given by

(30) $\displaystyle
F^{(1)}(x)=F(x)\beta_{n}(x),\quad\beta_{n}=\beta_{n}(x)=-\sum_{k=0}^{n}(x+2k+1)^{-1}.$

Repeated differentiations give us

$\displaystyle
F^{(2)}(x)=F(x)(\beta_{n}^{2}+\beta_{n}^{(1)}),\,\,F^{(3)}(x)=F(x)(\beta_{n}^{3}+3\beta_{n}\beta_{n}^{(1)}+\beta_{n}^{(2)}),$
(31) $\displaystyle
F^{(4)}(x)=F(x)(\beta_{n}^{4}+6\beta_{n}^{2}\beta_{n}^{(1)}+3(\beta_{n}^{(1)})^{2}+4\beta_{n}\beta_{n}^{(2)}+\beta_{n}^{(3)}),$
$\displaystyle
F^{(5)}(x)=F(x)\left(\beta_{n}^{5}+10\beta_{n}^{3}\beta_{n}^{(1)}+10\beta_{n}^{2}\beta^{(2)}+15\beta(\beta_{n}^{(1)})^{2}\right.$
$\displaystyle\left.\quad\quad\quad\quad\quad\quad\quad+5\beta_{n}\beta_{n}^{(3)}+10\beta_{n}^{(1)}\beta_{n}^{(2)}+\beta_{n}^{(4)}\right),\dots.$

In general, we let $F^{(s)}(x)=F(x)f_{n,s}(x)$, for $s\geq 1$. Then further
differentiation gives us

$\displaystyle
F^{(s+1)}(x)=F(x)(\beta_{n}(x)f_{n,s}(x)+\frac{d}{dx}f_{n,s}(x)).$

Thus we get the recurrence relation for
$\left\\{f_{n,s}(x)\right\\}_{s=1}^{\infty}$:

(32) $\displaystyle
f_{n,s+1}(x)=\beta_{n}(x)f_{n,s}(x)+\frac{d}{dx}f_{n,s}(x),\,\,(s\geq
1),\,\,f_{n,1}(x)=\beta_{n}(x).$

Now, from (28) and (29), we obtain

(33)
$\displaystyle(-1)^{r}r!\sum_{k=0}^{n}\binom{n}{k}(-1)^{k}\frac{1}{(x+2k+1)^{r+1}}=\frac{n!2^{n}}{\prod_{k=0}^{n}(x+2k+1)}f_{n,r}(x).$

Letting $x=0$ in (33), we get

(34)
$\displaystyle\sum_{k=0}^{n}\binom{n}{k}(-1)^{k}\frac{1}{(2k+1)^{r+1}}=\frac{(-1)^{r}}{r!}\frac{2^{2n}}{(2n+1)\binom{2n}{n}}f_{n,r}(0).$

Multiplying both sides of (34) by $\frac{1}{n}$ and summing over $n$, we get

$\displaystyle\frac{(-1)^{r}}{r!}\sum_{n=1}^{\infty}\frac{2^{2n}}{n(2n+1)\binom{2n}{n}}f_{n,r}(0)$
$\displaystyle=\sum_{n=1}^{\infty}\frac{1}{n}\sum_{k=0}^{n}\binom{n}{k}(-1)^{k}\frac{1}{(2k+1)^{r+1}}$
(35)
$\displaystyle=\sum_{n=1}^{\infty}\frac{1}{n}\int_{0}^{1}\cdots\int_{0}^{1}(1-x_{1}^{2}x_{2}^{2}\cdots
x_{r+1}^{2})^{n}dx_{1}dx_{2}\cdots dx_{r+1}$
$\displaystyle=-2\int_{0}^{1}\cdots\int_{0}^{1}(\log x_{1}+\cdots+\log
x_{r+1})dx_{1}\cdots dx_{r+1}$ $\displaystyle=2(r+1).$

Finally, we note that $f_{n,s}(x)$ is a polynomial in $x$ involving
$\beta_{n}(x),\beta_{n}^{(1)}(x),\cdots,\beta_{n}^{(s-1)}(x)$, as we can see
from the recurrence relation in (32). From (30), we see that

(36)
$\displaystyle\beta_{n}^{(j)}(0)=(-1)^{j-1}j!\sum_{k=0}^{n}(2k+1)^{-(j+1)}=(-1)^{j-1}j!\big{(}H_{2n+1}^{(j+1)}-\frac{1}{2^{j+1}}H_{n}^{(j+1)}\big{)}.$

Combining (32) and (34)-(36) altogether, we obtain our main result in Theorem
1.

## References

  * [1] Andrews, G. E.; Askey, R.; Roy, R. _Special functions, Encyclopedia of Mathematics and its Applications_ 71, Cambridge University Press, Cambridge, 1999.
  * [2] Araci, S.; Acikgoz, M.; Özel, C.; Srivastava, H. M.; Diagana, T. _Recent trends in special numbers and special functions and polynomials_ [Editorial], Int. J. Math. Math. Sci. 2015, Art. ID 573893, 1 p.
  * [3] Kim, T.; Kim, D. S. _Some Identities on Truncated Polynomials Associated with Degenerate Bell Polynomials,_ Russ. J. Math. Phys. 28 (2021), no. 3, 342–355.
  * [4] Magnus, W.; Oberhettinger, F. _Formulas and Theorems for the Special Functions of Mathematical Physics,_ Translated by John Wermer, Chelsea Publishing Company, New York, N.Y., 1949.
  * [5] Whittaker, E. T.; Watson, G. N. _A course of modern analysis, An introduction to the general theory of infinite processes and of analytic functions; with an account of the principal transcendental functions,_ Reprint of the fourth (1927) edition. Cambridge Mathematical Library. Cambridge University Press, Cambridge, 1996.