The Riemann problem for two-layer shallow water equations with bottom topography

Fan Wang; Wancheng Sheng; Yanbo Hu; Qinglong Zhang

doi:10.1515/anona-2024-0058

Article Open Access

The Riemann problem for two-layer shallow water equations with bottom topography

Fan Wang , Wancheng Sheng , Yanbo Hu and Qinglong Zhang

Published/Copyright: December 7, 2024

Published by

Become an author with De Gruyter Brill

Submit Manuscript Author Information Explore this Subject

From the journal Advances in Nonlinear Analysis Volume 13 Issue 1

Abstract

In this article, we propose a model describing laminar shallow water flows with two velocities to handle an uneven bottom. The model is established by taking the different velocities into account. The source terms generated from the bottom topography prevent us from solving the Riemann problem directly. We first derive the elementary waves including the stationary wave where the global entropy condition is used to ensure uniqueness. Then, we analyze the resonance phenomenon and coalescence of waves by classifying the initial data into different regions. Finally, the Riemann problem is resolved explicitly on a case-by-case basis.

Keywords: the Riemann problem; laminar shallow water equations; Riemann invariants; the stationary wave; Rankine-Hugoniot condition

MSC 2010: Primary: 35L40; 76L05; 76M12; 35L65; Secondary: 35F55

1 Introduction

Recently, the multilayer shallow water equations have been of great interest. On the one hand, it is very efficient in describing flows with different densities and varying velocities both theoretically and numerically. On the other hand, the multilayer models have wide applications in the two-phase and multiphase flows, which play a key role in the research of modern industry. The two-layer shallow water system is one of the representatives of the multilayer flows, it can be obtained by taking the vertical averaging of the layer depth [1,2,4]. In this work, we consider the two-layer shallow water equations with non-flat bottom topography as follows:

(1.1) ( h 1 ) t + ( h 1 u 1 ) x = m e , ( h 1 u 1 ) t + h 1 u 1 2 + g h 1 2 2 + g h 1 h 2 x = g h 2 ( h 1 ) x + u in m e − g h 1 B ′ ( x ) , ( h 2 ) t + ( h 2 u 2 ) x = − m e , ( h 2 u 2 ) t + h 2 u 2 2 + g h 2 2 2 x = − g h 2 ( h 1 ) x − u in m e − g h 2 B ′ ( x ) ,

where h i and u i ( i = 1 , 2 ) represent the approximation of the layer thickness and the horizontal velocity of the i th layer, m e denotes the mass exchange from the second layer to the first layer, u in is the velocity of the interface between two layers, B ( x ) denotes the bottom elevation. The system contains the source terms resulting from the bottom topography and mass exchange between the layers. For more details about the derivation of the model, we refer to [2,6] and references therein.

The reasons we are interested in the two-layer shallow water equations (1.1) are twofold: (i) It is equivalent to the bilayer version of the layerwise model introduced in [6], which was used to approximate the 3D hydrostatic equations; (ii) System (1.1) also takes into consideration of the mass exchange between the layers, which can be seen as an extension of the model proposed in [4]. There are a lot of research studies on the two-layer flows in the literature [7,11,20,23–25]. Audusse et al. [6] introduced a multilayer model with mass exchange between the neighboring layer from the Navier-Stokes system. The local well-posedness of the two-layer shallow water equations was obtained in [22]. Abgrall and Karni [1] introduced a relaxation approach where two artificial equations have been added so that the extended system is hyperbolic as the original two-layer flow is conditionally hyperbolic (see also in [14]). Numerical simulations are also available. A well-balanced scheme was introduced in [16] for both one-dimensional (1D) and two-dimensional (2D) two-layer shallow water equations. A time-splitting approach proposed in [15] is implemented without having the full eigenstructure of the system. Several upwind schemes were developed in the last decade, including the finite-volume and finite-element types [1,14]. More related works can be found in [5,9,10, 12,13,26].

Recently, an interesting approach dealing with the two-layer model with flat bottom has been proposed in [2]. For convenience, we post the 1D case here for later use,

(1.2) h t + ( h u ¯ ) x = 0 , ( h u ¯ ) t + h ( u ¯ 2 + u ˆ 2 ) + g 2 h 2 x = 0 , u ˆ t + ( u ˆ u ¯ ) x = 0 ,

where h denotes the total height of the water flows, the notations u ¯ and u ˆ represent the average of the vertical velocity and the oriented standard deviation, respectively, which are given by

(1.3) u ¯ = u 1 + u 2 2 , u ˆ = u 2 − u 1 2 .

In order to derive (1.2), the authors have made the assumption that all the layers have a homogeneous thickness, i.e., h 1 = h 2 , and the interface velocity is defined by u in = u ¯ . The detailed derivations can be found in [2]. One may note that system (1.2) is hyperbolic and can be used to describe laminar flows with two velocities. The major advantage of investigating (1.2) is that the Riemann problem has been fully solved for both 1D and 2D cases, which provides direct Riemann solvers for numerical simulations of the two-layer flows. Another reason we are focusing on (1.2) is that the wave interaction and numerical investigation on (1.2) are still lacking, which have already been studied in [8,17,19,30] for the single-layer shallow water equations. The comparison of system (1.2) with other two-layer models are also given in [2].

In this work, we are interested in extending the results obtained in [2] to the non-flat bottom case (1.1), which has a more practical application when considering the multilayer flows. Besides, it is also interesting to investigate the Riemann structures with the existence of the bottom function B ( x ) . We should emphasize here that although the extension is quite natural, the bottom topography introduces a non-conservative term in system (1.2) as we will see below, which is the major difficulty in solving the Riemann problem due to the non-strictly eigenstructure of the system. As indicated in [21], the weak solution can be obtained by defining a path connecting the left-hand and the right-hand states W L and W R , respectively. For the two-layer shallow water equations, it does not offer a clear definition of such a path. A different approach is proposed in [18] for duct flows by adding a trivial equation where a linear degenerate characteristic field was introduced. This approach is efficient and shows its robustness in the two-phase flows [3], the blood flows in arteries [27,28], traffic flow models [29] and so on. In this article, we take this approach and introduce a linear degenerate characteristics for the bottom function B ( x ) . We then analyze the characteristics of the system before solving the Riemann problem. The elementary waves include rarefaction wave, shock wave, and particularly, stationary wave. After we present the main property of the stationary wave, the Riemann problem is solved explicitly on a case-by-case basis. The resonance phenomenon and coalescence of waves are also considered. The solution may be useful for the design of numerical schemes in future study.

