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

HAL
open science

Analysis of the Riemann Problem for a shallow water model with two velocities

Nina Aguillon, Emmanuel Audusse, Edwige Godlewski, Martin Parisot

► To cite this version:

Nina Aguillon, Emmanuel Audusse, Edwige Godlewski, Martin Parisot. Analysis of the Riemann Problem for a shallow water model with two velocities. SIAM Journal on Mathematical Analysis, 2018, 10.1137/17M1152887 . hal-01618722v2

HAL Id: hal-01618722

https://inria.hal.science/hal-01618722v2

Submitted on 19 Oct 2017

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

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
Analysis of the Riemann Problem for a shallow
water model with two velocities

Nina Aguillon*¹,², Emmanuel Audusse†³, Edwige Godlewski‡¹,²,⁴,
and Martin Parisot§⁴,¹,²

¹Sorbonne Universités, UPMC Univ Paris 06, UMR 7598,
Laboratoire Jacques-Louis Lions, F-75005, Paris, France

²CNRS, UMR 7598, Laboratoire Jacques-Louis Lions, F-75005,
Paris, France

³Université Paris 13, Laboratoire d'Analyse, Géométrie et
Applications, 99 av. J.-B. Clément, F-93430 Villetaneuse, France

⁴INRIA Paris, ANGE Project-Team, 75589 Paris Cedex 12, France

October 18, 2017

Abstract

Some shallow water type models describing the vertical profile of the
horizontal velocity with several degrees of freedom have been recently pro-
posed. The question addressed in the current work is the hyperbolicity
of a shallow water model with two velocities. The model is written in a
nonconservative form and the analysis of its eigenstructure shows the pos-
sibility that two eigenvalues coincide. A definition of the nonconservative
product is given which enables us to analyse the resonance and coalescence
of waves. Eventually, we prove the well-posedness of the two dimensional
Riemann problem with initial condition constant by half-plane.

Introduction

In this work, we are interested in the solution of the planar Riemann problem
for an enriched two dimensional shallow water type model where the horizontal

<footnote>* nina.aguillon@upmc.fr</footnote>

<footnote>† audusse@math.univ-paris13.fr</footnote>

<footnote>‡edwige.godlewski@upmc.fr</footnote>

<footnote>§martin.parisot@inria.fr</footnote>
velocity is described by two variables. More precisely, the system reads