This article is organized as follows. Section 2 includes some preliminaries and notations for (1.1). Particularly, we give the detailed derivation of the two-layer shallow water equations with bottom topography. In Section 3, the elementary waves for the presented system are concerned, which includes the rarefaction wave, the shock wave, and the stationary wave. In particular, the global entropy condition is proposed to ensure the uniqueness of the stationary solution. Section 4 deals with the Riemann problem explicitly. The initial datum is classified in the phase plane and the Riemann solution is constructed case by case where the resonance phenomenon and coalescence of waves are concerned.

2 Preliminaries

2.1 Model derivation and characteristics analysis

In this section, we begin to derive the two-layer shallow water equations with bottom topography. Basically, we follow the assumption made in [2] that the two layers share the same thickness h 1 = h 2 = h ⁄ 2 . By adding the two mass equations of (1.1), it yields

(2.1) h t + h 2 ( u 1 + u 2 ) x = 0 .

Similarly, the summation of the momentum equations in (1.1) leads to

(2.2) h 2 ( u 1 + u 2 ) t + h 2 ( u 1 2 + u 2 2 ) + g h 2 2 x = − g h B ′ ( x ) .

Besides, in order to make the system completed, we subtract the two momentum equations, which yields

(2.3) ( u 2 − u 1 ) t + u 2 − u 1 2 ( u 1 + u 2 ) x = 0 ,

where we have used that

(2.4) m e = − ( h ( u 2 − u 1 ) ) x 2 ,

which is obtained by just subtracting the mass equations.

By using the notations u ¯ , u ˆ defined in (1.3), and combine (2.1), (2.2), and (2.3), we arrive at

(2.5) h t + ( h u ¯ ) x = 0 , ( h u ¯ ) t + h ( u ¯ 2 + u ˆ 2 ) + g 2 h 2 x = − g h B ′ ( x ) , u ˆ t + ( u ˆ u ¯ ) x = 0 .

We note that system (2.5) has a time-independent function B ( x ) . In order to investigate the Riemann problem of (2.5), it is convenient to introduce a trivial equation associated with B ( x ) , namely

(2.6) B t = 0 .

This equation will be our starting point to define the weak solution of the Riemann problem, which is also used in the duct flows and traffic flows [18,29].

Before we make characteristic analysis of system (2.5), we can compare (1.2) with the Euler equations. Recall in [2], the authors have made the analogous with (1.2) and the Euler equations by setting

(2.7) p = g 2 h 2 + h u ˆ 2 , S = u ˆ ⁄ h ,

where p and S are the pressure and the entropy in the gas dynamics, respectively. However, the main difference between system (2.5) and the Euler system is h cannot be expressed explicitly as a function of p , which brings difficulties in the theoretical and numerical analysis. For notation convenience, here we use the variables ( h , u ¯ , S ) instead of the original variables ( h , u ¯ , u ˆ ) to express the system. Of course, here p and S do not have the meanings of pressure and entropy as in gas dynamics. As is pointed out in [2], system (2.5) is analogous to the full Euler equations at least for smooth solutions under the transformation (2.7). With the help of (2.7), one can rewrite (2.5) and (2.6) as follows:

(2.8) U t + A ( U ) U x = 0 ,

with

(2.9) U = h u ¯ S B , A ( U ) = u ¯ h 0 0 g + 3 S 2 h u ¯ 2 S h 2 g 0 0 u ¯ 0 0 0 0 0 .

The matrix A ( U ) has four eigenvalues, namely

(2.10) λ 1 = u ¯ − c , λ 2 = u ¯ , λ 3 = u ¯ + c , λ 4 = 0 ,

where c is the celerity given by

(2.11) c = g h + 3 S 2 h 2 .

We can further compute the corresponding eigenvectors as follows:

(2.12) r → 1 = ( − c , S h , h , 0 ) T ; r → 2 = ( 1 , 0 , 0 , 0 ) T ; r → 3 = ( c , S h , h , 0 ) T ; r → 4 = − u ¯ , S h , h , u ¯ 2 − c 2 g T .

One can see that λ 1 , λ 2 , and λ 3 may coincide with λ 4 , which means that the system is not strictly hyperbolic. More precisely, if we set

(2.13) Γ ± : u ¯ = ± c , Γ 0 : u ¯ = 0 ,

then

(2.14) λ 1 = λ 4 on Γ + , λ 3 = λ 4 on Γ − .

In the ( u ¯ , p ) phase plane, the three curves Γ ± and Γ 0 separate the whole plane by four parts such that in each part, system (2.5) is strictly hyperbolic. We should point out here that the construction of the Riemann problem is also based on the separation. For convenience, we will view these regions as supercritical regions ( D 1 , D 4 ), and subcritical regions ( D 2 , D 3 ), namely,

(2.15) D 1 = { ( u ¯ , p ) ∣ u > c } , D 2 = { ( u ¯ , p ) ∣ 0 < u < c } D 3 = { ( u ¯ , p ) ∣ − c < u < 0 } , D 4 = { ( u ¯ , p ) ∣ u < − c } .

3 Elementary waves

In this section, we present the elementary waves of systems (2.5) and (2.6). Generally, the elementary waves include the rarefaction waves, the shock waves, the contact discontinuity, and the stationary wave. We will discuss them in the following section.

3.1 Rarefaction wave

The rarefaction solution is a continuous solution that depends on the self-similar variable ξ = x ⁄ t . The solutions are associated with the two nonlinear characteristic fields λ 1 and λ 3 . We compute the Riemann invariants ω i corresponding to λ i ( i = 1 , 3 ) from

(3.1) ∇ ω i ⋅ r → i = 0 , i = 1 , 3 .

Consequently, the Riemann invariants ω 1 and ω 3 are given by [2]

(3.2) ω 1 = u ¯ + g h + 3 S 2 h 2 + g 3 S log 1 + 3 g S 2 h + 3 g S 2 h , B , S , ω 3 = u ¯ − g h + 3 S 2 h 2 − g 3 S log 1 + 3 g S 2 h + 3 g S 2 h , B , S .

From (3.2), we see that for the rarefaction wave, the bottom function B ( x ) remains constant; thus, the system degenerates to the flat case in [2]. For a given left-hand state ( h 0 , u ¯ 0 , S 0 ) , the 1-rarefaction wave curve that can be connected on the right is given by

(3.3) R ← 1 : u ¯ = u ¯ 0 − c + c 0 − g 3 S log 1 + 3 S 2 h g + 3 S 2 h g + g 3 S 0 log 1 + 3 S 0 2 h 0 g + 3 S 0 2 h 0 g , S = S 0 , h < h 0 .

Similarly, for a given left-hand state ( h 0 , u ¯ 0 , S 0 ) , the 3-rarefaction wave curve that can be connected on the right is given by

(3.4) R → 3 : u ¯ = u ¯ 0 + c − c 0 + g 3 S log 1 + 3 S 2 h g + 3 S 2 h g − g 3 S 0 log 1 + 3 S 0 2 h 0 g + 3 S 0 2 h 0 g , S = S 0 , h > h 0 .

3.2 Shock wave

For the discontinuous solution with the propagating speed σ , we start with the Rankine-Hugoniot condition associated with (2.6), which is given by

(3.5) σ [ B ] = 0 ,

where [ B ] = B + − B − is the jump of the bottom function B ( x ) . From (3.5), we have that either the shock speed vanishes ( σ = 0 ) or the bottom remains constant ( [ B ] = 0 ) across the non-zero shock.

We first look at the shock wave with non-zero speed, the bottom function keeps constant in this case. Thus, system (2.5) degenerates to

(3.6) h t + ( h u ¯ ) x = 0 , ( h u ¯ ) t + h ( u ¯ 2 + u ˆ 2 ) + g 2 h 2 x = 0 , u ˆ t + ( u ˆ u ¯ ) x = 0 .

The Rankine-Hugoniot condition associated with the above system reads

(3.7) σ [ h ] = [ h u ¯ ] , σ [ h u ¯ ] = h ( u ¯ 2 + u ˆ 2 ) + g 2 h 2 , σ [ u ˆ ] = [ u ˆ u ¯ ] .

In order to ensure uniqueness solution, we assume that the shock wave satisfies the Lax entropy condition

(3.8) λ i ( U ) < σ i ( U , U 0 ) < λ i ( U 0 ) , i = 1 , 3 ,

where U 0 and U represent the left- and right-hand states of the shock wave, respectively. By using the Lax entropy condition, we determine the 1-shock wave curve as follows:

(3.9) S ← 1 : u ¯ = u ¯ 0 − ( h − h 0 ) S 2 h h 0 + h 0 h + 1 + g 2 1 h 0 + 1 h , S = S 0 , h > h 0 .

Similarly, for a given left-hand state ( h 0 , u ¯ 0 , S 0 ) , the 3-shock wave curve that can be connected on the right is given by

(3.10) S → 3 : u ¯ = u ¯ 0 + ( h − h 0 ) S 2 h h 0 + h 0 h + 1 + g 2 1 h 0 + 1 h , S = S 0 , h < h 0 .

In order to solve the Riemann problem of (2.1), it is necessary to know when the shock speeds in the nonlinear characteristic fields equal to zero. From (3.7), for the shock speed σ , we have that

(3.11) σ = h u ¯ − h 0 u ¯ 0 h − h 0 = h ( u ¯ − u ¯ 0 ) + u ¯ 0 ( h − h 0 ) h − h 0 = u ¯ 0 − h S 2 h h 0 + h 0 h + 1 + g 2 1 h 0 + 1 h for 1-shock wave , u ¯ 0 + h S 2 h h 0 + h 0 h + 1 + g 2 1 h 0 + 1 h for 3-shock wave .

We take the 1-shock wave as an example. For a given left-hand state ( h 0 , u ¯ 0 , S 0 ) , when the 1-shock speed vanishes, it requires that u ¯ 0 > 0 and

(3.12) u ¯ 0 2 = h 2 S 2 h h 0 + h 0 h + 1 + g 2 h 2 1 h 0 + 1 h = h h 0 S 2 h 2 + S 0 2 h 0 2 + S 2 h h 0 + g 2 h 3 + g 2 h 2 h 0 > g h 0 + 3 S 0 2 h 0 2 = c 0 2 ,

where h > h 0 is used here. We thus conclude that the left-hand state ( h 0 , u ¯ 0 , S 0 ) belongs to D 1 . Similarly, for a 3-shock wave from a given right-hand state to a left-hand state, the shock speed vanishes if the given state belongs to D 4 . Consequently, we have the following proposition.

Proposition 3.1

For the 1-shock wave connecting the left-hand state U 0 = ( h 0 , u ¯ 0 , S 0 , B 0 ) and the right-hand state U = ( h , u ¯ , S , B ) , the shock speed may change its sign along the 1-shock curve. To be precise,
1. If U 0 ∈ D 2 ∪ D 3 ∪ D 4 , then the shock speed σ ( U 0 , U ) remains negative.
2. If U 0 ∈ D 1 , then the shock speed σ ( U 0 , U ) vanishes at some point denoted by U ˜ 0 ∈ D 2 on the 1-shock curve such that
  (3.13) σ ( U 0 , U ˜ 0 ) = 0 , σ ( U 0 , U ) > 0 , h ∈ ( h 0 , h ˜ ) , σ ( U 0 , U ˜ 0 ) < 0 , h > h ˜ .
For the 3-shock wave connecting the left-hand state U 0 = ( h 0 , u ¯ 0 , S 0 , B 0 ) and the right-hand state U = ( h , u ¯ , S , B ) , the shock speed may change its sign along the 3-shock curve. To be precise,
1. If U 0 ∈ D 1 ∪ D 2 ∪ D 3 , then the shock speed σ ( U 0 , U ) remains positive.
2. If U 0 ∈ D 4 , then the shock speed σ ( U 0 , U ) vanishes at some point denoted by U ˜ 0 ∈ D 3 on the 3-shock curve such that
  (3.14) σ ( U 0 , U ˜ 0 ) = 0 , σ ( U 0 , U ) < 0 , h ∈ ( h 0 , h ˜ ) , σ ( U 0 , U ˜ 0 ) > 0 , h > h ˜ .