$$ (1) \quad \left\{ \begin{aligned} \partial_t h &+ \nabla \cdot (h \bar{\mathbf{U}}) & = 0, \\ \partial_t (h \bar{\mathbf{U}}) &+ \nabla \cdot \left(h \left(\bar{\mathbf{U}} \otimes \bar{\mathbf{U}} + \hat{\mathbf{U}} \otimes \hat{\mathbf{U}}\right) + \frac{g}{2} h^2 I_d\right) & = 0, \\ \partial_t \hat{\mathbf{U}} &+ (\bar{\mathbf{U}} \cdot \nabla) \hat{\mathbf{U}} + (\hat{\mathbf{U}} \cdot \nabla) \bar{\mathbf{U}} & = 0, \end{aligned} \right. $$

with the initial data $(h, \bar{\mathbf{U}}, \hat{\mathbf{U}})(0, x, y) = (h^0, \bar{\mathbf{U}}^0, \hat{\mathbf{U}}^0)(x, y)$. By planar Riemann problem, we mean the case of an initial data which is constant by half-plane, i.e.

$$ (2) \quad \forall (x, y) \in \mathbb{R}^2, \quad (h, \bar{\mathbf{U}}, \hat{\mathbf{U}})(0, x, y) = \begin{cases} (h_L, \bar{\mathbf{U}}_L, \hat{\mathbf{U}}_L) & \text{if } x < 0 \\ (h_R, \bar{\mathbf{U}}_R, \hat{\mathbf{U}}_R) & \text{if } x > 0. \end{cases} $$

Here $h(t, x, y) \in \mathbb{R}_+$ denotes the water depth of the flow and $\bar{\mathbf{U}}(t, x, y) \in \mathbb{R}^2$ and $\hat{\mathbf{U}}(t, x, y) \in \mathbb{R}^2$ refer to its horizontal velocity and denote respectively the vertical-average and the oriented standard deviation along the vertical axis (the introduction of these quantities will be explained below, in (5)). Setting the standard deviation to zero, one retrieves the classical shallow water model for which the solution of the planar Riemann problem is well known [25]. The solution of the planar Riemann problem (1)-(2) is a major improvement in the understanding of Model (1). In particular it is closely related to the implementation of the finite volume method on two dimensional meshes since it can be used to construct Godunov type schemes [13, 16, 25]. But the solution of the Riemann problem (1)-(2) is far from being obvious since System (1) is not in conservative form. The definition of the nonconservative products have to be precised for solutions that contain discontinuities. One can for example use the theory developed in [8]. Here it is possible to use the simpler approach developed in [1, 15] that gives a path-independent definition to the nonconservative products. Furthermore, it is well known that such a nonconservative system may lead to resonant solutions of the Riemann problem [7, 12, 14]. Here we will show that two kinds of such phenomena may arise: the first one is a resonance between two linearly degenerate fields and the second is the coalescence, in a sense that will be defined later, of a contact discontinuity with a nonlinear shock wave. These two situations are analyzed in detail and the solution of the Riemann problem is then entirely computed for all initial data.

The main reason we are interested in Model (1) is that it is equivalent, at least for continuous solution, to the bilayer version of the layerwise model introduced in [4]. Such layerwise models were introduced as a way to approximate the three dimensional hydrostatic Euler equations by using two dimensional models, but avoiding the shallow flow hypothesis. Different models can be derived, depending on the closure that is chosen for the definition of the layer, see [2, 21]. Other multilayer models was proposed in the literature for stratified flow [17, 20] but does not enter in the scope of this work. In the layerwise discretization [4] that we consider here, the layer thickness is assumed to be proportional to the total water depth. This model can be interpreted as an ALE
(for Arbitrary-Lagrangian-Eulerian) vertical discretization and allows large de-
formations of the free surface as long as there is no breaking, i.e. the free surface
remains a monovaluated function in space. Several generalizations of the model
have been proposed for variable density flow [3], for hydrostatic Navier-Stokes
equations [5] and for non hydrostatic flows [10]. In this context, the solution
of the Riemann problem (1)-(2) appears as an important step in a better un-
derstanding of this now widely used family of free surface fluid models. For
example, in 1D, System (1) becomes conservative and then can be seen as a
natural way to define the nonconservative products involved in the layerwise
model, introduced below (3). But Model (1) is also, as such, an improvement
in the modeling of free surface shallow flows since it involves a new degree of
freedom that takes into account the vertical shear in the flow. Similar models
have been recently proposed to analyze particular regimes such as roll waves or
hydraulic jumps in 1D [22, 24] and in 2D [11], see also [6] for an extension to non
hydrostatic flows. Up to our knowledge, it is the first time that the Riemann
problem is analyzed in details for a shallow water model with two velocities.

The paper is organized as follows. In section 1, we first detail similarities
between Model (1) and other models mentioned in the Introduction, in particu-
lar [4, 11, 22]. Then we exhibit its main properties. In particular we show that
the one dimensional Riemann problem for the normal velocity plays a role in the
solution of the two-dimensional planar Riemann problem. Then, in Section 2,
this one-dimensional Riemann problem is entirely solved. This part does not
present particular difficulties since the one-dimensional model is conservative.
The analogy with Euler equations is also discussed. In Section 3 we come back
to the two-dimensional case. We first analyze the hyperbolicity of the model
and characterize the nature of the wave associated to each eigenvalue. Then
we propose a constructive definition of the nonconservative products for solu-
tions containing (possibly nonlinear) discontinuities by considering a regularized
problem, yielding to jump conditions that are path independent [15]. Equiped
with these notions, we are able to detail the solution of the planar Riemann
problem, first without coalescence nor resonance, then by including coalescence
phenomenon, and finally by taking into account the possibility of two resonant
waves. In Section 4 some analytical solutions, including coalescence and reso-
nance phenomena, are presented and discussed.

# 1 The shallow water model with two velocities

1.1 Relations with other shallow water type models

Let us first precise the similarities between the shallow water model with two
velocity (1) and two other models mentioned in the Introduction, i.e. the model
for shear shallow flow introduced in [22] and extended in 2D in [11] and the
bilayer version of the layerwise model for hydrostatic flows introduced in [4].
Let us begin by considering this latter model [4] which is a layerwise discretized
Euler system. It is derived by vertical averaging of the Euler equations between
artificial interfaces. More precisely, let $t \in \mathbb{R}_+$ and $(x, y) \in \mathbb{R}^2$ be the space variables, the 2D homogeneous bilayer model reads

$$ (3) \quad \left\{ \begin{array}{l} \partial_t h_1 + \nabla \cdot (h_1 \mathbf{U}_1) = G_{3/2}, \\ \partial_t h_2 + \nabla \cdot (h_2 \mathbf{U}_2) = -G_{3/2}, \\ \partial_t (h_1 \mathbf{U}_1) + \nabla \cdot \left(h_1 \mathbf{U}_1 \otimes \mathbf{U}_1 + \frac{g}{2} h_1 h I_d\right) = \mathbf{U}_{3/2} G_{3/2}, \\ \partial_t (h_2 \mathbf{U}_2) + \nabla \cdot \left(h_2 \mathbf{U}_2 \otimes \mathbf{U}_2 + \frac{g}{2} h_2 h I_d\right) = -\mathbf{U}_{3/2} G_{3/2}, \end{array} \right. $$

where $h_i$ and $\mathbf{U}_i = (u_i, v_i)^t$ ($i \in \{1, 2\}$) are respectively the approximation of the layer thickness and the mean horizontal velocity in the $i^{th}$ layer. The mass exchanged from the second layer to the first is denoted $G_{3/2}$. Since the model contains four equations and six unknowns, some closure relations have to be added. First the layer thickness is assumed to be proportional to the total flow depth. In the current work we assume for simplicity an homogeneous vertical discretization, i.e. all the layer have the same thickness. In addition $\mathbf{U}_{3/2}$, the velocity at the interface, has to be defined. Here, it is approximated by a centered formula. The two closure relations read

$$ (4) \qquad h_1 = h_2 \quad \text{and} \quad \mathbf{U}_{3/2} = \frac{\mathbf{U}_1 + \mathbf{U}_2}{2}. $$

**Proposition 1.** For any smooth enough solution, the 2D bilayer model (3) with closure (4) is equivalent to the shallow water model with two velocities (1) with the following relations between the unknowns

$$ (5) \qquad (h, \bar{\mathbf{U}}, \hat{\mathbf{U}}) = \left( h_1 + h_2, \frac{\mathbf{U}_1 + \mathbf{U}_2}{2}, \frac{\mathbf{U}_2 - \mathbf{U}_1}{2} \right) $$

and conversely

$$ (h_1, h_2, \mathbf{U}_1, \mathbf{U}_2, \mathbf{U}_{3/2}, G_{3/2}) = \left( \frac{h}{2}, \frac{h}{2}, \bar{\mathbf{U}} - \hat{\mathbf{U}}, \bar{\mathbf{U}} + \hat{\mathbf{U}}, \bar{\mathbf{U}}, -\frac{1}{2} \nabla \cdot (h \hat{\mathbf{U}}) \right). $$

*Proof.* The computations are straightforward, we only detail the way to write the shallow water model with two velocities (1) from the bilayer model (3). Since $h_1 = h_2$ in (4), the definition of the water depth $h$ in (5) leads to $h_1 = h_2 = h/2$. Then the total mass conservation in (1) is obtained by adding the two local mass conservation laws in (3). In the same way, the momentum conservation for the mean velocity $\bar{\mathbf{U}}$ in (1) is obtained by summing the two local momentum conservations in (3) after some easy algebraic manipulations on the tensor products. The derivation of the last equation in (1) for the standard deviation $\hat{\mathbf{U}}$ requires a little more computation. One first subtracts the two local momentum equations in (3) to obtain

$$ (6) \qquad \partial_t (h \hat{\mathbf{U}}) + \nabla \cdot \left( h \left( \bar{\mathbf{U}} \otimes \hat{\mathbf{U}} + \hat{\mathbf{U}} \otimes \bar{\mathbf{U}} \right) \right) = - \mathbf{U}_{3/2} G_{3/2}. $$
But the difference between the two local mass equations in (3) leads to

$$2G_{3/2} = -\nabla \cdot (h\hat{\mathbf{U}}).$$

Developing the derivatives and performing some computations on the tensor products in Relation (6) then lead the result. □

Let us now consider the model introduced in [22], or, more precisely, its two dimensional extension presented in [11]. It reads

$$ (7) \qquad \left\{ \begin{array}{l} \partial_t h + \nabla \cdot (h \bar{\mathbf{U}}) = 0, \\ \partial_t (h \bar{\mathbf{U}}) + \nabla \cdot (h \bar{\mathbf{U}} \otimes \bar{\mathbf{U}} + h \mathbf{P} + \frac{g}{2} h^2 \mathbf{I}_d) = 0, \\ \partial_t \mathbf{P} + (\bar{\mathbf{U}} \cdot \nabla) \mathbf{P} + (\nabla \bar{\mathbf{U}}) \mathbf{P} + \mathbf{P} (\nabla \bar{\mathbf{U}})^T = 0, \end{array} \right. $$

where $\mathbf{P}$ is a 2 × 2 symmetric stress tensor that measures the distortion of the horizontal velocity profile in the vertical direction. The models (1) and (7) are both extension of the shallow water model but with different modeling assumptions. The model (7) corresponds to a closure of the momentum equations of the vertical-averaged Euler system, whereas Model (1) comes from a vertical piecewise constant discretization of the horizontal velocity in the same system. Note that, in 2D, there are clearly different since they do not have the same number of unknowns but in the 1D case, both models are equivalent for smooth solutions, see [22, eq. (2.25)-(2.27)] and (11).

## 1.2 Properties of the model

Here we mention three properties associated to the shallow water model with two velocities (1): energy conservation, rotational invariance, hyperbolity. The energy is defined by

$$E = \frac{h}{2} \left( \|\bar{\mathbf{U}}\|^2 + \|\hat{\mathbf{U}}\|^2 \right) + \frac{g}{2} h^2.$$

**Proposition 2.** Any smooth enough solution of (1) satisfies the energy conservation law

$$\partial_t E + \nabla \cdot \left( \left( \left( E + \frac{g}{2} h^2 \right) \mathbf{I}_d + h \hat{\mathbf{U}} \otimes \hat{\mathbf{U}} \right) \bar{\mathbf{U}} \right) = 0.$$

*Proof.* The proof relies on classical computations and is left to the reader. □

For the classical shallow water model, the mechanical energy is used to determine the admissible solution in the case of a shock, i.e. it acts as a mathematical entropy. To obtain the same admissible solution in the case where $\hat{\mathbf{U}} = 0$, the same argument is used, i.e. we will consider that a discontinuous solution is admissible if the energy is decreasing

$$ (8) \qquad \partial_t E + \nabla \cdot \left( \left( \left( E + \frac{g}{2} h^2 \right) \mathbf{I}_d + h \hat{\mathbf{U}} \otimes \hat{\mathbf{U}} \right) \bar{\mathbf{U}} \right) \le 0.$$
**Proposition 3.** System (1) is invariant by rotation.

*Proof.* The proof relies on classical rotational invariance properties of the derivative operators. □

Let us now consider the hyperbolicity of the shallow water model with two velocities (1). It can be written in quasi linear form

$$ \partial_t W + A(W)\partial_x W + B(W)\partial_y W = 0, $$

for a convenient choice of vector W and matrices A(W) and B(W). Since the system is invariant by rotation, it is sufficient to consider the direction of the x axis, thus the case where the unknowns are independent of the y-direction. In this case, System (1) reads, with the notations $\bar{U} = (\bar{u}, \bar{v})^t$ and $\hat{U} = (\hat{u}, \hat{v})^t$,

$$ (9) \qquad \left\{ \begin{aligned} \partial_t h &+ \partial_x (h\bar{u}) & = 0, \\ \partial_t (h\bar{u}) &+ \partial_x \left(h(\bar{u}^2 + \hat{u}^2) + \frac{g}{2}h^2\right) & = 0, \\ \partial_t \hat{u} &+ \partial_x (\hat{u}\bar{u}) & = 0, \\ \partial_t (h\bar{v}) &+ \partial_x (h(\bar{u}\bar{v} + \hat{u}\hat{v})) & = 0, \\ \partial_t \hat{v} &+ \hat{u}\partial_x \bar{v} + \bar{u}\partial_x \hat{v} & = 0, \end{aligned} \right. $$

and the initial data of the planar Riemann problem, given in the set of variables $V = (h, \bar{u}, \hat{u}, \bar{v}, \hat{v})^t$, is

$$ (10) \qquad V(0, x) = \begin{cases} V_L = (h_L, \bar{u}_L, \hat{u}_L, \bar{v}_L, \hat{v}_L)^t \in \mathbb{R}_+^* \times \mathbb{R}^4 & \text{if } x < 0, \\ V_R = (h_R, \bar{u}_R, \hat{u}_R, \bar{v}_R, \hat{v}_R)^t \in \mathbb{R}_+^* \times \mathbb{R}^4 & \text{if } x \ge 0. \end{cases} $$

Note that the quantities $(\bar{u}, \hat{u})$ can be seen as the components of the normal velocity while the quantities $(\bar{v}, \hat{v})$ can be seen as the components of the transverse velocity. More important, it appears that the three first equations of system (9) on $(h, \bar{u}, \hat{u})$, which correspond to the one-dimensional model, are independent of the two last ones on the transverse velocity $(\bar{v}, \hat{v})$. Thus we begin by analyzing the hyperbolicity and the Riemann problem for this simplified case.

## 2 The Riemann problem for the 1D model

The one dimensional shallow water model with two velocities reads

$$ (11) \qquad \left\{ \begin{aligned} \partial_t h &+ \partial_x (h\bar{u}) & = 0, \\ \partial_t (h\bar{u}) &+ \partial_x \left(h(\bar{u}^2 + \hat{u}^2) + \frac{g}{2}h^2\right) & = 0, \\ \partial_t \hat{u} &+ \partial_x (\hat{u}\bar{u}) & = 0, \end{aligned} \right. $$

and the initial data for the Riemann problem are given (in the set of variables $(h, \bar{u}, \hat{u})^t$) by

$$ (12) \qquad U(0, x) = \begin{cases} U_L = (h_L, \bar{u}_L, \hat{u}_L)^t \in \mathbb{R}_+^* \times \mathbb{R}^2 & \text{if } x < 0, \\ U_R = (h_R, \bar{u}_R, \hat{u}_R)^t \in \mathbb{R}_+^* \times \mathbb{R}^2 & \text{if } x \ge 0. \end{cases} $$
Contrary to the two dimensional case (10), the one dimensional shallow water model with two velocities (11) is conservative. The solution of the Riemann problem will then be much more simple to compute. The associated energy is defined by

$$ (13) \qquad E = \frac{g}{2}h^2 + \frac{h}{2}(\bar{u}^2 + \hat{u}^2), $$

and the admissible discontinuous solutions have to satisfy

$$ (14) \qquad \partial_t E + \partial_x \left( \left( \frac{\bar{u}^2 + 3\hat{u}^2}{2} + gh \right) \bar{u} \right) \le 0. $$

For smooth solutions, System (11) is equivalent to the full-Euler system for gas dynamics. More precisely, using the notation of [13, Chapter I.2] the full-Euler system can be written as

$$ (15) \qquad \begin{cases} \partial_t \rho + \partial_x (\rho u) = 0, \\ \partial_t (\rho u) + \partial_x (\rho u^2 + p) = 0, \\ \partial_t (\rho e) + \partial_x ((\rho e + p) u) = 0, \end{cases} $$

where $\rho$ denotes the density of the gas, $u$ its velocity and the pressure $p$ is defined as a function of internal energy $\varepsilon$ and density. Introducing the specific entropy $S$ such that $\rho^2 T dS = \rho^2 d\varepsilon + p d\rho$ with $T$ the temperature smooth solutions satisfy the relation

$$ \partial_t \eta + \partial_x (\eta u) = 0 $$

with $\eta = -\rho S$. Discontinuous solutions are admissible if they satisfy the entropy dissipation

$$ \partial_t \eta + \partial_x (\eta u) \le 0. $$

The analogy between the two systems is emphasized if we set

$$ \begin{pmatrix} \rho \\ \rho u \\ \rho e \end{pmatrix} = \begin{pmatrix} h \\ \frac{h\bar{u}}{E} \\ E \end{pmatrix} \quad \text{with} \quad p = \frac{g}{2}h^2 + h\hat{u}^2 \quad \text{and} \quad \eta = \hat{u}. $$

Thanks to this observation, the following hyperbolicity analysis is a corollary of the analysis of the full-Euler system except for shocks. More precisely, conversely to the full-Euler system, we look for solutions of System (11) where the mechanical energy $E$ is not conserved but acts as a mathematical entropy and is decreasing through an admissible shock (while the variable $\hat{u}$, corresponding to the entropy $\eta$ in Euler system, is conserved).

## 2.1 Hyperbolicity of the 1D model

The eigenstructure of the model can be studied with any set of variables. Thus, by analogy with the Euler system (15), we obtain easily the following result.
<figure><img src="image_2.png" /><figcaption>Figure 1: Solution of the 1D Riemann problem in the case of a $\lambda_L$-rarefaction and $\lambda_R$-shock.</figcaption></figure>

**Lemma 1.** Assuming $h > 0$, Model (11) is strictly hyperbolic. More precisely, the eigenvalues are given by

$$ \lambda_L = \bar{u} - \sqrt{gh + 3\hat{u}^2} < \lambda_* = \bar{u} < \lambda_R = \bar{u} + \sqrt{gh + 3\hat{u}^2}. $$

In addition, the $\lambda_L$-field and the $\lambda_R$-field are genuinely nonlinear, whereas the $\lambda_*$-field is linearly degenerate. A set of independent $\kappa$-Riemann invariants, denoted $\mathbb{I}_\kappa$ is given by

$$ \mathbb{I}_{\lambda_L} = \left\{ \bar{u} + \sqrt{gh + 3\hat{u}^2} + \frac{gh}{\sqrt{3\hat{u}}} \log \left( \sqrt{1 + \frac{3\hat{u}^2}{gh}} + \sqrt{\frac{3}{gh}\hat{u}} \right), \frac{\hat{u}}{h} \right\}, $$

$$ \mathbb{I}_{\lambda_*} = \left\{ \bar{u}, \frac{g}{2}h^2 + h\hat{u}^2 \right\}, $$

$$ \mathbb{I}_{\lambda_R} = \left\{ \bar{u} - \sqrt{gh + 3\hat{u}^2} - \frac{gh}{\sqrt{3\hat{u}}} \log \left( \sqrt{1 + \frac{3\hat{u}^2}{gh}} + \sqrt{\frac{3}{gh}\hat{u}} \right), \frac{\hat{u}}{h} \right\}. $$

## 2.2 Nonlinear waves in the 1D model

We now focus on the analysis of the genuinely nonlinear waves associated to the extreme eigenvalues. Let us first show that the dissipation of the mechanical energy is equivalent to the Lax entropy condition see [13, Chapter I.5]. Note that the following Proposition is not a direct application of [13, Chapter I.5, Theorem 5.3] since the mechanical energy $E$ (13) is not a convex function of the conservative state variables $W = (h, h\bar{u}, \hat{u})$.

**Proposition 4.** We denote by $\sigma_k$ the speed of the $\lambda_k$-shock. Assuming that the water depth $h$ is positive, the following properties are equivalent:

i) The mechanical energy $E$ (13) is decreasing through a shock, i.e. for $k \in \{L, R\}$

$$ -\sigma_k [E] + \left[ \left( \frac{\bar{u}^2 + 3\hat{u}^2}{2} + gh \right) h\bar{u} \right] < 0. $$
ii) the shock is compressive, i.e. for $k \in \{L, R\}$ we have

$$ (16) \qquad \begin{aligned} & -\sigma_k [h^2] + [h^2 \bar{u}] > 0 \\ \text{or} \quad & -\sigma_k [h \bar{u}^2] + [h \bar{u}^3] < 0 \\ \text{or} \quad & -\sigma_k [h \hat{u}^2] + [h \hat{u}^2 \bar{u}] > 0. \end{aligned} $$

iii) the Lax entropy condition is satisfied

$$ (17) \qquad \lambda_L (U_L) > \sigma_L > \lambda_L (U_{L*}) \quad \text{and} \quad \lambda_R (U_{R*}) > \sigma_R > \lambda_R (U_R) $$

where $U_L = (h_L, \bar{u}_L, \hat{u}_L)^t \in \mathbb{R}_+^* \times \mathbb{R}^2$ and $U_{L*} = (h_{L*}, \bar{u}_{L*}, \hat{u}_{L*})^t \in \mathbb{R}_+^* \times \mathbb{R}^2$ designate respectively the states at the left and at the right of the $\lambda_L$-shock and $U_{R*} = (h_{R*}, \bar{u}_{R*}, \hat{u}_{R*})^t \in \mathbb{R}_+^* \times \mathbb{R}^2$ and $U_R = (h_R, \bar{u}_R, \hat{u}_R)^t \in \mathbb{R}_+^* \times \mathbb{R}^2$ designate respectively the states at the left and at the right of the $\lambda_R$-shock, see Figure 1.

**Proof.** We focus on the $\lambda_L$-shock. Let us first prove that the three conditions (16) are equivalent. We write the problem in the framework of the shock, i.e. we set $\bar{w} = \bar{u} - \sigma_L$ and $\hat{w} = \hat{u}$. The Rankine-Hugoniot relations for Model (11) read

$$ \left\{ \begin{array}{l} \sigma_L [h] = [h \bar{u}], \\ \sigma_L [h \bar{u}] = [h (\bar{u}^2 + \hat{u}^2) + \frac{g}{2} h^2], \\ \sigma_L [\hat{u}] = [\hat{u} \bar{u}] \end{array} \right. $$

or equivalently in the framework of the shock

$$ (18) \qquad \left\{ \begin{array}{l} [h \bar{u}] = 0, \\ [h (\bar{u}^2 + \hat{w}^2) + \frac{g}{2} h^2] = 0, \\ [\hat{w} \bar{u}] = 0. \end{array} \right. $$

It obviously yields that $Q = h\bar{w}$ and $S = \hat{w}/h$ are constant through the shock. Fixing the left state $U_L$ we write

$$ (19) \qquad \bar{w}_{L*} = \frac{Q}{h_{L*}} \quad \text{and} \quad \hat{w}_{L*} = Sh_{L*}. $$

The function $h_{L*} \mapsto \bar{w}_{L*}^2$ is then decreasing while $h_{L*} \mapsto \hat{w}_{L*}^2$ is increasing. It follows that the three conditions

$$ Q[h] > 0 \quad \text{or} \quad Q[\bar{w}^2] < 0 \quad \text{or} \quad Q[\hat{w}^2] > 0 $$

are equivalent. These relations are nothing but the three relations (16) written in the framework of the shock.

Now we show the equivalence between properties i and ii. The left-hand side of the mechanical energy inequality reads

$$ \begin{aligned} -\sigma_L [E] + \left[ \left( \frac{\bar{u}^2 + 3\hat{u}^2}{2} + gh \right) h \bar{u} \right] &= Q \left[ \frac{\bar{w}^2 + 3\hat{w}^2}{2} + gh \right] \\ &= Q \left( g[h] + \frac{Q^2}{2} \left[ \frac{1}{h^2} \right] + \frac{3}{2} S^2 [h^2] \right). \end{aligned} $$
Using the second Rankine-Hugoniot condition (18) it yields

$$ (20) \qquad \frac{Q^2}{h_L h_{L*}} = g \frac{h_L + h_{L*}}{2} + S^2 (h_L^2 + h_L h_{L*} + h_{L*}^2) $$

and then the left-hand side of the mechanical energy inequality reads

$$ (21) \qquad -\sigma_L [E] + \left[ \left( \frac{\bar{u}^2 + 3\hat{u}^2}{2} + gh \right) h\bar{u} \right] = -\frac{Q[h]^3}{4h_L h_{L*}} (g + 2S^2 (h_L + h_{L*})) $$

It follows that the mechanical energy is decreasing through a shock if and only if relations (16) are satisfied.

Finally we show the equivalence between properties ii and iii. First we assume the Lax entropy condition iii holds, which reads in the framework of the $\lambda_L$-shock

$$ \bar{w}_{L*} - \sqrt{gh_{L*} + 3\hat{w}_{L*}^2} < 0 < \bar{w}_L - \sqrt{gh_L + 3\hat{w}_L^2} $$

From the right inequality, it follows that $\bar{w}_L > 0$ then $Q = h_L \bar{w}_L > 0$. Since $h_{L*} \mapsto \bar{w}_{L*} > 0$ is decreasing whereas $h_{L*} \mapsto \hat{w}_{L*}^2$ is increasing, it follows that $h_{L*} \mapsto (\lambda_L (U_{L*}) - \sigma_L)$ is decreasing and the condition $[\lambda_L - \sigma_L] < 0$ yields $[h] > 0$. Since $Q > 0$ at the $\lambda_L$-shock, property ii is satisfied. Conversely, we assume property ii holds, i.e. the $\lambda_L$-shock is compressible. By definition of a shock we have $[\lambda_L] < 0$ then $[\lambda_L - \sigma_L] < 0$ and since $h_{L*} \mapsto (\lambda_L (U_{L*}) - \sigma_L)$ is decreasing, we conclude that $[h] > 0$ thus $Q > 0$. Now using (20) we get

$$ \begin{aligned} \bar{w}_L - \sqrt{gh_L + 3\hat{w}_L^2} &= \frac{h_{L*}}{h_L} \left( \frac{Q^2}{h_L h_{L*}} - \frac{h_L}{h_{L*}} (gh_L + 3S^2 h_L^2) \right) \\ &= \frac{[h]}{h_L} \left( \frac{g}{2} (2h_L + h_{L*}) + S^2 (3h_L^2 + 2h_L h_{L*} + h_{L*}) \right) > 0. \end{aligned} $$

Similarly we have $\bar{w}_{L*} - \sqrt{gh_{L*} + 3\hat{w}_{L*}^2} < 0$. We conclude property iii holds. We proceed similarly with the $\lambda_R$-shock. $\square$

From now on, we designate by admissible a shock that satisfies the conditions of Proposition 4. We are now in position to determine the relation linking the states on each side of an admissible shock. Let us start with the $\lambda_L$-wave.

**Lemma 2.** The admissible right states $U_{L*} = (h_{L*}, \bar{u}_{L*}, \hat{u}_{L*})^t \in \mathbb{R}_+^* \times \mathbb{R}^2$ that can be reached from the left state $U_L = (h_L, \bar{u}_L, \hat{u}_L)^t \in \mathbb{R}_+^* \times \mathbb{R}^2$ through an admissible $\lambda_L$-shock verifies

$$ \begin{aligned} h_{L*} &> h_L, \\ \bar{u}_{L*} &= \bar{u}_L - (\mu_L^s(h_{L*}) - \mu_L^s(h_L)), \\ \hat{u}_{L*} &= \frac{h_{L*}}{h_L} \hat{u}_L, \end{aligned} $$

with $\mu_L^s(X) = \mu^s(X; U_L)$,

$$ (22) \qquad \mu^s(X; (h, \bar{u}, \hat{u})) = \frac{X-h}{X} c(X; (h, \bar{u}, \hat{u})), $$
$c_L(X) = c(X; U_L)$ and

$$ (23) \qquad c(X; (h, \bar{u}, \hat{u})) = \sqrt{\frac{X}{h} \left( \left(1 + \frac{X}{h}\right) \frac{gh}{2} + \left(1 + \frac{X}{h} + \left(\frac{X}{h}\right)^2\right) \hat{u}^2 \right)}, $$

The velocity of the shock is given by $\sigma_L = \bar{u}_L - c_L(h_{L*})$. Note that $h_{L*} \mapsto \bar{u}_{L*}$ is strictly decreasing and $h_{L*} \mapsto \hat{u}_{L*}^2$ is strictly increasing.

**Proof.** From Proposition 4, we deduce that $h_{L*} > h_L$, Relation (19) holds and $\bar{w}_L$ is positive. The value of $\hat{u}_{L*}$ is a direct consequence of Relation (19). The values of $\bar{u}_{L*}$ and $\sigma_L$ are consequences of the equality $\bar{w}_L = c_L(h_{L*})$. This latter relation is deduced from the second Rankine-Hugoniot relation (18), remember that $\bar{w}_L$ is positive. The function $h_{L*} \mapsto \hat{u}_{L*}$ is obviously increasing. It follows that the function $h_{L*} \mapsto c_L(h_{L*})$ is also increasing. Thus the function $h_{L*} \mapsto \bar{u}_{L*} = \bar{w}_{L*} + \sigma_L = \bar{w}_{L*} + \bar{u}_L - c_L(h_{L*})$ is strictly decreasing since the function $h_{L*} \mapsto \bar{w}_{L*}$ is decreasing, see Proposition 4. $\square$

Then we analyse the $\lambda_L$-rarefaction. Since the solution is continuous in this case, the following results can be deduced from the analysis of the rarefaction wave of the full-Euler system with a convenient pressure law, see [13, Chapter II.3]. Nevertheless we detail the resulting formulas since the results will be used to solve the two dimensional Riemann problem in the next section.

**Lemma 3.** Inside a $\lambda_L$-rarefaction, the water depth $h$ is strictly decreasing. More precisely, the left state $U_L = (h_L, \bar{u}_L, \hat{u}_L)^t \in \mathbb{R}_+^* \times \mathbb{R}^2$ and the right state $U_{L*} = (h_{L*}, \bar{u}_{L*}, \hat{u}_{L*})^t \in \mathbb{R}_+^* \times \mathbb{R}^2$ are linked by a $\lambda_L$-rarefaction if and only if

$$ h_{L*} \le h_L, \qquad \hat{u}_{L*} = \frac{h_{L*}}{h_L} \hat{u}_L \qquad \text{and} \qquad \bar{u}_{L*} = \bar{u}_L - (\mu_L^r(h_{L*}) - \mu_L^r(h_L)) $$

with $\mu_L^r(X) = \mu^r(X; U_L)$,

$$ (24) \qquad \mu^r(X; (h, \bar{u}, \hat{u})) = \sqrt{gX + \frac{3\hat{u}^2}{h^2}X^2 + \frac{gh}{\sqrt{3}\hat{u}}} \log \left( \sqrt{1 + \frac{3\hat{u}^2}{gh^2}X} + \frac{\hat{u}}{h}\sqrt{\frac{3X}{g}} \right). $$

In particular, $h_{L*} \mapsto \bar{u}_{L*}$ is strictly decreasing and $h_{L*} \mapsto \hat{u}_{L*}^2$ is strictly increasing. Moreover, the rarefaction can be parametrized by $\xi = x/t$. The solution $\xi \mapsto (h, \bar{u}, \hat{u})$ varies from $U_L$ to $U_{L*}$ in the set

$$ \xi \in \left[ \bar{u}_L - \sqrt{ghL + 3\hat{u}_L^2}, \bar{u}_{L*} - \sqrt{gh_{L*} + 3\hat{u}_{L*}^2} \right], $$

and the value at $\xi$ is uniquely defined by

$$ \begin{align} \bar{u}(\xi) - \sqrt{gh(\xi) + 3\hat{u}^2(\xi)} &= \xi, \tag{25} \\ \bar{u}(\xi) &= \bar{u}_L - (\mu_L^r(h(\xi)) - \mu_L^r(h_L)), \nonumber \\ \hat{u}(\xi) &= \frac{h(\xi)}{h_L} \hat{u}_L. \nonumber \end{align} $$
*Proof.* Since we deal with continuous solutions, the computations are obtained from the analogy with the Euler system (15). □

The two previous results can also be established for the $\lambda_R$-wave. We do not detail their proofs, as they are completely similar with those of Lemmas 2 and 3.

**Lemma 4.** The left states $U_{R*} = (h_{R*}, \bar{u}_{R*}, \hat{u}_{R*})^t$ that can be linked to the right state $U_R = (h_R, \bar{u}_R, \hat{u}_R)^t \in \mathbb{R}_+^* \times \mathbb{R}^2$ through a $\lambda_R$-wave, rarefaction or shock satisfying the Lax entropy condition (17), is defined by

$$ \bar{u}_{R*} = \bar{u}_R + \begin{cases} \mu_R^r(h_{R*}) - \mu_R^r(h_R) & \text{if } h_{R*} \le h_R, \quad (\text{rarefaction}) \\ \mu_R^s(h_{R*}) - \mu_R^s(h_R) & \text{if } h_{R*} > h_R, \quad (\text{shock}) \end{cases} \\ \hat{u}_{R*} = \frac{h_{R*}}{h_R} \hat{u}_R, $$

where $\mu_R^s(X) = \mu^s(X; U_R)$, $\mu_R^r(X) = \mu^r(X; U_R)$ are respectively defined by (22) and (24). In particular, $h_{R*} \mapsto \bar{u}_{R*}$ and $h_{R*} \mapsto \hat{u}_{R*}^2$ are strictly increasing. Eventually, the speed of the $\lambda_R$-shock is $\sigma_R = \bar{u}_R + c_R(h_{R*})$ where $c_R(X) = c(X; U_R)$ is defined in (23), and the rarefaction fan is parameterized for any

$$ \xi = \frac{x}{t} \in \left[ \bar{u}_{R*} + \sqrt{gh_{R*} + 3\hat{u}_{R*}^2}, \bar{u}_R + \sqrt{gh_R + 3\hat{u}_R^2} \right], $$

by

$$ \begin{aligned} \bar{u}(\xi) + \sqrt{gh(\xi) + 3\hat{u}^2(\xi)} &= \xi, \\ \bar{u}(\xi) &= \bar{u}_R + (\mu_R^r(h(\xi)) - \mu_R^r(h_R)), \\ \hat{u}(\xi) &= \frac{h(\xi)}{h_R} \hat{u}_R. \end{aligned} $$

## 2.3 Resolution of the Riemann problem

With the precise knowledge of the simple waves associated to the $\lambda_L$- and $\lambda_R$-waves, we are now in position to solve the one dimensional Riemann problem for (11) (12). Thanks to the analogy with the full-Euler system, we may follow closely the resolution of the Riemann problem for gas dynamics, see [13, Chapter II.3]. The four possible wave patterns are depicted in Figure 2.

**Theorem 1.** If the following condition is fulfilled,

$$ (26) \qquad \bar{u}_R - \bar{u}_L < \mu_L^r(h_L) + \mu_R^r(h_R) $$

there exists a unique selfsimilar solution $U \in (L^\infty(\mathbb{R}_+^* \times \mathbb{R}))^3$ to the 1D Riemann problem (11)-(12) satisfying the mechanical energy dissipation (14). In addition the water depth $h$ is strictly positive for all $(t, x) \in \mathbb{R}_+ \times \mathbb{R}$.

*Proof.* Following the classical approach for Riemann problems, the solution is sought with bounded variation, self-similar and composed by four constant
<figure><img src="image_3.png" /><figcaption>Figure 2: Projection of the Riemann invariants $I_{\lambda*}$ onto the two-dimensional $(h_{R*}, h_{L*})$ plane.</figcaption></figure>

states, $U_L$, $U_{L*}$, $U_{R*}$, $U_R$, separated by three waves, see Figure 1. By Lemmas 2 and 3, the state $U_{L*}$ on the right of the $\lambda_L$-wave is uniquely determined by the knowledge of $h_{L*}$. Similarly by Lemma 4, the state $U_{R*}$ on the left of the $\lambda_R$-wave is parametrized by $h_{R*}$. Thus the proof amounts to find a unique pair $h_{L*}$ and $h_{R*}$ is such that the Riemann invariants $I_{\lambda*}$ are preserved, i.e. $\bar{u}_{L*} = \bar{u}_{R*}$ and $p_{L*} := \frac{g}{2}h_{L*}^2 + h_{L*}\hat{u}_{L*}^2 = p_{R*} := \frac{g}{2}h_{R*}^2 + h_{R*}\hat{u}_{R*}^2$.

Since $h_{L*} \mapsto \bar{u}_{L*}$ is strictly decreasing and $h_{R*} \mapsto \bar{u}_{R*}$ is strictly increasing, the set $\{(h_{L*}, h_{R*}) \in (\mathbb{R}_+^*)^2 : \bar{u}_{L*}(h_{L*}) = \bar{u}_{R*}(h_{R*})\}$ is represented by a decreasing curve in the $(h_{L*}, h_{R*})$-plane. Similarly, since $h_{L*} \mapsto p_{L*}$ and $h_{R*} \mapsto p_{R*}$ are strictly increasing, the curve of the equation $p_{L*} = p_{R*}$ is strictly increasing in the $(h_{L*}, h_{R*})$-plane, see Figure 2.

We conclude that there exists at most one pair $(h_{R*}, h_{L*}) \in \mathbb{R}^2$ such that the Riemann invariants $I_{\lambda*}$ are constant. It remains to show that the solution is physically relevant, i.e. $(h_{R*}, h_{L*}) \in (\mathbb{R}_+^*)^2$. It is clear that the curve of the second Riemann invariant $p$ passes through $(0, 0)$ and tends to $+\infty$ as $h_{L*}$ goes to $+\infty$. Thus a solution without dry area ($h > 0$) exists if and only if the curve of the first Riemann invariant $\bar{u}$ crosses the axis $h_{R*} = 0$ in the upper half plane $h_{L*} > 0$. The function $h_{R*} \to \bar{u}_{L*}(0) - \bar{u}_{R*}(h_{R*})$ being decreasing and tending to $-\infty$ when $h_{R*}$ tends to $+\infty$, this is the case if and only if $\bar{u}_{L*}(0) - \bar{u}_{R*}(0)$ is strictly positive. It yields the condition (26).

Eventually, we check that the waves are strictly separated even for large initial data. It is trivial for rarefaction waves, and for shocks we obtained in Lemma 2 that $\sigma_L = \bar{u}_{L*} - \frac{h_L}{h_{L*}}c_L(h_{L*}) < \bar{u}_{L*}$ and similarly $\sigma_R > \bar{u}_{R*}$. $\square$

From now on, we denote by $\bar{u}_* := \bar{u}_{L*} = \bar{u}_{R*}$ and $p_* := p_{L*} = p_{R*}$ the mean velocity and the mean “pressure” (using the full-Euler terminology, here $p = \frac{g}{2}h^2 + h\hat{u}^2$) in the intermediate states $U_{L*}$ and $U_{R*}$.

**Remark 1.** A direct corollary of the monotonicity of the pressure $h_{k*} \mapsto p_{k*}$ is
<figure><img src="image_4.png" /><figcaption>Figure 3: Solution of the Riemann problem in the case of a $\lambda_L$-rarefaction and $\lambda_R$-shock without coalescence or resonance.</figcaption></figure>

that the pressure increases through the admissible shock, i.e. the relation

$$-\sigma_k [hp] + [hp \bar{u}] > 0$$

is equivalent to the conditions of Proposition 4.

# 3 The Riemann problem for the 2D model

Now that we have exhibited the solution of the 1D Riemann problem, we consider the planar Riemann problem for the 2D model (9)-(10). Constructing the solution of this problem is the main result of this work.

## 3.1 Hyperbolicity and waves patterns

Let us first study the eigenstructure of System (9). In the following, we will note $\lambda_L$, $\lambda_*$ and $\lambda_R$ the eigenvalues given in Lemma 1.

**Lemma 5.** If $\hat{u} \neq 0$, the 2D bilayer model (9) is strictly hyperbolic. More precisely, the eigenvalues are given by

$$\lambda_L < \gamma_L \leq \lambda_* \leq \gamma_R < \lambda_R$$

where we have noted $\gamma_L = \bar{u} - |\hat{u}|$ and $\gamma_R = \bar{u} + |\hat{u}|$. In addition, the $\lambda_L$-field and the $\lambda_R$-field are genuinely nonlinear whereas the $\gamma_L$-field, the $\lambda_*$-field and the $\gamma_R$-field are linearly degenerate. A set of independent $\kappa$-Riemann invariants, denoted $J_\kappa$, reads

$$J_{\lambda_L} = \mathbb{I}_{\lambda_L} \cup \left\{ \bar{v} - \frac{\hat{u}\hat{v}}{\frac{g}{2}h + \hat{u}^2} \bar{u}, \frac{\hat{u}\hat{v}}{\frac{g}{2}h + \hat{u}^2} \right\},$$

$$J_{\gamma_L} = \{h, \bar{u}, \hat{u}, \bar{v} + s\hat{v}\},$$

$$J_{\lambda_*} = \mathbb{I}_{\lambda_*} \cup \{\bar{v}, h\hat{u}\},$$

$$J_{\gamma_R} = \{h, \bar{u}, \hat{u}, \bar{v} - s\hat{v}\},$$

$$J_{\lambda_R} = \mathbb{I}_{\lambda_R} \cup \left\{ \bar{v} - \frac{\hat{u}\hat{v}}{\frac{g}{2}h + \hat{u}^2} \bar{u}, \frac{\hat{u}\hat{v}}{\frac{g}{2}h + \hat{u}^2} \right\},$$
with $s = \operatorname{sgn}(\hat{u})$ and where $\mathbb{I}_{\lambda_L}$, $\mathbb{I}_{\lambda_*}$ and $\mathbb{I}_{\lambda_R}$ are defined in Lemma 1.

Before proving the result, let us remark that the two new eigenvalues $\gamma_L, \gamma_R$ are clearly inherited from System (3) since they correspond to $u_1, u_2$. More precisely, we can easily check from formulas (5) that $\gamma_L = \min(u_1, u_2)$ and $\gamma_R = \max(u_1, u_2)$.

*Proof.* In order to identify the structure of the system, we write it under the quasi linear form

$$ \partial_t V + \mathcal{A}(V) \partial_x V = 0 $$

with $V = (h, \bar{u}, \hat{u}, \bar{v}, \hat{v})^t$, where the matrix $\mathcal{A}(V)$ is defined by

$$ (27) \qquad \mathcal{A}(V) = \begin{pmatrix} \bar{u} & h & 0 & 0 & 0 \\ g + \frac{\hat{u}^2}{h} & \bar{u} & 2\hat{u} & 0 & 0 \\ 0 & \hat{u} & \bar{u} & 0 & 0 \\ \frac{\hat{u}\hat{v}}{h} & 0 & \hat{v} & \bar{u} & \hat{u} \\ 0 & 0 & 0 & \hat{u} & \bar{u} \end{pmatrix}. $$

It is easy to check that $\lambda_L, \lambda_*, \lambda_R, \gamma_L$ and $\gamma_R$ are the eigenvalues of $\mathcal{A}(V)$, with the respective right eigenvectors

$$ \begin{aligned} \mathcal{R}_{\lambda_L} &= \begin{pmatrix} -h \\ \sqrt{gh + 3\hat{u}^2} \\ -\hat{u} \\ \frac{\hat{u}\hat{v}\sqrt{gh + 3\hat{u}^2}}{\frac{g}{2}h + \hat{u}^2} \\ -\frac{\hat{u}^2\hat{v}}{\frac{g}{2}h + \hat{u}^2} \end{pmatrix}, & \mathcal{R}_{\lambda_*} &= \begin{pmatrix} 2h\hat{u} \\ 0 \\ -(gh + \hat{u}^2) \\ 0 \\ (gh - \hat{u}^2)\frac{\hat{v}}{\hat{u}} \end{pmatrix}, & \mathcal{R}_{\lambda_R} &= \begin{pmatrix} h \\ \sqrt{gh + 3\hat{u}^2} \\ \hat{u} \\ \frac{\hat{u}\hat{v}\sqrt{gh + 3\hat{u}^2}}{\frac{g}{2}h + \hat{u}^2} \\ \frac{\hat{u}^2\hat{v}}{\frac{g}{2}h + \hat{u}^2} \end{pmatrix}, \\ \mathcal{R}_{\gamma_L} &= \begin{pmatrix} 0 \\ 0 \\ 0 \\ -s \\ 1 \end{pmatrix}, & \mathcal{R}_{\gamma_R} &= \begin{pmatrix} 0 \\ 0 \\ 0 \\ s \\ 1 \end{pmatrix}. \end{aligned} $$

Finally, we check that for each choice of eigenvalue $\kappa$ in $\{\lambda_L, \gamma_L, \lambda_*, \gamma_R, \lambda_R\}$, any component $\phi_\kappa \in \mathbb{J}_\kappa$ verifies $\nabla_V \phi_\kappa \cdot \kappa = 0$ and is thus a $\kappa$-Riemann invariant. $\square$

Let us now come back on the important point that was previously mentioned and motivated the study of the Riemann problem for the 1D model: the three first equations of (9) are independent of the two last ones since they do not involve $(\bar{v}, \hat{v})$, and they correspond to the 1D problem. It follows that the quantities $(h, \bar{u}, \hat{u})$ are the unique solution of the Riemann problem (11-12). Note that this result is only true if the definition of the admissible shocks are unchanged, i.e. the Lax entropy condition (17) and the dissipation of the 2D mechanical energy (8) are equivalent, which will be proved later on. This fact has two main consequences.
* The wave speeds are entirely determined by the 1D problem. This is in particular true for the two new internal $\gamma_L$- and $\gamma_R$-contact discontinuities, as their speeds are equal to the eigenvalues which are respectively $\bar{u}_* - |\hat{u}_{L*}|$ and $\bar{u}_* + |\hat{u}_{R*}|$.

* The quantities $(h, \bar{u}, \hat{u})$ are entirely determined by the 1D problem. In particular they do not jump through the $\gamma_L$- and $\gamma_R$-contact, and $\bar{u}$ does not jump across the $\lambda_*$-contact either. Let us assume the five waves associated to the five eigenvalues are separated and denote, from the left to the right, the intermediate states by $V_{L1}$, $V_{L2}$, $V_{R2}$ and $V_{R1}$, see Figure 3, $V_{kj} = (h_{kj}, \bar{u}_{kj}, \hat{u}_{kj}, \bar{v}_{kj}, \hat{v}_{kj})^t$ with $k \in \{L, R\}$ and $j \in \{1, 2\}$. We note as previously $(h_{L*}, \bar{u}_*, \hat{u}_{L*})$ and $(h_{R*}, \bar{u}_*, \hat{u}_{R*})$ the intermediate states of Theorem 1. It follows that we have

$$V_{L1} = (h_{L*}, \bar{u}_*, \hat{u}_{L*}, \bar{v}_{L1}, \hat{v}_{L1})^t,$$

$$V_{L2} = (h_{L*}, \bar{u}_*, \hat{u}_{L*}, \bar{v}_{L2}, \hat{v}_{L2})^t,$$

$$V_{R2} = (h_{R*}, \bar{u}_*, \hat{u}_{L*}, \bar{v}_{R2}, \hat{v}_{R2})^t,$$

$$V_{R1} = (h_{R*}, \bar{u}_*, \hat{u}_{L*}, \bar{v}_{R1}, \hat{v}_{R1})^t,$$

where the only unknown variables at this step are the eight components of the velocities in the direction $y$: $\bar{v}_{L1}, \hat{v}_{L1}; \bar{v}_{L2}, \hat{v}_{L2}; \bar{v}_{R2}, \hat{v}_{R2}; \bar{v}_{R1}, \hat{v}_{R1}$.

Let us now investigate the wave pattern of the solution of the Riemann problem.

* If $V_L$ and $V_R$ are close enough to a constant state $V = (h, \bar{u}, \hat{u}, \bar{v})$ with $\hat{u} \neq 0$, the strict hyperbolicity of Lemma 5 ensures that the five waves are separated and the solution contains five waves, as depicted on Figure 3. The solution of the Riemann problem is given in Section 3.2.

* When $V_L$ and $V_R$ are not close and if the $\lambda_R$-wave is a shock, it may happen that $\sigma_R$ is smaller than $\bar{u}_* + |\hat{u}_{R*}|$ the speed of the $\gamma_R$-contact, see Figure 4, left. It occurs only if $\hat{u}_R \neq 0$ and $h_{R*}$ is large enough, see Proposition 5 below. In that case, the shock and the contact coalesce and we have to define new jump conditions through this merged wave. This phenomena is not a resonance since the system is still strictly hyperbolic and the internal contact discontinuities will never cross the nonlinear waves to become the external waves. The analysis in this case is realized in Section 3.3 and numerical examples are given in Table 1 cases 5, 6 and 7.

* If $\hat{u}_R = 0$ it is clear that the $\gamma_R$-contact has the same speed as the $\lambda_*$-contact, namely $\bar{u}_*$. This is depicted on Figure 4, right. This case is a resonance in the classical sense since, when $\hat{u} = 0$, the matrix $\mathcal{A}$ defined in (27) is not diagonalizable. Note that two (if $\hat{u}_L = 0$ and $\hat{u}_R \neq 0$ or $\hat{u}_L \neq 0$ and $\hat{u}_R = 0$) or three (if $\hat{u}_L = \hat{u}_R = 0$) contact discontinuities can be in resonance. The analysis in this case is realized in Section 3.4 and numerical examples are given in Table 1 cases 3, 4 and 6.

The $\lambda_L$-rarefaction wave is always separated from the $\gamma_L$-contact. More precisely, the right velocity of the $\lambda_L$-rarefaction $\bar{u}_* - \sqrt{gh_{L*}} + 3\hat{u}_{L*}^2$ is clearly
<figure><img src="image_5.png" /><figcaption>Figure 4: Coalescence of a contact with a shock, case $\gamma_R \ge \sigma_R$ (left) and resonance of two contacts, case $\hat{u}_R = 0$ (right).</figcaption></figure>

smaller than the $\gamma_L$-contact velocity $\bar{u}_* - |\hat{u}_{L*}|$. The same holds for the $\gamma_R$-contact. Since the wave velocities only depend on the three first unknowns which are already estimated from the 1D analysis, we are in position to determine when the $\gamma_L$-contact wave coalesces with the $\lambda_L$-shock wave, see Figure 5.

**Proposition 5.** Let us define the following polynomial function for $X > h$

$$P(X; (h, \hat{u})) := \frac{1}{2} \left(1 + \frac{X}{h}\right) gh + \left(1 + \frac{X}{h} + \left(\frac{X}{h}\right)^2 - \left(\frac{X}{h}\right)^3\right) \hat{u}^2.$$

If $\hat{u} \ne 0$, $P$ has a unique real root larger than $h$ noted $\eta(h, \hat{u})$. The coalescence of the $\lambda_L$-shock and the $\gamma_L$-contact discontinuity occurs if and only if $h_{L*} \ge \eta_L := \eta(h_L, \hat{u}_L)$. Similarly, the coalescence of the $\lambda_R$-shock and the $\gamma_R$-contact discontinuity occurs if and only if $h_{R*} \ge \eta_R := \eta(h_R, \hat{u}_R)$.

*Proof.* We prove the result for the left waves as the proof for the right wave is exactly similar. The coalescence occurs if and only if $\bar{u}_* - |\hat{u}_{L*}| \le \sigma_L$, see Figure 3. Using Lemma 2, the condition reads

$$\bar{u}_{L*} - \frac{h_{L*}}{h_L} |\hat{u}_L| \le \bar{u}_L - c_L(h_{L*}) = \bar{u}_{L*} - \frac{h_L}{h_{L*}} c_L(h_{L*}).$$

After squaring and using (22), we obtain the equivalent condition $P(h_{L*}; (h_L, \hat{u}_L)) \le 0$. A trivial analysis of the polynomial function $P$ shows that it is increasing for $X \in [h, X_0]$ and decreasing for $X > X_0$ with

$$X_0 = \frac{\sqrt{1 + 3 \left(1 + \frac{gh}{2\hat{u}^2}\right) - 1}}{3} h.$$

In addition, the polynomial function $P$ is nonnegative when $X = h$ and tends to $-\infty$ when $X$ goes to $+\infty$. Thus $P$ has a unique real root larger than $h$. $\square$
<figure><img src="image_6.png" /><figcaption>Figure 5: Projection of the $I_*$-Riemann invariant curves in the $(h_{R*}, h_{L*})$ plane. The coalescence of a shock and a contact is indicated by parallel plain and dashed lines.</figcaption></figure>

In practice, it is not necessary to compute the critical water depth $\eta_L$. Once the intermediate water depth $h_{L*}$ is known, the coalescence occurs if and only if $P(h_{L*}; (h_L, \hat{u}_L)) \le 0$. Coalescences are not artificial and can appear in practice. More precisely on Figure 5, for any $(h_L, \hat{u}_L, h_R, \hat{u}_R)$ fixed, the curve $p_{L*} = p_{R*}$ is fixed. On the other hand the curve $\bar{u}_{L*} = \bar{u}_{R*}$ can be shifted to the top increasing $\bar{u}_L - \bar{u}_R$. Since the two curves are strictly monotonous, we deduce that there is no upper bound to the intermediate water depths $h_{L*}$ and $h_{R*}$. Numerical examples are given in Table 1 cases 5, 6 and 7.

## 3.2 Riemann problem without coalescence or resonance

We look for the solution of the Riemann problem (9) with (10) in the case without coalescence or resonance, i.e. if $\hat{u}_L \hat{u}_R \ne 0$, $h_{L*} < \eta_L$ and $h_{R*} < \eta_R$. Let us start by the analysis of the simple waves for Model (9) and more precisely with the rarefaction waves.

**Lemma 6.** The right states $V_{L2} = (h_{L*}, \bar{u}_*, \hat{u}_{L*}, \bar{v}_{L2}, \hat{v}_{L2})^t \in \mathbb{R}_+^* \times \mathbb{R}^4$ can be linked to the left state $V_L = (h_L, \bar{u}_L, \hat{u}_L, \bar{v}_L, \hat{v}_L)^t \in \mathbb{R}_+^* \times \mathbb{R}^4$ with a $\lambda_L$-rarefaction and a $\gamma_L$-contact discontinuity if and only if $(h, \bar{u}, \hat{u})$ satisfies the Lemma 3 and the following conditions hold

$$ \bar{v}_{L1} = \bar{v}_L - \frac{\hat{u}_L (\bar{u}_L - \bar{u}_*)}{\frac{g}{2} h_L + \hat{u}_L^2} \hat{v}_L \quad \text{and} \quad \hat{v}_{L1} = \frac{\frac{g}{2} h_L^2 + h_{L*} \hat{u}_L^2}{\frac{g}{2} h_L^2 + h_L \hat{u}_L^2} \hat{v}_L $$

and the rarefaction fan is parameterized
for any $\xi = \frac{x}{t} \in [\bar{u}_L - \sqrt{gh_L + 3\hat{u}_L^2}, \bar{u}_* - \sqrt{gh_{L*} + 3\hat{u}_{L*}^2}]$ by (25) and

$$ (\bar{v}, \hat{v})(\xi) = \left( \bar{v}_L - \frac{\hat{u}_L (\mu_L^r (h(\xi)) - \mu_L^r (h_L))}{\frac{g}{2} h_L + \hat{u}_L^2} \hat{v}_L, \frac{\frac{g}{2} h_L^2 + h(\xi) \hat{u}_L^2}{\frac{g}{2} h_L^2 + h_L \hat{u}_L^2} \hat{v}_L \right). $$

**Proof.** The part of the result concerning $(h, \bar{u}, \hat{u})$ has already been proven in Lemma 3. We recall that the $\gamma_L$-contact discontinuity is necessarily separated from the $\lambda_L$-rarefaction. The expressions of the mean velocity $\bar{v}_{L1}$ and of the standard deviation $\hat{v}_{L1}$ inside the rarefaction are directly obtained from the third and the fourth $\lambda_L$-Riemann invariant. $\square$

The case where the $\lambda_L$-wave is a shock is much more difficult because of the nonconservative products of the last equation of (9). Through the shock waves, the products are really nonconservative, i.e. the standard deviation $\hat{v}$ jumps simultaneously with the mean velocity $\bar{u}$ and the mean velocity $\bar{v}$ jumps simultaneously with the standard deviation $\hat{u}$, otherwise the Rankine-Hugoniot jump condition of the fourth equation of (9) may not be satisfied. Note that the $\lambda_L$-shock may also coalesce with the $\gamma_L$-contact but it is not treated in this section. We refer to Section 3.3 in this case. We are able to propose a definition of the nonconservative products based on a regularization of the water depth inside the shock, see [1, 15]. In the case of an isolated $\lambda_L$-shock wave, the water depth is given by

$$ h(t, x) = h_L + (h_{L*} - h_L) H(x - \sigma_L t) \quad \text{with} \quad H(x) = \begin{cases} 0 & \text{if } x \le 0, \\ 1 & \text{if } x > 0. \end{cases} $$

In order to define the nonconservative products we introduce a regularization of the Heaviside function $H_\epsilon(z)$ defined as follow. Let the regularization parameter $\epsilon > 0$ be fixed, and the regularization satisfies

$$ (28) \quad H_{\epsilon} \in C^{0}(\mathbb{R}) \cap C^{1}([0, \epsilon]), \quad H_{\epsilon} \text{ monotonic and } H_{\epsilon}(x) = \begin{cases} 0 & \text{if } x \le 0, \\ 1 & \text{if } x \ge \epsilon. \end{cases} $$

We defined the regularized water depth $h_\epsilon(z)$ as

$$ (29) \quad h_{\epsilon}(z) = h_{L} + (h_{L*} - h_{L}) H_{\epsilon}(z) \quad \text{with} \quad z = x - \sigma t. $$

We now look for a regularized profiles of the velocities $(\bar{u}_\epsilon, \hat{u}_\epsilon, \bar{v}_\epsilon, \hat{v}_\epsilon)$, also only depending on $z$, such that $(h_\epsilon, \bar{u}_\epsilon, \hat{u}_\epsilon, \bar{v}_\epsilon, \hat{v}_\epsilon)(0) = (h_L, \bar{u}_L, \hat{u}_L, \bar{v}_L, \hat{v}_L)$ and $(h_\epsilon, \bar{u}_\epsilon, \hat{u}_\epsilon, \bar{v}_\epsilon, \hat{v}_\epsilon)$ is a solution of (9). It appears that the states in $z = 0$ and $z = \epsilon$ are always linked by the same relations independently of the choice of the regularized Heaviside function $H_\epsilon(z)$.

**Lemma 7.** Let the left state $(h_L, \bar{u}_L, \hat{u}_L)$ and the right state $(h_{L*}, \bar{u}_*, \hat{u}_{L*})$ be two states linked by a $\lambda_L$-shock without coalescence, i.e. $h_L < h_{L*} < \eta_L$ with $\sigma_L$ the speed of this shock. Let $h_\epsilon$ defined by (29). Consider $V_\epsilon = (h_\epsilon, \bar{u}_\epsilon, \hat{u}_\epsilon, \bar{v}_\epsilon, \hat{v}_\epsilon)$ a piecewise $C^1$-solution of (9), depending only on $z = x - \sigma_L t$, such that $V_\epsilon(0) =$
$(h_L, \bar{u}_L, \hat{u}_L, \bar{v}_L, \hat{v}_L)$. Then $(h_\epsilon, \bar{u}_\epsilon, \hat{u}_\epsilon)(\epsilon) = (h_{L*}, \bar{u}_{L*}, \hat{u}_{L*})$ and $(\bar{v}_\epsilon, \hat{v}_\epsilon)(\epsilon)$ belong to the following set $\mathbb{V}_L$, independent of $\epsilon$:

$$ \mathbb{V}_L = \left\{ \begin{array}{l} (\bar{v}_{L*}, \hat{v}_{L*}) \in \mathbb{R}^2 \text{ such that} \\ \bullet h_{L*}((\bar{u}_{L*} - \sigma_L)\bar{v}_{L*} + \hat{u}_{L*}\hat{v}_{L*}) = h_L((\bar{u}_L - \sigma_L)\bar{v}_L + \hat{u}_L\hat{v}_L) \\ \bullet h_{L*}\hat{v}_{L*}\sqrt{(\bar{u}_{L*} - \sigma_L)^2 - \hat{u}_{L*}^2} = h_L\hat{v}_L\sqrt{(\bar{u}_L - \sigma_L)^2 - \hat{u}_L^2} \end{array} \right\}. $$

Conversely, if $(\bar{v}_{L*}, \hat{v}_{L*})$ belong to $\mathbb{V}_L$, there exists a unique piecewise $C^1$-solution of (9), depending only on $z = x - \sigma t$, that takes value $(h_{L*}, \bar{u}_{L*}, \hat{u}_{L*}, \bar{v}_{L*}, \hat{v}_{L*})$ at $z = \epsilon$.

*Proof.* As it was done in Lemma 2, we consider the velocities in the shock framework setting $(\bar{w}_L, \bar{w}_{L*}) = (\bar{u}_L, \bar{u}_{L*}) - \sigma_L$ and $(\hat{w}_L, \hat{w}_{L*}) = (\hat{u}_L, \hat{u}_{L*})$. The first and third equations of (9) yield to $h_\epsilon(z)\bar{w}_\epsilon(z) = h_L\bar{w}_L$ and $\hat{w}_\epsilon(z)\bar{w}_\epsilon(z) = \hat{w}_L\bar{w}_L$. It is clear in the smooth part of the solution and also a point of discontinuity by the Rankine–Hugoniot relations. Thus $\bar{w}_\epsilon$ and $\hat{w}_\epsilon$ verify

$$ \bar{w}_\epsilon(z) = \frac{h_L}{h_\epsilon(z)} \bar{w}_L, \qquad \hat{w}_\epsilon(t, z) = \frac{h_\epsilon(z)}{h_L} \hat{w}_L, $$

and are $C^1$-regular on the whole interval $[0, \epsilon]$ since $h_\epsilon > 0$ is $C^1$-regular. We recall that $\bar{w}_L$ is positive see Lemma 2, and thus that $\bar{w}_\epsilon$ is positive. Using the second equation of (9), it yields that $\frac{g}{2}h_\epsilon^2 + h_\epsilon(\bar{w}_\epsilon^2 + \hat{w}_\epsilon^2) = \frac{g}{2}h_L^2 + h_L(\bar{w}_L + \hat{w}_L)$, then we obtain that $(h_\epsilon, \bar{u}_\epsilon, \hat{u}_\epsilon)(\epsilon) = (h_{L*}, \bar{u}_{L*}, \hat{u}_{L*})$ exactly as in Lemma 2. The fourth equation of (9) yields that $h_\epsilon(\bar{w}_\epsilon\bar{v}_\epsilon + \hat{w}_\epsilon\hat{v}_\epsilon) = h_L(\bar{w}_L\bar{v}_L + \hat{w}_L\hat{v}_L)$ which is the first relation of $\mathbb{V}_L$.

Now let us focus on a point $z_0$ in $[0, \epsilon]$ such that $(\bar{v}_\epsilon, \hat{v}_\epsilon)$ is discontinuous. Thanks to the regularization $h_\epsilon$, $\bar{w}_\epsilon$ and $\hat{w}_\epsilon$ are continuous then the two last equations of (9) at $z_0$ reads

$$ (30) \qquad \begin{cases} [h_\epsilon(\bar{w}_\epsilon\bar{v}_\epsilon + \hat{w}_\epsilon\hat{v}_\epsilon)] = 0, \\ \hat{w}_\epsilon[\bar{v}_\epsilon] + \bar{w}_\epsilon[\hat{v}_\epsilon] = 0, \end{cases} \iff \begin{cases} \bar{w}_\epsilon(z_0)[\bar{v}_\epsilon] + \hat{w}_\epsilon(z_0)[\hat{v}_\epsilon] = 0, \\ \bar{w}_\epsilon(z_0)[\hat{v}_\epsilon] + \hat{w}_\epsilon(z_0)[\bar{v}_\epsilon] = 0. \end{cases} $$

The only solution is $[\bar{v}_\epsilon] = [\hat{v}_\epsilon] = 0$, unless $|\hat{w}_\epsilon(z_0)| = \bar{w}_\epsilon(z_0)$, i.e. $\sigma_L = \bar{u}_\epsilon(z_0) - |\hat{u}_\epsilon(z_0)|$. The study of Proposition 5 shows that it happens only if $h_\epsilon(z_0) = \eta_L$. Since $h_{L*} < \eta_L$ it yields that $(\bar{v}_\epsilon, \hat{v}_\epsilon)$ are $C^1$-regular. Using the fact that $h_\epsilon\bar{w}_\epsilon$ is a nonnegative constant, the fourth and fifth equations on (9) yield

$$ (31) \qquad \begin{cases} (h_\epsilon(\bar{w}_\epsilon\bar{v}_\epsilon + \hat{w}_\epsilon\hat{v}_\epsilon))' = 0, \\ \hat{w}_\epsilon\bar{v}_\epsilon' + \bar{w}_\epsilon\hat{v}_\epsilon' = 0, \end{cases} \iff \begin{cases} \left(\bar{v}_\epsilon + \frac{\hat{w}_\epsilon}{\bar{w}_\epsilon}\hat{v}_\epsilon\right)' = 0, \\ \frac{\hat{w}_\epsilon}{\bar{w}_\epsilon}\bar{v}_\epsilon' + \hat{v}_\epsilon' = 0. \end{cases} $$

Combining those two ODEs we obtain

$$ \left(1 - \left(\frac{\hat{w}_\epsilon}{\bar{w}_\epsilon}\right)^2\right) \hat{v}_\epsilon' - \frac{\hat{w}_\epsilon}{\bar{w}_\epsilon} \left(\frac{\hat{w}_\epsilon}{\bar{w}_\epsilon}\right)' \hat{v}_\epsilon = 0. $$
The hypothesis $h_{L*} < \eta_L$ implies that $|\hat{w}_\epsilon|$ remains smaller than $\bar{w}_\epsilon$. It follows that

$$ \hat{v}_\epsilon(z) = \hat{v}_L \sqrt{\frac{\bar{w}_\epsilon(z)^2}{\bar{w}_L^2} \frac{\bar{w}_L^2 - \hat{w}_L^2}{\bar{w}_\epsilon(z)^2 - \hat{w}_\epsilon(z)^2}} = \hat{v}_L \frac{h_L}{h_\epsilon(z)} \sqrt{\frac{\bar{w}_L^2 - \hat{w}_L^2}{\bar{w}_\epsilon(z)^2 - \hat{w}_\epsilon(z)^2}} $$

and leads to the second condition of $\mathbb{V}_L$. The converse is clear in that case. $\square$

Remark that applying the same regularization strategy in the case of the $\gamma_k$-contact discontinuity ($k \in \{L, R\}$), it corresponds to the case with a discontinuity point $z_0$ since $\bar{w}_{L*} = \bar{u}_* - \gamma_L = |\hat{u}_{L*}| = |\hat{w}_{L*}|$. We obtain from (30) that $\bar{v}_\epsilon + s_L \hat{v}_\epsilon = \bar{v}_L + s_L \hat{v}_L$ with $s_L = \operatorname{sgn}(\hat{u}_L) = \operatorname{sgn}(\hat{u}_{L*})$, which corresponds to the fourth $\mathbb{J}_{\gamma_L}$-Riemann invariant. In the case of the $\lambda_*$-contact discontinuity, we have $\bar{w}_\epsilon = 0$ and the regularized problem simply reads

$$ \begin{cases} (h_\epsilon \hat{w}_\epsilon \hat{v}_\epsilon)' = 0, \\ \hat{w}_\epsilon \hat{v}_\epsilon' = 0, \end{cases} $$

and we easily check that, when $\hat{w} \neq 0$, the relation verified by the regularized solution on $z = 0$ and $z = \epsilon$ are nothing but the $\mathbb{J}_{\lambda_*}$-Riemann invariants.

The regularization procedure of the shocks yields to relations between the left state and the right state that are independent of the choice of $\epsilon$ and of the regularization. This encourage us to define the nonconservative product as follow.

**Definition 1.** Through a $\lambda_k$-shock without coalescence ($k \in \{L, R\}$), i.e. $h_k < h_{k*} < \eta_k$, the following jump conditions hold:

$$ \left\{ \begin{array}{l} \sigma_k[h] = [h\bar{u}], \\ \sigma_k[h\bar{u}] = [h(\bar{u}^2 + \hat{u}^2) + \frac{g}{2}h^2], \\ \sigma_k[\bar{u}] = [\hat{u}\bar{u}], \\ \sigma_k[h\bar{v}] = [h(\bar{u}\bar{v} + \hat{u}\hat{v})], \\ \left[h\sqrt{(\bar{u} - \sigma_k)^2 - \hat{u}^2}\right] \bar{v} = 0. \end{array} \right. $$

with $\sigma_k$ the speed of the shock.

The natural question to address is once the jump conditions are defined, does the 2D mechanical energy (acting as a mathematical entropy) decrease through a $\lambda_k$-shock without coalescence (8). More precisely, does the Lax condition Proposition 4.ii define the admissible shock. To answer to this point, we highlight the following property.

**Proposition 6.** We denote by $\sigma_k$ the speed of the $\lambda_k$-shock. Assuming that the water depth $h$ is positive and there is no coalescence, i.e. $h_k < h_{k*} < \eta_k$, the following properties are equivalent:
i) The mechanical energy $E = \frac{g}{2}h^2 + \frac{h}{2}(\bar{u}^2 + \hat{u}^2 + \bar{v}^2 + \hat{v}^2)$ is decreasing through a shock, i.e. for $k \in \{L, R\}$

$$-\sigma_k [E] + \left[ \left( \frac{\bar{u}^2 + 3\hat{u}^2 + \bar{v}^2 + \hat{v}^2}{2} + gh \right) h\bar{u} + h\hat{v}\bar{v}\hat{u} \right] < 0.$$

ii) the shock is compressive, i.e. for $k \in \{L, R\}$ the inequalities of Proposition 4.ii hold and

$$-\sigma_k [h\hat{v}^2] + [h\hat{v}^2\bar{u}] > 0.$$

iii) the Lax entropy condition Proposition 4.iii is satisfied.

*Proof.* Let us first study the variation of the velocities $\bar{v}$ and $\hat{v}$ through a shock. From the last jump condition of Definition 1, we have

$$[\sqrt{Q^2 - S^2 h^4} \hat{v}] = 0$$

Where we have used the constants through the shock $Q = h\bar{w}$ and $S = \frac{\hat{w}}{h}$. It follows that $h_{L*} \mapsto \hat{v}_{L*}^2$ is increasing and using Proposition 4.ii we get that ii and iii are equivalent.

Let us now focus on the energy jump relation. Using the fourth jump relation of Definition 1 we set the constant $J = h(\bar{w}\bar{v} + \hat{w}\hat{v})$ and we get the relation

$$\begin{array}{l} J^2 = (Q\bar{v} + Sh^2\hat{v})^2 = Q^2\bar{v}^2 + S^2h^4\hat{v}^2 + 2QSh^2\bar{v}\hat{v} \\ \text{or} \\ v^2 = \frac{J^2 - S^2h^4\hat{v}^2 - 2QSh^2\bar{v}\hat{v}}{Q^2}. \end{array}$$

We study the jump of the energy flux coming from the 2D variable in the framework of the shock, i.e.

$$(32) \quad \left[ \frac{\bar{v}^2 + \hat{v}^2}{2} h\bar{w} + h\hat{v}\bar{v}\hat{w} \right] = \left[ \frac{J^2 - S^2h^4\hat{v}^2 - 2QSh^2\bar{v}\hat{v}}{2Q} + \frac{Q}{2}\hat{v}^2 + Sh^2\bar{v}\hat{v} \right] \\ = \frac{1}{2Q} [(Q^2 - S^2h^4)\hat{v}^2].$$

Using the last jump condition of Definition 1, we conclude that this term vanishes. Using Proposition 4, we conclude that i and iii are equivalent. $\square$

With Definition 1, we are in position to complete the description of the $\sigma_L$-shock wave in the case without coalescence.

**Lemma 8.** Assuming the jump conditions of Definition 1, the admissible right states $V_{L1} = (h_{L*}, \bar{u}_*, \hat{u}_{L*}, \bar{v}_{L1}, \hat{v}_{L1})^t$ that can be reached from the left state $V_L = (h_L, \bar{u}_L, \hat{u}_L, \bar{v}_L, \hat{v}_L)^t \in \mathbb{R}_+^4 \times \mathbb{R}^4$ with a $\lambda_L$-shock without coalescence, i.e. $h_L < h_{L*} < \eta_L$, satisfy the condition of Lemma 2 and

$$\begin{align*} \bar{v}_{L1} &= \bar{v}_L - \left( \frac{h_{L*}}{h_L} \sqrt{\frac{(\bar{u}_L - \sigma_L)^2 - \hat{u}_L^2}{(\bar{u}_* - \sigma_L)^2 - \hat{u}_{L*}^2}} - 1 \right) \frac{\hat{u}_L}{\bar{u}_L - \sigma_L} \hat{v}_L, \\ \hat{v}_{L1} &= \frac{h_L}{h_{L*}} \sqrt{\frac{(\bar{u}_L - \sigma_L)^2 - \hat{u}_L^2}{(\bar{u}_* - \sigma_L)^2 - \hat{u}_{L*}^2}} \hat{v}_L. \end{align*}$$
*Proof.* Using the Definition 1, the fifth condition directly gives the relation between $\hat{v}_{L1}$ and $\hat{v}_L$. Injecting this result in the fourth jump condition of Definition 1 and using that $h\bar{w}$ and $\frac{\hat{w}}{h}$ are constant, it yields the result. □

The previous results can as well be established for the right wave.

**Lemma 9.** Assuming the jump conditions of Definition 1, the admissible left states $V_{R1} = (h_{R*}, \bar{u}_*, \hat{u}_{R*}, \bar{v}_{R1}, \hat{v}_{R1})^t$ that can be reached from the right state $V_R = (h_R, \bar{u}_R, \hat{u}_R, \bar{v}_R, \hat{v}_R)^t \in \mathbb{R}_+^4 \times \mathbb{R}^4$ by a $\lambda_R$-wave (rarefaction or shock without coalescence) satisfying is characterized by

* If $h_{R*} < h_R$, the $\lambda_R$-wave is a rarefaction parametrized by $\xi = \frac{x}{t}$. The evolution of $(h, \bar{u}, \hat{u})(\xi)$ on $\xi \in [\bar{u}_{R*} + \sqrt{gh_{R*} + 3\hat{u}_{R*}^2}, \bar{u}_R + \sqrt{gh_R + 3\hat{u}_R^2}]$ is given by Lemma 4, and

$$ \begin{aligned} \bar{v}(\xi) &= \bar{v}_R - \frac{\hat{u}_R (\mu_R^r (h(\xi)) - \mu_R^r (h_R))}{\frac{g}{2} h_R + \hat{u}_R^2} \hat{v}_R, \\ \hat{v}(\xi) &= \frac{\frac{g}{2} h_R^2 + h(\xi) \hat{u}_R^2}{\frac{g}{2} h_R^2 + h_R \hat{u}_R^2} \hat{v}_R. \end{aligned} $$

* If $h_R < h_{R*} < \eta_R$, the intermediate state $V_{R1} = (h_{R*}, \bar{u}_*, \hat{u}_{R*}, \bar{v}_{R1}, \hat{v}_{R1})^t$ between the $\lambda_R$-shock and the $\gamma_R$-contact discontinuity satisfies Lemma 4 and

$$ \begin{aligned} \bar{v}_{R1} &= \bar{v}_R - \left( \frac{h_{R*}}{h_R} \sqrt{\frac{(\bar{u}_R - \sigma_R)^2 - \hat{u}_R^2}{(\bar{u}_* - \sigma_R)^2 - \hat{u}_{R*}^2}} - 1 \right) \frac{\hat{u}_R}{\bar{u}_R - \sigma_R} \hat{v}_R, \\ \hat{v}_{R1} &= \frac{h_R}{h_{R*}} \sqrt{\frac{(\bar{u}_R - \sigma_R)^2 - \hat{u}_R^2}{(\bar{u}_* - \sigma_R)^2 - \hat{u}_{R*}^2}} \hat{v}_R. \end{aligned} $$

*Proof.* The proof is similar to Lemma 6 and Lemma 8. □

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

We are now in position to solve the Riemann problem (9) with (10) in the case without coalescence or resonance. The following theorem will be improved in Section 3.3 and Section 3.4 to take into account coalescence and resonance.

**Theorem 2.** Assume the initial condition of the Riemann problem (10) satisfies the condition (26) and there is no-coalescence ($h_{L*} < \eta_L$ and $h_{R*} < \eta_R$) and no resonance ($\hat{u}_L \hat{u}_R \neq 0$). Then there exists a unique self-similar solution $V \in (L^\infty(\mathbb{R}_+^4 \times \mathbb{R}))^5$ of the Riemann problem (9) with (10), and the water depth is strictly positive.

*Proof.* The structure of the solution is depicted on Figure 3. The three first unknowns $(h, \bar{u}, \hat{u})$ are solution to the 1D-Riemann problem (1)-(12). In particular, under the condition (26) these three unknowns are unique and the water
depth is strictly positive. In addition, $(h, \bar{u}, \hat{u})$ do not vary through the second and fourth waves. Thus, the wave pattern contains five waves separated six constant states $V_L, V_{L1}, V_{L2}, V_{R2}, V_{R1}$ and $V_R$. The last two unknowns $(\bar{v}, \hat{v})$ in each state remain to be determined.

According to Lemma 6 and by Lemma 8, the state $V_{L1}$ is entirely determined by the knowledge of $V_L$ and the three first unknowns $(h_{L*}, \bar{u}_*, \hat{u}_{L*})$. Similarly, according to Lemma 9, the state $V_{R1}$ is entirely determined by the knowledge of $V_R$ and the three first unknowns $(h_{R*}, \bar{u}_*, \hat{u}_{R*})$. Thus the Riemann problem (9) is reduced to a Riemann problem of three contact discontinuities ($\gamma_L$-contact $\lambda_*$-contact and $\gamma_R$-contact discontinuities). Using the fourth $\mathbb{J}_{\gamma_L}$-Riemann invariant (resp. $\mathbb{J}_{\gamma_R}$-Riemann invariant), the relation between $V_{L1}$ and $V_{L2}$ (resp. $V_{R1}$ and $V_{R2}$) are of the form

$$ (33) \qquad \begin{aligned} \bar{v}_{L2} &= \bar{v}_{L1} + s_L (b_L \hat{v}_{L1} - a_L \hat{v}_{L2}), \\ \bar{v}_{R2} &= \bar{v}_{R1} - s_R (b_R \hat{v}_{R1} - a_R \hat{v}_{R2}). \end{aligned} $$

with $s_k = \operatorname{sgn}(\hat{u}_k) = \operatorname{sgn}(\hat{u}_{k*})$ for $k \in \{L, R\}$ and $a_L = b_L = a_R = b_R = 1$ in this case. The parameters $a_L, b_L, a_R$ and $b_R$ are introduced for homogeneity with the case with coalescence, see Section 3.3.

Finally, we focus on the $\lambda_*$-contact. The two last $\mathbb{J}_{\lambda_*}$-Riemann invariants, namely the continuity of $\bar{v}$ and of $h\hat{u}\hat{v}$, together with (33), leads to a linear system of 4 unknowns and 4 equations. Replacing $\bar{v}_{L2}$ and $\bar{v}_{R2}$ thanks to (33), $\hat{v}_{L2}$ and $\hat{v}_{R2}$ are solution of the following linear system

$$ (34) \qquad \begin{pmatrix} s_L a_L & s_R a_R \\ -h_{L*} \hat{u}_{L*} & h_{R*} \hat{u}_{R*} \end{pmatrix} \begin{pmatrix} \hat{v}_{L2} \\ \hat{v}_{R2} \end{pmatrix} = \begin{pmatrix} (\bar{v}_{L1} + s_L b_L \hat{v}_{L1}) - (\bar{v}_{R1} - s_R b_R \hat{v}_{R1}) \\ 0 \end{pmatrix} $$

whose determinant reads (we recall that $\hat{u}_L \hat{u}_{L*} > 0$ and $\hat{u}_R \hat{u}_{R*} > 0$)

$$ s_L a_L h_{R*} \hat{u}_{R*} + s_R a_R h_{L*} \hat{u}_{L*} = s_L s_R (a_L h_{L*} |\hat{u}_{L*}| + a_R h_{R*} |\hat{u}_{R*}|) \neq 0. $$

We conclude the uniqueness of $(\hat{v}_{L2}, \hat{v}_{R2})$ then of $(\bar{v}_{L2}, \bar{v}_{R2})$ by (33). $\square$

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

### 3.3 Coalescence of the shock and the contact discontinuity

We now focus on the case with coalescence of the $\lambda_L$-shock and the $\gamma_L$-contact discontinuity. More precisely, we look for an improvement of Lemma 7 in the case of coalescence, i.e. $h_{L*} \ge \eta_L$ to define the nonconservative product. As in Section 3.2, we use the regularization of the Heaviside function $H_\epsilon(z)$ satisfying (28) and we look for a regularized profiles of the velocities $(\bar{u}_\epsilon, \hat{u}_\epsilon, \bar{v}_\epsilon, \hat{v}_\epsilon)(z)$ solution of (9) with $(h_\epsilon, \bar{u}_\epsilon, \hat{u}_\epsilon, \bar{v}_\epsilon, \hat{v}_\epsilon)(0) = (h_L, \bar{u}_L, \hat{u}_L, \bar{v}_L, \hat{v}_L)$.

**Lemma 10.** Let the left state $(h_L, \bar{u}_L, \hat{u}_L)$ and the right state $(h_{L*}, \bar{u}_*, \hat{u}_{L*})$ be two states linked by a $\lambda_L$-shock with coalescence, i.e. $h_{L*} > \eta_L$ with $\sigma_L$ the speed of this shock. Let $h_\epsilon$ defined by (29) and $[z_0, z_1]$ the interval such on
which $h_{\varepsilon}(z \in [z_0, z_1]) = \eta_L$. Consider $V_{\varepsilon} = (h_{\varepsilon}, \bar{u}_{\varepsilon}, \hat{u}_{\varepsilon}, \bar{v}_{\varepsilon}, \hat{v}_{\varepsilon})$ a solution of (9) piecewise $C^1(\mathbb{R} - [z_0, z_1])$, depending only on $z = x - \sigma_L t$, such that $V_{\varepsilon}(0) = (h_L, \bar{u}_L, \hat{u}_L, \bar{v}_L, \hat{v}_L)$. Then $(h_{\varepsilon}, \bar{u}_{\varepsilon}, \hat{u}_{\varepsilon})(\varepsilon) = (h_{L*}, \bar{u}_{*}, \hat{u}_{L*})$ and $(\bar{v}_{\varepsilon}, \hat{v}_{\varepsilon})(\varepsilon)$ belong to the following set $V_L^c$, independent of $\varepsilon$:

$$V_L^c = \left\{ (\bar{v}_{L*}, \hat{v}_{L*}) \in \mathbb{R}^2 \text{ such that } \begin{array}{l} \bullet h_{L*}((\bar{u}_{L*} - \sigma_L)\bar{v}_{L*} + \hat{u}_{L*}\hat{v}_{L*}) = h_L((\bar{u}_L - \sigma_L)\bar{v}_L + \hat{u}_L\hat{v}_L) \end{array} \right\}.$$

Conversely, if $(\bar{v}_{L*}, \hat{v}_{L*})$ belong to $V_L^c$, there exists a unique piecewise $C^1(\mathbb{R} - [z_0, z_1])$ solution of (9), depending only on $z = x - \sigma t$, that takes value $(h_{L*}, \bar{u}_{L*}, \hat{u}_{L*}, \bar{v}_{L*}, \hat{v}_{L*})$ at $z = \varepsilon$.

**Proof.** Similarly to the proof of Lemma 7, we obtain the ODEs given by System (31) out of the intervals $z \in [z_0, z_1]$. It follows that there exists a constant $\hat{v}_{L*}$ such that the solution reads

$$\hat{v}_{\varepsilon}(z) = \begin{cases} \frac{h_L}{h_{\varepsilon}(z)} \sqrt{\frac{\bar{w}_L^2 - \hat{w}_L^2}{\bar{w}_{\varepsilon}(z)^2 - \hat{w}_{\varepsilon}(z)^2}} |\hat{v}_L| & \text{for } z \in [0, z_0], \\ \frac{h_L}{h_{\varepsilon}(z)} \sqrt{\frac{\bar{w}_L^2 - \hat{w}_L^2}{\bar{w}_{\varepsilon}(z)^2 - \hat{w}_{\varepsilon}(z)^2}} |\hat{v}_{L*}| & \text{for } z \in [z_1, \varepsilon], \end{cases}$$

which concludes the proof. □

Thanks to Lemma 10 we extend the definition of the nonconservative product Definition 1 in the case of a shock with coalescence.

**Definition 2.** Through a $\lambda_k$-shock with coalescence ($k \in \{L, R\}$), i.e. $h_{k*} \ge \eta_k$, the following jump conditions hold:

$$\left\{ \begin{array}{l} \sigma_k[h] = [h\bar{u}] \\ \sigma_k[h\bar{u}] = [h(\bar{u}^2 + \hat{u}^2) + \frac{g}{2}h^2] \\ \sigma_k[\hat{u}] = [\bar{u}\hat{u}] \\ \sigma_k[h\bar{v}] = [h(\bar{u}\bar{v} + \hat{u}\hat{v})]. \end{array} \right.$$

with $\sigma_k$ is the speed of the shock.

We now focus on the mechanical energy dissipation through a shock with coalescence.

**Proposition 7.** We denote by $\sigma_k$ the speed of the $\lambda_k$-shock with coalescence, i.e. $h_{L*} \ge \eta_L$ or $h_{R*} \ge \eta_R$. Then the lax condition is satisfied Proposition 4.iii and the mechanical energy is decreasing, i.e. Proposition 6.i.

**Proof.** We focus on the $\lambda_L$-shock with coalescence. Similar arguments hods for the $\lambda_R$-shock with coalescence. Since $h_{L*} \ge \eta_L$, we have $Q[h] > 0$ and thank to Proposition 4, we conclude that Proposition 4.iii holds.
On the other hands, similarly to the proof of Proposition 6 the jump of the flux of mechanical energy in the framework of the shock is given by (32). However, this term is not vanishing since the last jump relation of Definition 1 does not occur in the case of coalescence. Using (21) and (32), the jump of energy in the framework of the shock reads

$$-\sigma_L [E] + \left[ \left( \frac{\bar{u}^2 + 3\hat{u}^2 + \bar{v}^2 + \hat{v}^2}{2} + gh \right) h\bar{u} + h\bar{v}\hat{u}\hat{v} \right] = -Q[h] \left( \frac{g + 2S^2(h_L + h_{L*})}{4h_Lh_{L*}} [h]^2 - \frac{[(Q^2 - S^2h^4)\hat{v}^2]}{2Q^2[h]} \right).$$

Remark that in the case of coalescence, we have $\bar{w}_L^2 - \hat{w}_L^2 > 0$ while $\bar{w}_{L*}^2 - \hat{w}_{L*}^2 < 0$. It follows that the jump $\frac{[(Q^2 - S^2h^4)\hat{v}^2]}{[h]} < 0$ and we conclude. □

With Definition 2, we are in position to complete the description of the $\sigma_L$-shock wave in the case of a coalescence $\gamma_L$-contact discontinuity. In this case the state $V_{L1}$ disappears and the left state $V_L$ is directly linked with the state $V_{L2}$, see Figure 4.

**Lemma 11.** Assuming the jump conditions of Definition 2, the admissible right state $V_{L2} = (h_{L*}, \bar{u}_*, \hat{u}_{L*}, \bar{v}_{L2}, \hat{v}_{L2})^t$ that can be reached from the left state $V_L = (h_L, \bar{u}_L, \hat{u}_L, \bar{v}_L, \hat{v}_L)^t \in \mathbb{R}_+^* \times \mathbb{R}^4$ with a $\lambda_L$-shock which coalesce with a $\gamma_L$-contact discontinuity, i.e. $h_{L*} \ge \eta_L$, satisfies the condition of Lemma 2 and

$$\bar{v}_{L2} = \bar{v}_L + \left( \frac{\hat{u}_L}{\bar{u}_L - \sigma_L} \hat{v}_L - \frac{\hat{u}_{L*}}{\bar{u}_* - \sigma_L} \hat{v}_{L2} \right).$$

*Proof.* The result directly yields from the fourth jump condition of Definition 2 using the first and the third jump conditions. □

The previous results can as well be established for the right wave.

**Lemma 12.** Assuming the jump conditions of Definition 2, the admissible left state $V_{R2} = (h_{R*}, \bar{u}_*, \hat{u}_{R*}, \bar{v}_{R2}, \hat{v}_{R2})^t$ that can be reached from the right state $V_R = (h_R, \bar{u}_R, \hat{u}_R, \bar{v}_R, \hat{v}_R)^t \in \mathbb{R}_+^* \times \mathbb{R}^4$ with a $\lambda_L$-shock which coalesce with a $\gamma_L$-contact discontinuity, i.e. $h_{R*} \ge \eta_R$, satisfies the condition of Lemma 2 and

$$\bar{v}_{R2} = \bar{v}_R - \left( \frac{\hat{u}_R}{\bar{u}_R - \sigma_R} \hat{v}_R - \frac{\hat{u}_{R*}}{\bar{u}_* - \sigma_R} \hat{v}_{R2} \right).$$

*Proof.* The result directly yields from the fourth jump condition of Definition 2 using the first and the third jump conditions. □

We are now in position to solve the Riemann problem (9) with (10) in the case with coalescence.
**Theorem 3.** The result of Theorem 2 holds even with coalescence, i.e. $h_{L*} \ge \eta_L$ and/or $h_{R*} \ge \eta_R$.

*Proof.* For the sake of completeness, for $k \in \{L, R\}$, we give a compact expression of $\bar{v}_{k1}$, gathering the results obtain in Section 3.2 with the results obtain in Section 3.3.

$$ \bar{v}_{k1} = \begin{cases} \bar{v}_k - \frac{\hat{u}_k (\bar{u}_k - \bar{u}_*)}{\frac{g}{2} h_k + \hat{u}_k^2} \hat{v}_k & \text{if } h_{k*} \le h_k \quad \text{(rarefaction)}, \\ \bar{v}_k - \left( \frac{h_{k*}}{h_k} \sqrt{\frac{(\bar{u}_k - \sigma_k)^2 - \hat{u}_k^2}{(\bar{u}_* - \sigma_k)^2 - \hat{u}_{k*}^2}} - 1 \right) \frac{\hat{u}_k}{\bar{u}_k - \sigma_k} \hat{v}_k & \text{if } h_k < h_{k*} < \eta_k \quad \text{(shock)}, \\ \bar{v}_k & \text{if } \eta_k \le h_{k*} \quad \text{(coalescence)}. \end{cases} $$

and similarly for the standard deviation velocity $\hat{v}_{k1}$ reads

$$ (35) \qquad \hat{v}_{k1} = \begin{cases} \frac{g}{2} h_k^2 + h_{k*} \hat{u}_k^2 \hat{v}_k & \text{if } h_{k*} \le h_k \quad \text{(rarefaction)}, \\ \frac{g}{2} h_k^2 + h_k \hat{u}_k^2 & \text{if } h_k < h_{k*} < \eta_k \quad \text{(shock)}, \\ \frac{h_k}{h_{k*}} \sqrt{\frac{(\bar{u}_k - \sigma_k)^2 - \hat{u}_k^2}{(\bar{u}_* - \sigma_k)^2 - \hat{u}_{k*}^2}} \hat{v}_k & \text{if } \eta_k \le h_{k*} \quad \text{(coalescence)}. \end{cases} $$

In the case of coalescence the definition of $\bar{v}_{k1}$ is extended for homogeneity but do not corresponds to the solution of the Riemann problem in a part of the domain. Eventually, following the proof of Theorem 2, we obtain (33) with the parameters $a_k$ and $b_k$ given by

$$ (36) \qquad \begin{aligned} a_k &= \begin{cases} 1 & \text{if } h_{k*} \le \eta_k \quad \text{(rarefaction or shock)}, \\ \frac{|\hat{u}_{k*}|}{\bar{u}_* - \sigma_k} & \text{if } \eta_k \le h_{k*} \quad \text{(coalescence)}, \end{cases} \\ \text{and } b_k &= \begin{cases} 1 & \text{if } h_{k*} \le \eta_k \quad \text{(rarefaction or shock)}, \\ \frac{|\hat{u}_k|}{\bar{u}_k - \sigma_k} & \text{if } \eta_k \le h_{k*} \quad \text{(coalescence)}. \end{cases} \end{aligned} $$

The two last $J_{\lambda_*}$-Riemann invariants leads to the linear system (34) and we conclude similarly to the proof of Theorem 2. $\square$

## 3.4 Resonance of the contact discontinuities

The last point to study is the case of resonance, i.e. $\hat{u} = 0$. If $\hat{u} = 0$ (and $\hat{v} \ne 0$), the system is not hyperbolic as a Jordan block appears in the Jacobian matrix $A(V)$, see (27). First of all, when $\hat{u}_L = \hat{u}_R = 0$, the three internal contacts are in resonance. However, the two first Riemann invariants of $J_{\gamma_L}$, $J_{\lambda_*}$ and $J_{\gamma_R}$ are the same namely $\bar{u}$ and $h$, while $\hat{u}$ is null everywhere. It follows from (9) that $\bar{v}$ and $\hat{v}$ are simply advected with the velocity $\bar{u}$ thus $(\bar{v}_{L1}, \hat{v}_{L1}) = (\bar{v}_L, \hat{v}_L)$ and $(\bar{v}_{R1}, \hat{v}_{R1}) = (\bar{v}_R, \hat{v}_R)$ and the states $V_{L2}$ and $V_{R2}$ disappear. This results is consistent with formulae (35) even if the linear system (34) is singular in this
case. This case can be compare to the classical shallow water model with two passive pollutants namely $\bar{v}$ and $\hat{v}$.

Now assume that $\hat{u}_L \neq 0$ and $\hat{u}_R = 0$ (the other case can be treated in a similar way) then the $\gamma_R$ and the $\lambda_*$-contact discontinuities are in resonance and the state $V_{R2}$ disappears. In order to identify the jump condition through the resulting contact discontinuity, it is more convenient to use the energy conservation (Proposition 2) instead of the nonconservative equation on $\hat{v}$ (last equation of (9)). This formulation is equivalent since we focus on the simple contact discontinuity wave. The Rankine-Hugoniot jump condition reads

$$ (37) \qquad \begin{aligned} [\bar{u}] &= 0, \\ [h\hat{u}^2 + \frac{g}{2}h^2] &= 0, \\ [h\hat{u}\hat{v}] &= 0, \\ [h\bar{v}\hat{u}\hat{v}] &= 0. \end{aligned} $$

Using Lemma 4 we deduce that $\hat{u}_{R*} = 0$ then from the two last jump conditions of (37) it yields to $\hat{v}_{L2} = 0$. Since the $\gamma_L$-contact discontinuities is not in resonance with the $\lambda_*$-contact discontinuities (but can coalesce with the $\lambda_L$-shock), the $\bar{v}_{L2}$ is given by (33) and yields to $\bar{v}_{L2} = \bar{v}_{L1} + s_L b_L \hat{v}_{L1}$ with $b_L$ given by (36). Not that this result is consistent with the solution of (33) with (34) when $\hat{u}_{R*} = 0$ and $\hat{u}_{L*} \neq 0$. More precisely, even if System (34) is singular (we can not determine $\hat{v}_{R2}$), we recover that $\hat{v}_{L2}$ vanishes then using (33) we obtain the previous result. Moreover, we are able to estimate the limit of the state $V_{R2}$ when $\hat{u}_{R*}$ goes to zero (different from zero). More precisely, when $\hat{u}_R$ goes to zero $\eta_R$ goes to infinity thus the $\lambda_R$ and $\gamma_R$-waves does not coalesce, then $a_R = b_R = 1$. Using (34), we get

$$ \begin{aligned} & \hat{v}_{R2} \underset{\hat{u}_R \to 0}{\longrightarrow} s_R \bar{v}_{L1} + s_L s_R b_L \hat{v}_{L1} - s_R \bar{v}_{R1} + \hat{v}_{R1} \\ & \bar{v}_{R2} \underset{\hat{u}_R \to 0}{\longrightarrow} \bar{v}_{L1} + s_L b_L \hat{v}_{L1} \end{aligned} $$

with $(\bar{v}_{L1}, \hat{v}_{L1})$ given by Lemma 6 or Lemma 8 and $(\bar{v}_{R1}, \hat{v}_{R1})$ given by Lemma 9. We conclude that the velocities $(\bar{v}_{R2}, \hat{v}_{R2})$ remain bounded. Remark that the limit of $\hat{v}_{R2}$ is not unique, i.e. it is not the same if $\hat{u}_R \to 0_+$ or $\hat{u}_R \to 0_-$, which can be see as a discontinuity of the solution with respect to the initial condition. In addition the velocities $\bar{v}_{L2}$ and $\bar{v}_{R2}$ tends to the same values, only depending on the left side. Eventually, we conclude the following results.

**Theorem 4.** The result of Theorem 3 holds even with resonance, i.e. $\hat{u}_L = 0$ and/or $\hat{u}_R = 0$.

# 4 Some analytical solutions

In the Table 1, the approximated values of the solution of the Riemann problem for several initial conditions are given.
<table><caption>Table 1: Some analytical solutions of the Riemann problem</caption><thead><tr><th></th><th>V<sub>L</sub></th><th>V<sub>L1</sub></th><th>V<sub>L2</sub></th><th>V<sub>R2</sub></th><th>V<sub>R1</sub></th><th>V<sub>R</sub></th><th>wave pattern</th></tr></thead><tbody><tr><td>case 1</td><td>(2<br>0<br>2<br>2<br>1)</td><td>(1.58<br>1.24<br>1.58<br>2.18<br>0.94)</td><td>(1.58<br>1.24<br>1.58<br>3.92<br>-0.80)</td><td>(1.27<br>1.24<br>-2.54<br>3.92<br>0.62)</td><td>(1.27<br>1.24<br>-2.54<br>1.16<br>3.38)</td><td>(1<br>0<br>-2<br>2<br>3)</td><td></td></tr><tr><td>case 2</td><td>(1<br>0<br>2<br>2<br>1)</td><td>(1<br>0<br>2<br>2<br>1)</td><td>(1<br>0<br>2<br>4<br>-1)</td><td>(1<br>0<br>-2<br>4<br>1)</td><td>(1<br>0<br>-2<br>2<br>3)</td><td>(1<br>0<br>-2<br>2<br>3)</td><td></td></tr><tr><td>case 3</td><td>(2<br>0<br>0<br>2<br>1)</td><td>(1.45<br>1.31<br>0<br>2<br>1)</td><td></td><td></td><td>(1.45<br>1.31<br>0<br>2<br>3)</td><td>(1<br>0<br>0<br>2<br>3)</td><td></td></tr><tr><td>case 4</td><td>(2<br>0<br>0<br>2<br>1)</td><td>(1.54<br>0.54<br>0<br>2<br>1)</td><td></td><td>(1.12<br>0.54<br>-2.23<br>4.80<br>0)</td><td>(1.12<br>0.54<br>-2.23<br>1.64<br>3.16)</td><td>(1<br>0<br>-2<br>2<br>3)</td><td></td></tr><tr><td>case 5</td><td>(2<br>6<br>2<br>2<br>1)</td><td>(3.56<br>1.82<br>3.56<br>1.34<br>1.31)</td><td>(3.56<br>1.82<br>3.56<br>2.56<br>0.86)</td><td>(2.64<br>1.82<br>-5.27<br>2.56<br>-0.78)</td><td></td><td>(1<br>-6<br>-2<br>2<br>3)</td><td></td></tr><tr><td>case 6</td><td>(2<br>6<br>0<br>2<br>1)</td><td>(4.47<br>1.35<br>0<br>2.50<br>0.50)</td><td></td><td>(2.55<br>1.35<br>-5.09<br>2.50<br>0)</td><td></td><td>(1<br>-6<br>-2<br>2<br>3)</td><td></td></tr><tr><td>case 7</td><td>(2<br>13<br>2<br>2<br>1)</td><td></td><td>(5.69<br>2.87<br>5.69<br>2.19<br>-0.61)</td><td>(4.04<br>2.87<br>-8.08<br>2.19<br>0.60)</td><td></td><td>(1<br>-13<br>-2<br>2<br>3)</td><td></td></tr></tbody></table>
* The case 1 corresponds to a classical situation without coalescence or resonance. The five waves are distinct with a rarefaction wave at the left and a shock wave at the right. Remark that the standard deviation $\hat{u}$, which can be interpreted as the shear, can be larger (in absolute value) in the internal states. In this case we have $|\hat{u}_{R*}| > \max(|\hat{u}_L|, |\hat{u}_R|)$. More precisely, according to Proposition 4.ii, the standard deviation increases (in absolute value) through a shock.

* The case 2 is relevant in the sense that it is a steady state of the 1D Riemann problem (11). The velocity $v_1 = \bar{v} - \hat{v}$ is advected by the velocity $u_1 = \bar{u} - \hat{u}$ and the velocity $v_2 = \bar{v} + \hat{v}$ is advected by the velocity $u_2 = \bar{u} + \hat{u}$ except at the centered contact discontinuity where the tangential velocities switch from one layer to the other.

* In the cases 3 and 4, two or three internal waves are in resonance. The water depth $h$ and the mean velocity $\bar{u}$ of the case 3 are exactly those of the classical shallow water model and the velocities $(\bar{v}, \hat{v})$ are simply advected with the velocity $\bar{u}$. In the case 4 by $\hat{u}_L = 0$, the velocity $\hat{v}_{R2}$ vanishes as it was expected.

* In the cases 5, 6 and 7, the variation of mean velocity $\bar{u}_L - \bar{u}_R$ is large enough to lead to coalescence. In the case 5, the velocity $\hat{u}_L$ is small enough and the coalescence occurs only at the right. Here the left wave is a shock, however it might be a rarefaction by reducing the variation of mean velocity $\bar{u}_L - \bar{u}_R$ and by increasing the velocity $\hat{u}_R$, while maintaining the coalescence of the right waves. In the case 6, the velocity $\hat{u}_L$ vanishes such that the $\gamma_L$ and the $\lambda_*$-contact discontinuity are in resonance while the right waves still coalesce. Eventually in the case 7, choosing $\hat{u}_L$ large enough, coalescences occur on both sides of the $\lambda_*$-wave.

Conclusion

In the current work, we have established the global well-posedness of the Riemann problem for a two dimensional enriched shallow water model, that involves not only the mean velocity but also a correction that can be interpreted either as an estimation of the standard deviation to the vertical averaged value or through the description of the flow with two layers of different velocities. This result makes our work a first step in the analysis of the hyperbolicity of the layerwise models proposed in [4]. A new formulation of the system derived in [4] is studied, it is conservative for all the equations except the standard deviation on the transverse velocity. Then a definition of the nonconservative product which is path independent is proposed by using the tools developed in [1, 15]. The well posedness is first proved assuming the waves are separated (Theorem 2). Indeed, for large initial data, a coalescence phenomenon can appear at the extreme waves and is analyzed, leading to Theorem 3. Resonance can also occur in the internal waves, however, thanks to the linear degeneracy of these waves
it does yield nonuniqueness and Theorem 4 completes the previous results.

As already mentioned, the bilayer case is not satisfying if the objective is an approximation of the Euler equations. A natural extension of the current work would be the analysis for a larger, ideally arbitrary, number of layers. This goal seems out of reach in general regime but the analysis could be performed in some asymptotic regime. In particular, the weak shear regime, i.e. $u_{i+1} - u_i \ll \sqrt{gh}$ is relevant for application in oceanography. Let us mention that an analysis of another type of layerwise model, without mass exchange, was performed in such regime in [18, 23, 9, 19]. Another possible improvement of the current work will be to consider the layerwise discretization of a non-homogeneous density flow [3]. The analysis of the Riemann problem, or even the hyperbolicity of the model, is not trivial even assuming a particular layer thickness and velocity at the interface.

## References

[1] AGUILLON, N. Riemann problem for a particle-fluid coupling. *Mathematical Models and Methods in Applied Sciences* 25, 01 (2017/03/27 2014), 39–78.

[2] AUDUSSE, E. A multilayer Saint-Venant model: derivation and numerical validation. *Discrete Contin. Dyn. Syst. Ser. B* 5 (2005), 189–214.

[3] AUDUSSE, E., BRISTEAU, M. O., PELANTI, M., AND SAINTE-MARIE, J. Approximation of the hydrostatic Navier-Stokes system for density stratified flows by a multilayer model: Kinetic interpretation and numerical solution. *Journal of Computational Physics* 230, 9 (5 2011), 3453–3478.

[4] AUDUSSE, E., BRISTEAU, M.-O., PERTHAME, B., AND SAINTE-MARIE, J. A multilayer Saint-Venant system with mass exchanges for shallow water flows. Derivation and numerical validation. *ESAIM: M2AN* 45, 1 (2011), 169–200.

[5] BRISTEAU, M.-O., GUICHARD, C., DI MARTINO, B., AND SAINTE-MARIE, J. Layer-averaged Euler and Navier-Stokes equations. working paper or preprint, hal-01202042.

[6] CASTRO, A., AND LANNES, D. Fully nonlinear long-wave models in the presence of vorticity. *Journal of Fluid Mechanics* 759 (2014), 642–675.

[7] CHINAYA, A., LEROUX, A., AND SEGUIN, N. A well-balanced numerical scheme for the approximation of the shallow-water equations with topography: the resonance phenomenon. *International Journal on Finite Volume* 1 (2004), 1–33.
[8] DAL MASO, G., LEFLOCH, P. G., AND MURAT, F. Definition and weak stability of nonconservative products. *Journal de mathématiques pures et appliquées* **74**, 6 (1995), 483–548.

[9] DUCHÊNE, V. Asymptotic shallow water models for internal waves in a two-fluid system with a free surface. *SIAM J. Math. Anal.* **42** (2010), 2229–2260.

[10] FERNANDEZ-NIETO, E. D., PARISOT, M., PENEL, Y., AND SAINTE-MARIE, J. Layer-averaged approximation of Euler equations for free surface flows with a non-hydrostatic pressure. working paper or preprint, hal-01324012.

[11] GAVRILYUK, S., IVANOVA, K., AND FAVRIE, N. Multi-dimensional shear shallow water flows : problems and solutions. working paper or preprint, hal-01529497.

[12] GOATIN, P., AND LEFLOCH, P. G. The riemann problem for a class of resonant hyperbolic systems of balance laws. *Annales de l'Institut Henri Poincare (C) Non Linear Analysis* **21**, 6 (2004), 881 – 902.

[13] GODLEWSKI, E., AND RAVIART, P.-A. *Numerical approximation of hyperbolic systems of conservation laws*, vol. 118 of *Applied Mathematical Sciences*. Springer-Verlag, New York, 1996.

[14] ISAACSON, E., AND TEMPLE, B. Convergence of the 2x2 Godunov Method for a General Resonant Nonlinear Balance Law. *SIAM Journal on Applied Mathematics* **55**, 3 (1995), 625–640.

[15] LAGOUTIÈRE, F., SEGUIN, N., AND TAKAHASHI, T. A simple 1d model of inviscid fluid-solid interaction. *Journal of Differential Equations* **245**, 11 (2008), 3503 – 3544.

[16] LEVEQUE, R. J. *Finite volume methods for hyperbolic problems*, vol. 31. Cambridge university press, 2002.

[17] LISKA, R., AND WENDROFF, B. Analysis and computation with stratified fluid models. *J. Comput. Phys.* **137**, 1 (1997), 212–244.

[18] LONG, R. R. Long waves in a two-fluid system. *Journal of Meteorology* **13**, 1 (1956), 70–74.

[19] MONJARRET, R. Local well-posedness of the two-layer shallow water model with free surface. *SIAM Journal on Applied Mathematics* **75**, 5 (2016/04/11 2015), 2311–2332.

[20] OVSYANNIKOV, L. Two-layer “shallow water” model. *J. Appl. Mech. Tech. Phys.* **20**, 2 (1979), 127–135.
[21] RAMBAUD, A. *Modélisation, analyse mathématique et simulations numériques de quelques problèmes aux dérivées partielles multi-échelles*. PhD thesis, 2011.

[22] RICHARD, G. L., AND GAVRILYUK, S. L. A new model of roll waves: comparison with brock's experiments. *Journal of Fluid Mechanics* **698** (005 2012), 374–405.

[23] SCHIJF, J., AND SCHÖNFLED, J. Theoretical considerations on the motion of salt and fresh water. IAHR, 1953.

[24] TESHUKOV, V. M. Gas-dynamic analogy for vortex free-boundary flows. *Journal of Applied Mechanics and Technical Physics* **48**, 3 (2007), 303–309.

[25] TORO, E. F. *Riemann solvers and numerical methods for fluid dynamics: a practical introduction*. Springer Science & Business Media, 2013.