3.3 Stationary wave

We first look at the stationary contact discontinuity with zero speed. In other words, we look for the time independent solution of (2.5), which is equivalent to

(3.15) ( h u ¯ ) ′ = 0 , h ( u ¯ 2 + u ˆ 2 ) + g 2 h 2 ′ = − g h B ′ ( x ) , ( u ˆ u ¯ ) ′ = 0 ,

where ( ) ′ means d d x . The solution of the ordinary differential equation is given in the following lemma.

Lemma 3.2

For smooth solutions, system (3.15) is equivalent to

(3.16) h u ¯ = h 0 u ¯ 0 , u ¯ 2 2 + 3 2 S 2 h 2 + g ( h + B ) = u ¯ 0 2 2 + 3 2 S 0 2 h 0 2 + g ( h 0 + B 0 ) , S = S 0 ,

where ( h 0 , u ¯ 0 , S 0 , B 0 ) and ( h , u ¯ , S , B ) represent the left- and right-hand state of the stationary wave, respectively.

Proof

For a given left-hand state ( h 0 , u ¯ 0 , S 0 , B 0 ) , the first equation of (3.15) is equivalent to

(3.17) h u ¯ = h 0 u ¯ 0 ,

and the third equations is equivalent to

(3.18) S = S 0 .

For the second equation, one can substitute u ˆ = S h into it and rewrite as

(3.19) h u ¯ 2 + S 2 h 3 + g 2 h 2 ′ = − g h B ′ ( x ) ,

which is equivalent to

(3.20) h u ¯ u ¯ ′ + 3 S 2 h 2 h ′ + g h h ′ = − g h B ′ ( x ) ,

where (3.17) and (3.18) are used. After eliminating h in (3.20), we finally obtain

(3.21) u ¯ 2 2 + 3 2 S 2 h 2 + g ( h + B ) ′ = 0 .

Thus, we have proved the lemma.□

We are now ready to solve system (3.16) for a given left-hand state ( h 0 , u ¯ 0 , S 0 , B 0 ) . To achieve this goal, we substitute u ¯ = h 0 u ¯ 0 h and S = S 0 into the second equation of (3.16), which yields the following function:

(3.22) Φ ( h , B ) = 3 2 S 0 2 h 4 + g h 3 + g ( B − B 0 ) − 1 2 u ¯ 0 2 − 3 2 S 0 2 h 0 2 − g h 0 h 2 + 1 2 ( h 0 u ¯ 0 ) 2 ,

where Φ ( h , B ) depends on h and the parameter B .

In order to find the roots of (3.22), we have

(3.23) d Φ ( h , B ) d h = 6 S 0 h 3 + 3 g h 2 + 2 h g ( B − B 0 ) − 1 2 u ¯ 0 2 − 3 2 S 0 2 h 0 2 − g h 0 .

So that

(3.24) d Φ ( h , B ) d h < 0 , h < h min , d Φ ( h , B ) d h > 0 , h > h min ,

where

(3.25) h min = − 3 g + 9 g 2 − 48 S 0 2 g ( B − B 0 ) − 1 2 u ¯ 0 2 − 3 2 S 0 2 h 0 2 − g h 0 12 S 0 2 ,

provided that

(3.26) g ( B − B 0 ) − 1 2 u ¯ 0 2 − 3 2 S 0 2 h 0 2 − g h 0 < 0 .

As a result, from (3.22), we know that Φ ( h = 0 , B ) > 0 , thus (3.22) has a root if and only if

(3.27) Φ ( h min , B ) ≤ 0 .

Besides, if Φ ( h min , B ) < 0 , then there are exactly two roots denoted by h * < h min < h * such that

(3.28) Φ ( h * , B ) = Φ ( h * , B ) = 0 .

Moreover, we have the following lemma.

Lemma 3.3

There exists a stationary contact discontinuity from a given state U 0 = ( h 0 , u ¯ 0 , S 0 , B 0 ) connecting the state U = ( h , u ¯ , S , B ) if and only if Φ ( h min , B ) ≤ 0 . More precisely,

If Φ ( h min , B ) > 0 , there are no stationary contacts.
If Φ ( h min , B ) ≤ 0 , there are exactly two points U * = ( h * , u ¯ * , S = S 0 , B ) and U * = ( h * , u ¯ * , S = S 0 , B ) that satisfying
(3.29) Φ ( h * , B ) = Φ ( h * , B ) = 0 ,
with U * belonging to the supercritical regions and U * belonging to the subcritical regions.

Proof

From the above discussion, here we only prove that u * > c * and u * < c * with c defined in (2.11). From (3.23), one has

(3.30) 6 S 0 h min 2 + 3 g h min + 2 g ( B − B 0 ) − 1 2 u ¯ 0 2 − 3 2 S 0 2 h 0 2 − g h 0 = 0 .

By using (3.16), we obtain that

(3.31) 6 S 0 h min 2 + 3 g h min = u ¯ * 2 + 3 S 0 2 h * 2 + 2 g h * .

Thus, it holds

(3.32) u ¯ * 2 = 6 S 0 h min 2 + 3 g h min − 3 S 0 2 h * 2 − 2 g h * > 3 S 0 2 h * 2 + g h * = c * 2 ,

where h min > h * is used. Similarly, we can prove that ( u ¯ * ) 2 < ( c * ) 2 . Consequently, we have proved the lemma.□

As we have seen in Lemma 3.3 that for a given state U 0 on one side of the stationary contact, there may admit up to a one-parameter family of solutions. In order to ensure the uniqueness of the solutions, we follow [18,30] and impose the global entropy condition on the stationary contact discontinuity.

Global entropy condition. Along the stationary contact curve S 0 ( U 0 , U ) in the ( u ¯ , h ) plane, the bottom function B obtained from (3.16) is a monotone function of h .

Under the global entropy condition, the stationary contact discontinuity can be called as the stationary wave. Moreover, we have the following lemma.

Lemma 3.4

The global entropy condition is equivalent to the statement that any stationary wave has to remain in the closure of the same domain D i , i = 1 , 2 , 3 , 4 . Moreover, the stationary wave can be viewed as a parameterized curve denoted by S 0 ( U 0 , U ( a ) ) , which is strictly increasing (decreasing) in u ¯ if u ¯ < 0 ( u ¯ > 0 ).

Proof

To prove this lemma, for a given state U 0 = ( h 0 , u ¯ 0 , S 0 , B 0 ) , we differentiate (3.16) on both sides and have that

(3.33) h d u ¯ + u ¯ d h = 0 , u ¯ d u ¯ + 3 S 2 h d h + g d h + g h d B = 0 .

We insert the first equation into the second one and have that

(3.34) d u ¯ d B = − g u ¯ u ¯ 2 − c 2 , d h d B = g h u ¯ 2 − c 2 , d u ¯ d h = − u ¯ h .

As we assume that the bottom B is a monotonic decreasing function, thus d B < 0 . From (3.34), if u ¯ 0 2 > c 0 2 , then u ¯ is increasing and h is decreasing, i.e., u ¯ > u ¯ 0 and h < h 0 . Thus, it has u ¯ 2 > c 2 . Similarly, if u ¯ 0 2 < c 0 2 , we derive that u ¯ 2 < c 2 . We thus complete the proof.□

4 Riemann problem for (2.5)

In this section, we establish the existence of the Riemann solution. Some notations are used and given below for convenience.

We use W k ( U i , U j ) to denote the k th wave (namely, the k th shock wave or the k th rarefaction wave) that connecting the left-hand state U i to the right-hand state U j .
We use W m ( U i , U j ) ⊕ W n ( U j , U k ) to denote the solution including an m th wave from U i to U j , then followed by an n th wave from U j to U k .
The notation U ¯ is used to denote the state obtained by a contact discontinuity from a given state U , which is associated with the second characteristic field λ 2 = u ¯ .

4.1 Construction of the Riemann solutions

In order to construct the Riemann solutions, we first classify the initial data by the following cases and then discuss the Riemann solutions case by case.

Case 1 U L ∈ D 1 , U R ∈ D 1 ∪ Γ + and some part of D 2 . Denote U L * ∈ D 1 as the state obtained from U L by a stationary jump from B L to B R . Then, if W 1 ( U L * ) ∩ W 3 ( U R ) ≠ ∅ , then there exists a solution. Denote

(4.1) { U 1 } = W 1 ( U L * ) ∩ W 3 ( U R ) .

Thus, the solution is given by

(4.2) W 4 ( U L , U L * ) ⊕ W 1 ( U L * , U 1 ) ⊕ J ( U 1 , U ¯ 1 ) ⊕ W 3 ( U ¯ 1 , U R ) .

See Figure 1. In particular, if W 1 ( U L * ) ∩ W 3 ( U R ) = ∅ , then there is a vacuum in this case,

(4.3) W 4 ( U L , U L * ) ⊕ W 1 ( U L * , N ) ⊕ V acuum ⊕ W 3 ( N , U R ) .

$Figure 1 Case 1. U L ∈ D 1 {U}_{L}\in {D}_{1} , U R ∈ D 1 ∪ Γ + {U}_{R}\in {D}_{1}\cup {\Gamma }_{+} and some part of D 2 {D}_{2} .$

Figure 1

Case 1. U L ∈ D 1 , U R ∈ D 1 ∪ Γ + and some part of D 2 .

Case 2 U L ∈ D 1 , U R ∈ D 2 ∪ D 3 ∪ Γ ± and some part of D 1 . First, denote U ˜ L ∈ D 2 as the state that the 1-shock speed vanishes, i.e., σ ( U L , U ˜ L ) = 0 . U ˜ L is the solution of the following equation:

(4.4) g 2 h 4 + S L 2 + g 2 h L h 3 + S L 2 h L 2 h 2 + S L 2 h L 2 h − u ¯ L 2 h L = 0 .

From each point U on the backward shock curve S 1 ( U L , U ) with p > p ˜ L , U will jump to a state U * ∈ D 2 by a stationary wave from B L to B R . Such state U * can be formed to a curve which we denote as W 1 B ( U L ) ,

(4.5) W 1 B ( U L ) ≔ { U * ∣ U * ∈ W 4 ( U , U * ) from B L to B R with U ∈ W 1 ( U L ) , p ≥ p ˜ L } .

See Figure 2. If W 3 ( U R ) ∩ W 1 B ( U L ) ≠ ∅ , denote as U 2 * = W 3 ( U R ) ∩ W 1 B ( U L ) , where U 2 * is jumped from U 2 by a stationary wave, then there exists a Riemann solution. Furthermore, if u 2 ≥ 0 , the solution is

(4.6) S 1 ( U L , U 2 ) ⊕ W 4 ( U 2 , U 2 * ) ⊕ J ( U 2 * , U ¯ 2 * ) ⊕ W 3 ( U ¯ 2 * , U R ) .

If u 2 < 0 , then the solution is

(4.7) S 1 ( U L , U 2 ) ⊕ J ( U 2 , U ¯ 2 ) ⊕ W 4 ( U ¯ 2 , U ¯ 2 * ) ⊕ W 3 ( U ¯ 2 * , U R ) .

$Figure 2 Case 2. U L ∈ D 1 {U}_{L}\in {D}_{1} , U R ∈ D 2 ∪ D 3 ∪ Γ ± {U}_{R}\in {D}_{2}\cup {D}_{3}\cup {\Gamma }_{\pm } and some part of D 1 {D}_{1} .$

Figure 2

Case 2. U L ∈ D 1 , U R ∈ D 2 ∪ D 3 ∪ Γ ± and some part of D 1 .

Case 3 We analyze the coalescence of waves in this case. Namely, there are three waves with the same zero speeds coincide with each other in the Riemann solution. To be precise, we construct the Riemann solution as follows. U L first jumps to the state U L m with the bottom function B ( x ) shifting from B L to an intermediate state B ∈ [ B L , B R ] . Then, U L m jumps to U 3 ∈ D 2 by a standing shock wave with σ ( U L m , U 3 ) = 0 . Finally, the solution follows by another stationary wave jumping from U 3 to U 4 with the bottom function shifting from B to B R (Figure 3). In fact, the three waves have the same zero speeds and thus the waves coincide in the physical plane. Now, we are left to determine the middle bottom function B , which varies from B L to B R . It is natural to denote the following curve:

(4.8) L ( U L , B L , B R ) = { U ( B ) ∣ B ∈ [ B L , B R ] } .

Thus, if

(4.9) W 3 ( U R ) ∩ L ≠ ∅ ,

we can construct the solution as follows. The solution begins from a stationary from U L to U L m , follows by a standing shock wave from U L m to U 3 , then follows by another stationary wave to reach U 4 , U 4 jumps to U ¯ 4 by a contact discontinuity, and finally the solution gets to U R by a forward wave W 3 . Precisely,

(4.10) S 0 ( U L , U L m ) ⊕ S 1 ( U L m , U 3 ) ⊕ S 0 ( U 3 , U 4 ) ⊕ J ( U 4 , U ¯ 4 ) ⊕ W 3 ( U ¯ 4 , U R ) .

Figure 3

Case 3. The coalescence of waves.

Case 4 This case holds when U L belongs to the whole plane and U R ∈ D 4 . We may construct the solution from the right-hand side to the left-hand side. U R jumps to a state U R * ∈ D 4 by a stationary wave with the bottom function shifting from B R to B L , then if

(4.11) W 3 ( U R * ) ∩ W 1 ( U L ) ≠ ∅ ,

we construct the Riemann solution as follows, Figure 4:

(4.12) S 1 ( U L , U 5 ) ⊕ J ( U 5 , U ¯ 5 ) ⊕ S 0 ( U ¯ 5 , U R * ) ⊕ W 3 ( U R * , U R ) .

$Figure 4 Case 4. U R ∈ D 4 {U}_{R}\in {D}_{4} .$

Figure 4

Case 4. U R ∈ D 4 .

Case 5 This case holds when U R belongs to some part of D 4 , and there also admits the coalescence of waves (Figure 5). Precisely, the solution begins from U L to a state U 6 by a backward shock wave S 1 , then followed by a contact discontinuity from U 6 to U ¯ 6 , followed by a stationary wave from U ¯ 6 to U 7 ∈ D 3 with the bottom function shift from B L to an intermediate state B ∈ [ B L , B R ] , then followed by a standing forward shock wave S 3 from U 7 to U R m ∈ D 4 . Finally, U R m jumps to U R by a stationary wave with B ( x ) shifting from B L to B R , i.e.,

(4.13) S 1 ( U L , U 6 ) ⊕ J ( U 6 , U ¯ 6 ) ⊕ S 0 ( U ¯ 6 , U 7 ) ⊕ S 1 ( U 7 , U R m ) ⊕ S 0 ( U R m , U R ) .

$Figure 5 Case 5. U R {U}_{R} belongs to some part of D 4 {D}_{4} .$

Figure 5

Case 5. U R belongs to some part of D 4 .

Case 6 U L ∈ D 2 ∪ D 3 ∪ D 4 ∪ Γ ± , and U R ∈ D 1 and some part of D 2 . In this case, the solution begins with a backward rarefaction wave R 1 from U L to a sonic state U C ∈ Γ + , then U C jumps to U C * ∈ D 1 , as far as W 1 ( U C * ) ∩ W 3 ( U R ) ≠ ∅ , we have a Riemann solution here (Figure 6). Denote { U 8 } = W 1 ( U C * ) ∩ W 3 ( U R ) . If p 8 ≤ p ˜ c * with σ ( U C * , U ˜ C * ) = 0 , such that the backward wave W 1 ( U C * , U 8 ) has a positive speed, then the Riemann solution is

(4.14) R 1 ( U L , U C ) ⊕ S 0 ( U C , U C * ) ⊕ W 1 ( U C * , U 8 ) ⊕ J ( U 8 , U ¯ 8 ) ⊕ W 3 ( U ¯ 8 , U R ) .

$Figure 6 Case 6. U L ∈ D 2 ∪ D 3 ∪ D 4 ∪ Γ ± {U}_{L}\in {D}_{2}\cup {D}_{3}\cup {D}_{4}\cup {\Gamma }_{\pm } , and U R ∈ D 1 {U}_{R}\in {D}_{1} and some part of D 2 {D}_{2} .$

Figure 6

Case 6. U L ∈ D 2 ∪ D 3 ∪ D 4 ∪ Γ ± , and U R ∈ D 1 and some part of D 2 .

Case 7 U R ∈ D 2 ∪ D 3 ∪ Γ ± and some part of D 1 . In this case, when U L touches U C ∈ Γ + through a rarefaction wave, we use the subcritical point U C * ∈ D 2 instead of the U C * in case 6. Similarly, to construct the Riemann solution, it is convenient to define a stationary curve S 0 ( U , U * ) with U belonging to W 1 ( U L ) , p ≥ p c . More precisely,

(4.15) W 1 B ( U L ) ≔ { U * ∣ U * ∈ W 4 ( U , U * ) from B L to B R with U ∈ W 1 ( U L ) , p ≥ p c } .

As far as W 1 ( U L B ) ∩ W 3 ( U R ) ≠ ∅ , we have a Riemann solution (Figure 7). Denote { U 9 * } = W 1 ( U L B ) ∩ W 3 ( U R ) , where U 9 * is jumped by the U 9 ∈ W 1 ( U L ) by a stationary wave. Then, the Riemann solution is

(4.16) W 1 ( U L , U 9 ) ⊕ J ( U 9 , U ¯ 9 ) ⊕ S 0 ( U ¯ 9 , U ¯ 9 * ) ⊕ W 3 ( U ¯ 9 * , U R ) .

$Figure 7 Case 7. U R ∈ D 2 ∪ D 3 ∪ Γ ± {U}_{R}\in {D}_{2}\cup {D}_{3}\cup {\Gamma }_{\pm } and some part of D 1 {D}_{1} .$

Figure 7

Case 7. U R ∈ D 2 ∪ D 3 ∪ Γ ± and some part of D 1 .

Case 8 This construction is similar to case 3. We also consider the coalescence of waves, and three waves coincide with the same zero speeds (Figure 8). We can compare it with case 6, when U L gets to U C on Γ + , instead of jumping to U C * ∈ D 1 with the bottom function shifting from B L to B R , the solution jumps to U C m ∈ D 1 with B ( x ) shifting from B L to an intermediate state B ∈ [ B L , B R ] , then U C m ∈ D 1 jumps to U 10 by a standing shock wave, followed by another stationary wave from U 10 to U 10 * with the bottom function shifting from B to B R . In order to determine the middle bottom function B , it is similar to denote the curve

(4.17) L ( U L , B L , B R ) = { U ( B ) ∣ B ∈ [ B L , B R ] } .

Thus, if

(4.18) W 3 ( U R ) ∩ L ≠ ∅ ,

we have a Riemann solution here. The solution begins from a rarefaction wave from U L to U C , followed by a stationary from U C to U C m , followed by a standing shock wave from U C m to U 10 , then followed by another stationary wave to arrive at U 10 * , U 10 * jumps to U ¯ 10 * by a contact discontinuity, and finally the solution gets to U R by a forward wave W 3 , precisely,

(4.19) R 1 ( U L , U C ) ⊕ S 0 ( U C , U C m ) ⊕ S 1 ( U C m , U 10 ) ⊕ S 0 ( U 10 , U 10 * ) ⊕ J ( U 10 , U ¯ 10 * ) ⊕ W 3 ( U ¯ 10 * , U R ) .

$Figure 8 Case 8. U L ∈ D 2 ∪ D 3 {U}_{L}\in {D}_{2}\cup {D}_{3} , U R ∈ D 1 ∪ D 2 ∪ D 3 {U}_{R}\in {D}_{1}\cup {D}_{2}\cup {D}_{3} and some part of D 4 {D}_{4} .$

Figure 8

Case 8. U L ∈ D 2 ∪ D 3 , U R ∈ D 1 ∪ D 2 ∪ D 3 and some part of D 4 .

Case 9 U L belongs to D 4 and some part of D 3 , U R ∈ D 2 ∪ D 3 . In this case, we construct the Riemann solution as follows. Define a curve S 0 ( U , U * ) formed by a stationary jump from U ∈ Γ − to U * ∈ D 3 . If S 0 ( U , U * ) ∩ W 3 ( U R ) ≠ ∅ , then denote { U 12 * } = S 0 ( U , U * ) ∩ W 3 ( U R ) , which is jumped from U 12 ∈ Γ − by a stationary wave. Whenever W 1 ( U L ) ∩ W 3 ( U 12 ) ≠ ∅ , there is a Riemann solution (Figure 9),

(4.20) W 1 ( U L , U 13 ) ⊕ J ( U 13 , U ¯ 13 ) ⊕ R 1 ( U ¯ 13 , U 12 ) ⊕ S 0 ( U 12 , U 12 * ) ⊕ W 3 ( U 12 * , U R ) .

$Figure 9 Case 9. U L {U}_{L} belongs to D 4 {D}_{4} and some part of D 3 {D}_{3} , U R ∈ D 2 ∪ D 3 {U}_{R}\in {D}_{2}\cup {D}_{3} .$

Figure 9

Case 9. U L belongs to D 4 and some part of D 3 , U R ∈ D 2 ∪ D 3 .

In summary, we have established the model describing the two-layer shallow water equations with the bottom topography. We further consider the elementary waves of the model under some proper assumptions. The Riemann solution is then solved explicitly by classifying the initial data case by case. The resonance phenomenon and coalescence of waves are also dealt with by using the phase plane method. The Riemann solution in the 1D case is appealing at least for the following two reasons. On the one hand, it can be used to design the Riemann solver for numerical scheme. On the other hand, in many practical applications, we need to extend the model to the 2D or 3D case. When we consider the Riemann problem for 2D case with initial data being constant in each quadrant [31], the elementary waves are quite similar to that in 1D case, and our results can also be applied there, this will be left for future consideration.

Acknowledgements

The authors are very grateful for the anonymous reviewers’ valuable comments that improved the manuscript.

Funding information: This work is supported by NSFC (12171305, 12201323, 12171130, 12071106, 12272059) and Zhejiang Provincial Natural Science Foundation of China (LQ22A010001).
Author contributions: Fan Wang: Writing, Editing, and Reviewing. Wancheng Sheng: Writing, Editing, and Supervision. Yanbo Hu: Writing, Editing, and Supervision. Qinglong Zhang: Writing, Editing, and Reviewing.
Conflict of interest: The authors state no conflict of interest.

References

[1] R. Abgrall and S. Karni, Two-layer shallow water systems: A relaxation approach, SIAM J. Sci. Comput. 31 (2009), 1603–1627. 10.1137/06067167XSearch in Google Scholar

[2] N. Aguillon, E. Audusse, E. Godlewski, and M. Parisot, Analysis of the Riemann problem for a shallow water model with two velocities, SIAM J. Math. Anal. 5 (2018), 4861–4888. 10.1137/17M1152887Search in Google Scholar

[3] N. Andrianov and G. Warnecke, The Riemann problem for the Baer-Nunziato two-phase flow model, J. Comput. Phys. 195 (2004), 434–464. 10.1016/j.jcp.2003.10.006Search in Google Scholar

[4] E. Audusse, A multilayer Saint-Venant model: Derivation and numerical validation, Discrete Contin. Dyn. Syst. Ser. B 5 (2005), 189–214. 10.3934/dcdsb.2005.5.189Search in Google Scholar

[5] E. Audusse, M. Bristeau, and A. Decoene, Numerical simulations of 3d free surface flows by a multilayer Saint-Venant model, Int. J. Numer. Methods Fluids 56 (2008), 331–350. 10.1002/fld.1534Search in Google Scholar

[6] E. Audusse, M. Bristeau, B. Perthame, and J. Sainte-Marie, A multilayer Saint-Venant system with mass exchanges for shallow water flows. Derivation and numerical validation, ESAIM: Math. Model. Numer. Anal. 45 (2011), 169–200. 10.1051/m2an/2010036Search in Google Scholar

[7] B. Barros, Conservation laws for one-dimensional shallow water models for one and two-layer flows, Math. Models Methods Appl. Sci. 16 (2006), 119–137. 10.1142/S0218202506001078Search in Google Scholar

[8] M. Ben-Artzi and J. Q. Li, Hyperbolic balance laws: Riemann invariants and the generalized Riemann problem, Numer. Math. 106 (2007), 369–425. 10.1007/s00211-007-0069-ySearch in Google Scholar

[9] F. Berger and J. F. Colombeau, Numerical solutions of one-pressure models in multifluid flows, SIAM J. Numer. Anal. 32 (1995), 1139–1154. 10.1137/0732052Search in Google Scholar

[10] F. Bouchut and T. Morales de Luna, An entropy satisfying scheme for two-layer shallow water equations with uncoupled treatment, M2AN Math. Model. Numer. Anal. 42 (2008), 683–698. 10.1051/m2an:2008019Search in Google Scholar

[11] M. O. Bristeau, B. Di-Martino, C. Guichard, and J. Sainte-Marie, Layer-averaged Euler and Navier-Stokes equations, Commun. Math. Sci. 15 (2017), 1221–1246. 10.4310/CMS.2017.v15.n5.a3Search in Google Scholar

[12] M. Castro, J. Maciás, and C. Parés, A Q-scheme for a class of system of coupled conservation laws with source term. Application to a two-layer 1D shallow water system, M2AN Math. Model. Numer. Anal. 35 (2001), 107–127. 10.1051/m2an:2001108Search in Google Scholar

[13] D. M. Causon, D. M. Ingram, and C. G. Minghan, A Cartesian cut cell method for shallow water flows with moving boundaries, Adv. Water Resources 24 (2001), 899–911. 10.1016/S0309-1708(01)00010-0Search in Google Scholar

[14] A. Chiapolino and R. Saurel, Models and methods for two-layer shallow water flows, J. Comput. Phys. 371 (2018), 1043–1066. 10.1016/j.jcp.2018.05.034Search in Google Scholar

[15] S. Gavrilyuk, K. Ivanova, and N. Favrie, Multi-dimensional shear shallow water flows: Problems and solutions, J. Comput. Phys. 366 (2018), 252–280. 10.1016/j.jcp.2018.04.011Search in Google Scholar

[16] A. Kurganov and G. Petrova, Central-upwind schemes for two-layer shallow water equations, SIAM J. Sci. Comput. 31 (2009), 1742–1773. 10.1137/080719091Search in Google Scholar

[17] J. Li and G. Chen, The generalized Riemann problem method for the shallow water equations with bottom topography, Int. J. Numer. Methods Eng. 65 (2006), 834–862. 10.1002/nme.1471Search in Google Scholar

[18] P. G. LeFloch and M. D. Thanh, The Riemann problem for fluid flows in a nozzle with discontinuous cross-section, Commun. Math. Sci. 1 (2003), 763–797. 10.4310/CMS.2003.v1.n4.a6Search in Google Scholar

[19] P. G. LeFloch, M. D. Thanh, and A Godunov-type method for shallow water equations with discontinuous topography in the resonant regime, J. Comput. Phys. 230 (2011), 7631–7660. 10.1016/j.jcp.2011.06.017Search in Google Scholar

[20] R. R. Long, Long waves in a two-fluid system, J. Met. 13 (1956), 70–74. 10.1175/1520-0469(1956)013<0070:LWIATF>2.0.CO;2Search in Google Scholar

[21] G. Dal Maso, P. G. LeFloch, and F. Murat, Definition and weak stability of nonconservative products, J. Math. Pure Appl. 74 (1995), 483–548. Search in Google Scholar

[22] R. Monjarret, Local well-posedness of the two-layer shallow water model with free surface, SIAM J. Appl. Math. 75 (2015), 2311–2332. 10.1137/140957020Search in Google Scholar

[23] P. J. Montgomery and T. G. Moodie, On the number of conserved quantities for the two-layer shallow water equations, Stud. Appl. Math. 106 (2001), 229–259. 10.1111/1467-9590.00166Search in Google Scholar

[24] L. Ovsyannikov, Two-layer “shallow wate” model, J. Appl. Mech. Tech. Phys. 20 (1979), 127–135. 10.1007/BF00910010Search in Google Scholar

[25] M. Petcu and R. Temam, An interface problem: the two-layer shallow water equations, Discrete Contin. Dyn. Syst. 33 (2013), 5327–5345. 10.3934/dcds.2013.33.5327Search in Google Scholar

[26] L. Qian, D. M. Causon, C. G. Mingham, and D. M. Ingram, A free-surface capturing method for two fluid flows with moving bodies, Proc. R. Soc. A 462 (2006), 21–42. 10.1098/rspa.2005.1528Search in Google Scholar

[27] W. C. Sheng, Q. L. Zhang, and Y. X. Zheng, The Riemann problem for a blood flow model in arteries, Commun. Comput. Phys. 27 (2020), 227–250. 10.4208/cicp.OA-2018-0220Search in Google Scholar

[28] W. C. Sheng, Q. L. Zhang, and Y. X. Zheng, A direct Eulerian GRP scheme for a blood flow model in arteries, SIAM J. Sci. Comput. 43 (2021), A1975–A1996. 10.1137/19M1284476Search in Google Scholar

[29] W. C. Sheng and Q. L. Zhang, The Riemann problem for a traffic flow model on a road with variable widths, IMA J. Appl. Math. 87 (2022), 757–785. 10.1093/imamat/hxac020Search in Google Scholar

[30] Q. L. Zhang, W. C. Sheng, and Y. X. Zheng, Interaction of the elementary waves for shallow water equations with discontinuous topography, Commun. Math. Sci. 19 (2021), 1381–1402. 10.4310/CMS.2021.v19.n5.a9Search in Google Scholar

[31] T. Zhang and Y. X. Zheng, Conjecture on the structure of solutions of the Riemann problem for two-dimensional gas dynamics systems, SIAM J. Math. Anal. 21 (1990), 593–630, 10.1137/0521032Search in Google Scholar

Received: 2023-05-01

Revised: 2023-10-20

Accepted: 2024-11-13

Published Online: 2024-12-07

This work is licensed under the Creative Commons Attribution 4.0 International License.

Articles in the same Issue

https://doi.org/10.1515/anona-2024-0058

Keywords for this article

the Riemann problem; laminar shallow water equations; Riemann invariants; the stationary wave; Rankine-Hugoniot condition

Creative Commons

BY 4.0