The Riemann problem for the shallow water equations with discontinuous topography: The wet–dry case

Accepted Manuscript The Riemann Problem for the Shallow Water Equations with Discontinuous Topography: The Wet-Dry Case

Carlos Parés, Ernesto Pimentel

PII: DOI: Reference:

S0021-9991(18)30745-9 https://doi.org/10.1016/j.jcp.2018.11.019 YJCPH 8372

To appear in:

Journal of Computational Physics

Received date: Accepted date:

9 October 2018 12 November 2018

Please cite this article in press as: C. Parés, E. Pimentel, The Riemann Problem for the Shallow Water Equations with Discontinuous Topography: The Wet-Dry Case, J. Comput. Phys. (2018), https://doi.org/10.1016/j.jcp.2018.11.019

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Highlights • • • • •

Wet-dry Riemann problems for the shallow water system are considered. 0, 1, or 2 solutions are found depending on the wet state. Solutions based on a reinterpretation are proposed in the case of 0 solutions. This analysis is useful to better understand the behaviour of the numerical methods. The correct simulation of wet-dry fronts is crucial in applications.

The Riemann Problem for the Shallow Water Equations with Discontinuous Topography: The Wet-Dry Case Carlos Par´es, Ernesto Pimentel University of M´alaga (Spain) November 15, 2018 Abstract In this paper we consider Riemann problems for the shallow water equations with discontinuous topography whose initial conditions correspond to a wet-dry front: at time t = 0 there is vacuum on the right or on the left of the step. Besides the theoretical interest of this analysis, the results may be useful to design numerical methods and/or to produce reference solutions to compare different schemes. We show that, depending on the state at the wet side, 0, 1, or 2 self-similar solutions can be constructed by composing simple waves. In problems with 0 solutions, the step acts as an obstacle for the fluid and physically meaningful solutions can be constructed by interpreting the problem as a partial Riemann problems for the homogeneous shallow water system. Some numerical results are shown where different numerical methods are compared. In particular, it is shown that, in the non-uniqueness cases, the numerical solutions can converge to one or to the other solution, what is the reason that explains the huge differences observed when different numerical methods are applied to the shallow water system with abrupt changes in the bottom. Moreover, problems with zero solutions will be reinterpreted as Partial Riemann problems for the homogeneous system what will allow us to build a physically solution. when one side of the step is wet and the other one is dry. We will specify the regions where we can find zero, one or two solutions, giving the form of the solution when it is possible and giving an alternative when it is not possible.

Keywords: Shallow Water model, well-balanced methods, finite volume methods, approximate Riemann solvers, high order methods. Acknowledgments. This research has been partially supported by the Spanish Government and FEDER through Research project MTM2015-70490-C2-1-R.

1

Introduction

We consider the shallow water system of Partial Differential Equations ⎧ ∂t h + ∂x (hu)= 0, ⎪ ⎪  ⎨ gh2 ∂t (hu) + ∂x hu2 + = gh∂x a, ⎪ 2 ⎪ ⎩ ∂t a = 0, that governs the flow of a shallow layer of fluid, with the following notation: 1

(1)

• h = h(x, t) ≥ 0 is the height of the water from the bottom to the surface; • u = u(x, t) is the depth-averaged horizontal velocity of the water; • g is the intensity of the gravitational field; • a = a(x) is the depth of the bottom from a reference level; • η = η(x, t) is the elevation of the surface of the water (h(x, t) = η(x, t) + a(x)). In the homogenous case, i.e. when the bottom is flat, the equations have the structure of a system of conservation laws and standard numerical methods can be used (see [21]). When the bottom is not flat, they have the structure of a system of balance laws and new difficulties arise. On the one hand, nontrivial stationary solutions appear that have to be correctly handled in order to be able to simulate waves that disturb a steady state: the design of well-balanced methods, i.e. methods that preserve all or some of the stationary solutions has been a very active front of research in the last years: see [9] and the references therein. On the other hand, due to the initial condition or to the interaction of the fluid with the varying topography, wet-dry fronts can develop. These fronts appear very frequently in practical applications: floods, dam-breaks, breaking of waves on beaches, etc. and the design of numerical methods handling correctly with them is challenging: see [3, 6, 7, 11, 13, 18, 22]. This article focuses on this specific difficulty. Godunov-type methods are based on the exact or approximate solution of Riemann problems at the intercells. A Riemann problem associated with (1) is a Cauchy Problem with initial condition:  (hl , ul , al ), x < 0, (h, u, a)(x, 0) = (2) (hr , ur , ar ), x > 0. (Notice that, in a Riemann problem, the bottom is assumed to be a step linking two regions where the bottom is flat). A good knowledge of the solution of the Riemann problem is then necessary to design good numerical Godunov-type methods. The goal of this work is to study the solution of Riemann problems corresponding to wet-dry fronts, i.e. Riemann problems whose initial condition is such that hl = 0 or hr = 0. Besides the theoretical interest of this analysis, the results may be useful to design numerical methods and/or to produce reference solutions to compare different schemes. As far as we know, a complete analysis of the solutions of this problem has never been carried out. The solution of the Riemann problem in the homogeneous case, including wet-dry situations, is discussed in [21]. In the case of a varying topography, different studies of the solution of the Riemann problem with hl > 0, hr > 0 have been performed: [1], [5], [20], [16] and [17], etc. In [7] the Riemann problem corresponding to a wet-dry situation over a step has been considered, but its solution has been discussed only in some cases in which the step acts as an obstacle for the fluid motion. In those cases, the Riemann problem is interpreted as a Partial Riemann problem (in the sense proposed in [14]) associated to the homogeneous system. When a is discontinuous (and this is the case for Riemann problems) the source term gh∂x a is a nonconservative product that, in general, cannot be defined in the distributions sense. Nevertheless, this term can be interpreted in the sense of measures using, for instance, the theory developed in [12], but this interpretation is not unique. The definition of weak solutions of the system depends thus on the interpretation of the nonconservative products as measures. The notion of weak solution considered in the articles where the Riemann problem has been studied is not always the same. We follow here the one adopted in [16] and [17] which is, in our view, the most natural: according to this definition, weak solutions develop stationary waves over the step across which Riemann invariants 2

are preserved. Since these stationary waves are associated to a linearly degenerate characteristic field, as it will be seen in next section, they can be understood as contact discontinuities, so that the preservation of Riemann invariants is natural. According to this definition, in [16] and [17] it is shown that the wet-wet Riemann problem over a step may have 1, 2, or 3 solutions depending on the initial data. The conclusions are different for other definitions of weak solution: for instance, in [5] 0 or 1 solutions are found. The paper has been structured as follows: in next section the simple waves of the system will be described. The solutions of the Riemann problems are built in Section 3 by composing these simple waves. Depending on the initial conditions, we find zero, one, or two solutions: the data sets leading to one or another situations are specified. Moreover, following [7], problems with zero solutions will be reinterpreted as Partial Riemann problems associated to the homogeneous system what will allow us to build a solution. Some numerical results will be shown in the Section 4, where different numerical methods are compared. In particular, the behavior of the numerical methods in the non-uniqueness cases will be studied: as it will be seen, they can converge to one or to the other, what is the reason of the huge differences discussed in [18]. Finally some conclusions are drawn.

2

Simple waves

System (1) can be written in nonconservative form as follows: ∂t u + A(u)ux = 0, ⎛

⎞ h u = ⎝ u ⎠ , a

where:

⎛

u A(u) = ⎝ g 0

(3) ⎞ h 0 u g ⎠. 0 0

The eigenvalues of A are: λ1 (u) = u − and:

⎞ h √ r1 (u) = ⎝ − gh ⎠ , 0 ⎛

gh,

λ2 (u) = u +

gh,

⎞ √h r2 (u) = ⎝ gh ⎠ , 0 ⎛

λ3 (u) = 0,

(4)

⎛

⎞ gh r3 (u) = ⎝ −gu ⎠ , u2 − gh

(5)

are associated eigenvectors. We note that λ1 (u) = λ3 (u) in the surface: C + = {(h, u, a) : u = λ2 (u) = λ3 (u) in the surface:

gh},

C − = {(h, u, a) : u = −

gh},

(6)

(7)

and λ1 (u) = λ2 (u) when h = 0. Therefore, system (1) is non-strictly hyperbolic. The characteristic fields associated to the eigenvalues λ1 and λ2 are genuinely nonlinear when h > 0, while the one associated to λ3 is linearly degenerate.

3

The surface C := C + ∪ C − contains the critical states and it divides the half-space h ≥ 0 of the (h, u, a)-space in the following regions: A1 := {u = (h, u, a) ∈ R+ × R × R : λ2 (u) > λ1 (u) > λ3 (u)}, A2 := {u = (h, u, a) ∈ R+ × R × R : λ2 (u) > λ3 (u) > λ1 (u)}, A+ u = (h, u, a) ∈ A2 : u > 0}, 2 := { A− u = (h, u, a) ∈ A2 : u < 0}, 2 := { A3 := {u = (h, u, a) ∈ R+ × R × R : λ3 (u) > λ2 (u) > λ1 (u)}

(8)

(see Figure 1). System (1) is strictly hyperbolic in the interior of these regions. Moreover, states belonging to A1 ∪ A3 are supercritical and those belonging to A2 subcritical.

Figure 1: Projection of the regions Ai , i = 1, 2, 3 on the (h, u)-plane Following [16], the simple waves for system (1) are: • Rarefaction waves, which are smooth solutions of (1) with constant value of a associated to the genuinely nonlinear fields. These waves depend only on the self-similar variable x/t, and satisfy: – Riemann invariants preservation: the Riemann invariants

u + 2 gh, u − 2 gh

(9)

are constant along 1-rarefactions and 2-rarefactions, respectively. – Divergence of the characteristics: λi (ul ) < λi (ur ), where ul and ur are the left and the right states of the i-rarefaction, i = 1, 2. 4

(10)

• Shock waves, which are discontinuous solutions of (1) with constant value of a associated to the genuinely nonlinear fields. These waves satisfy: – Rankine Hugoniot conditions: σi [h] = [hu], σi [hu] = [u2 h + gh2 /2],

(11)

where σi is the speed of the shock associated to the λi -field and [w] represents the jump at a discontinuity of the variable w. – Lax entropy conditions: λi (ul ) > σi > λi (ur ).

(12)

• Stationary contact discontinuities, which are discontinuous solutions of (1) with discontinuous value of a, associated to the linearly degenerated field corresponding to the null eigenvalue that satisfy: [hu]  2 = 0,  (13) u /2 + g(h − a) = 0. Remark 1. The first equality in (13) can be understood as the preservation of the mass-flow through stationary contact discontinuities and the second one as the preservation of the mechanical energy: u2 /2 is the kinetic energy and g(h − a) the potential one. As it has been mentioned in the Introduction, the definition of weak-solution is not unique and different choices may lead to different jump conditions for these stationary discontinuities: this is the case in [5] or [20]. The jump conditions chosen by these authors imply the dissipation of mechanical energy through stationary contact discontinuities. Given a left-hand state ul , the 1-shock S1 (ul ) and the 2-shock S2 (ul ) consisting of all right-hand states u that can be connected to ul by a shock associated with λ1 and λ2 , respectively, are:   g 1 1 (h − hl ) + , h > hl , (14) S1 (ul ) : u = ul − 2 h hl   g 1 1 (h − hl ) + , h < hl . (15) S2 (ul ) : u = ul + 2 h hl Given a right-hand state ur , the backward 1-shock curve S1B (ur ) and the backward 2-shock curve S2B (ur ), consisting of all left-hand states u that can be connected to ur by a shock associated with λ1 and λ2 , respectively, are:   g 1 1 B S1 (ur ) : u = ur − (h − hr ) + , h < hr , (16) 2 h hr   g 1 1 B (h − hr ) + , h > hr . (17) S2 (ur ) : u = ur + 2 h hr Moreover, given two states u0 and u connected by a 1-shock wave or by a 2-shock wave, the speed of that shock will be, respectively:    g h2 h+ , (18) σ1 (u0 , u) = u0 − 2 h0 5

   g h2 . h+ σ2 (u0 , u) = u0 + 2 h0

(19)

The following result, whose proof can be found in [16], gives information about the sign of the speed of shocks: ˆ 0, u Proposition 1. 1. If u0 = (h0 , u0 , a0 ) ∈ A1 , then there exists u ˆ 0 = (h ˆ0 , a0 ) ∈ S1 (u0 ) ∩ A+ 2 ˆ 0 > h0 such that: with h (a) σ1 (u0 , u ˆ0 ) = 0, ˆ 0 ), (b) σ1 (u0 , u) > 0 for all u ∈ S1 (u0 ) such that h ∈ (h0 , h ˆ 0 , +∞). (c) σ1 (u0 , u) < 0 for all u ∈ S1 (u0 ) such that h ∈ (h 2. If u0 ∈ A2 ∪ A3 , then σ1 (u0 , u) < 0 for all u ∈ S1 (u0 ). ˆ 0, u ˆ 3. If u0 = (h0 , u0 , a0 ) ∈ A3 , then there exists u ˆ 0 = (h ˆ0 , a0 ) ∈ S2B (u0 ) ∩ A− 2 with h0 > h0 such that: ˆ0 ) = 0, (a) σ2 (u0 , u ˆ 0 ), (b) σ2 (u0 , u) < 0 for all u ∈ S2B (u0 ) such that h ∈ (h0 , h ˆ 0 , +∞). (c) σ2 (u0 , u) > 0 for all u ∈ S B (u0 ) such that h ∈ (h 2

4. If u0 ∈ A1 ∪ A2 , then σ2 (u0 , u) > 0 for all u ∈ S2B (u0 ). Given a left-hand state ul , the 1-rarefaction R1 (ul ) and the 2-rarefaction R2 (ul ) consisting of all right-hand states u that can be connected to ul by a rarefaction associated with λ1 and λ2 , respectively, are:

(20) R1 (ul ) : u = ul − 2( gh − ghl ), h ≤ hl ,

R2 (ul ) : u = ul + 2( gh − ghl ), h ≥ hl . (21) Given a right-hand state ur , the backward 1-rarefaction curve R1B (ur ) and the backward 2rarefaction curve R2B (ur ), consisting of all left-hand states u that can be connected to ur by a rarefaction associated with λ1 and λ2 , respectively, are:

(22) R1B (ur ) : u = ur − 2( gh − ghr ), h ≥ hr ,

B R2 (ur ) : u = ur + 2( gh − ghr ), h ≤ hr . (23) Moreover, given two states ul and ur connected by a 1-rarefaction or by a 2-rarefaction, the speed of the head and the tail of these rarefactions are, respectively: √ SH1 = λ1 (ul ) = ul − √ghl , (24) ST1 = λ1 (ur ) = ur − ghr , √ SH2 = λ2 (ur ) = ur +√ ghr , (25) ST2 = λ2 (ul ) = ul + ghl .

6

Using again the notation in [16], we define the curves: W1 (ul )

=

S1 (ul ) ∪ R1 (ul ),

(26)

W1B (ur )

=

S1B (ur )

(27)

W2 (ul )

=

S2 (ul ) ∪ R2 (ul ),

(28)

=

S2B (ur )

(29)

W2B (ur )

∪ ∪

R1B (ur ), R2B (ur ).

The four curves can be parameterized in the form u = u(h), h > 0, where the function is: • strictly convex and strictly decreasing for W1 (ul ) and W1B (ul ); • strictly concave and strictly increasing for W2 (ul ) and W2B (ul ). An example of W1 (ul ) and W2 (ul ) is shown in Figure 2 and another of W1B (ur ) and W2B (ur ) in Figure 3.

Figure 2: W1 and W2 curves for a state ul Let us consider now the set W3 (u0 ) of all the states u that can be connected to u0 by an admissible 3-stationary wave. Using the jump conditions (13), we obtain: (h, u, a) ∈ W3 (u0 )

u=

=⇒

where ϕ : (0, ∞) → R is given by: u2 ϕ(h) = a0 − a + 0 2g



u0 h0 , h

ϕ(h) = 0,

 h20 − 1 + h − h0 . h2

(30)

(31)

Therefore, given a and a0 , in order to find the states that can be linked through a 3-stationary wave to a state u0 , one has to look for the roots of the function ϕ. The following results hold (see [16]): 7

Figure 3: W1B and W2B curves for a state ur Lemma 1. Suppose that u0 = (h0 , u0 , a0 ) and a are given with u0 = 0. Let us define:  2 2  13 u 0 h0 hmin (u0 ) = , g  2  h0 u2 amin (u0 ) = a0 + 0 − 1 + hmin − h0 . 2g h2min Then

(32) (33)

• if a > amin the function ϕ has two roots h∗ (u0 ), h∗ (u0 ) with h∗ ≤ hmin ≤ h∗ ; • if a = amin it has only one root hmin ; • if a < amin it has no roots. Proposition 2. Given a left-hand state u0 = (h0 , u0 , a0 ) and a right-hand bottom level a: 1. If u0 = 0 and a > amin (u0 ), there are two distinct right-hand states u∗ = (h∗ , u∗ , a) and u∗ = (h∗ , u∗ , a) that can be connected to u0 by a 3-stationary wave. Here, h∗ , h∗ are the roots of ϕ and h0 u0 u∗∗ = ∗ . h∗ ∗ Moreover, u∗ is subcritical and u is supercritical. 2. If u0 = 0 and a = amin (u0 ), there is only a state that can be connected to u0 by a 3-stationary wave: (hmin , umin , amin ) with u 0 h0 . umin = hmin This state is critical. 8

3. If u0 = 0 and a < amin (u0 ), there is no stationary wave from u0 to a state with level a. 4. If u0 = 0 and a ≥ a0 − h0 , there is only a state that can be connected to u0 by a 3-stationary wave: u = u0 = 0, h = h0 + a − a0 . This state is subcritical. Remark 2. According to Remark 1, items 3 and 4 of the proposition can be understood as follows: if u0 = 0 and a < amin (u0 ) or if u0 = 0 and a < a0 − h0 , the mechanical energy of the water is not enough to go up the step. Following [16] not all the possible stationary waves are considered to be admissible. The following criterion is imposed: Monotonicity criterion (MC): Along W3 (u0 ), the bottom level a is a monotone function of h. Accordingly, the 3-stationary waves connecting a supercritical and a subcritical state are not admissible. Given a state u0 = (h0 , u0 , a0 ) ∈ Ai , i = 1, 2, 3, and a level a ≥ amin (u0 ), according to (MC) there is only one state that can be connected to u0 with an admissible 3-stationary wave: as in [16], it will be represented by SW (u, a). And given a critical state u0 ∈ C and a level a > a0 there are two states that can be connected to u0 with an admissible 3-stationary wave: we will represent by SWsup (u0 , a) (resp. SWsub (u0 , a)) the supercritical (resp. subcritical) one. The following notation will be used to represent the structure of the solution of the Riemann problems: the symbol Wi1 (u0 , u1 ) ⊕ · · · ⊕ Wik (uk−1 , uk ) will indicate that the solution of the Riemann problems is composed by k simple waves: u0 is the left state; uk , the right state; and uj , j = 1, . . . , k − 1 the intermediate states. Finally, the indexes ij ∈ {1, 2, 3} indicate the characteristic field to which the jth wave is associated, i.e. uj+1 ∈ Wij (uj ),

j = 0, . . . , k − 1.

Moreover, in the case of 1 and 2-waves, the type of simple wave (rarefaction or shock) will be specified by replacing Wi by Ri or Si .

3

The wet-dry Riemann Problem

We consider the shallow water system (1) with initial conditions:  ul = (hl , ul , al ) if x < 0, u(x, 0) = ur ∈ {(0, u, ar ) : u ∈ R} if x > 0, 

or u(x, 0) =

ul ∈ {(0, u, al ) : u ∈ R} ur = (hr , ur , ar )

if x < 0, if x > 0,

(34)

(35)

with hl , hr > 0. These initial conditions correspond to a situation in which there is a step at x = 0, the bottom is flat to the left and to the right of the step, and there is no water to the left or to the right of the step. Although from the physical point of view the value of u at a dry state is meaningless, from the mathematical point of view, a dry state can be represented by any point of the plane h = 0 in the (h, u, a)-space. From this point of view, the considered problem is a Partial 9

Riemann Problem (see [14]), as only the belonging to a given set is imposed at the right or at the left of x = 0. Moreover, as it will be seen, the value of u at the side that initially is dry will be the limit of the velocity at the wet-dry front. In the case of a flat bottom (i.e. if al = ar ) the following result holds (see [21]): Theorem 1. Let us suppose that al = ar . Then, the partial Riemann problem corresponding to the homogeneous shallow water system with initial conditions (34) (resp. (35)) has a unique solution consisting of a 1-rarefaction (resp. a 2-rarefaction) connecting ul (resp. ur ) to vacuum. Without loss of generality we will consider that al < ar , i.e. the right side of the step is deeper than the left one.

3.1

Case 1: initial condition (34)

To find a solution of the Riemann problem, we have to connect the state ul to a state ur belonging to the line {h = 0} of the plane a = ar through admissible simple waves. The following result holds: Theorem 2. Suppose that the Riemann problem with initial condition (34) has a self-similar solution u composed by admissible simple waves. Let us denote by u− u+ u to 0 and  0 the limits of  + − the left and to the right of x = 0. Then, necessarily u0 ∈ A1 and u0 ∈ A1 ∪ C + . Proof. Since the source term vanishes at (0, ∞), the function  + if x < 0, u0 v (x, t) = u(x, t) if x > 0, has to be a solution of the homogeneous shallow water system. Moreover, v has to be the selfsimilar solution of the Partial Riemann problem linking u+ 0 to vacuum. Then, due to Theorem 2, − v necessarily consists of a 1-rarefaction linking u+ u+ 0 to the line h = 0. If  0 ∈ A2 ∪ A3 ∪ C , the head of the 1-rarefaction SH1 would be negative, but then the 1-rarefaction could not follow the u+ u+ stationary wave linking u− 0 and  0 in the solution of the Riemann problem. Therefore  0 belongs + to C ∪ A1 . u+ The stationary wave linking u− 0 and  0 can be a 3-stationary wave or the composition of 3+ stationary wave and stationary shocks. Due to Theorem 1, since u+ 0 ∈ C ∪ A1 , it cannot be + the right state of a stationary shock in a = ar (it should belong to A2 to be the right state of a u∗0 1-shock or to A3 to be the right state of a 2-shock). Therefore, u+ 0 has to be linked to a state  in a = al through a 3-stationary wave. Due to (MC) the only possibility for having an admissible u∗0 cannot be the right state of a stationary 3-stationary wave is u∗0 ∈ A1 ∪ C + and u+ 0 ∈ A1 . Again,  − − ∗ shock in a = al . Therefore u0 = u0 and thus u0 ∈ A1 ∪ C + and u+ 0 ∈ A1 , as we wanted to prove. We are going to show that, if al < ar , the Partial Riemann problem with initial condition (34) has always one solution. To construct it, let us consider the following three regions of the plane a = al (see Fig 4): • Region I: A1 .

√ • Region II: A2 ∪ C − ∪ A3 ∩ {(h, u) : u > −2 gh}. √ • Region III: {(h, u) : u < −2 gh}. 10

Figure 4: Regions of the plane a = al for the partial Riemann problem with initial conditions (34). Let us study the Riemann problem with initial condition (34) for a left state ul belonging to any of these 3 regions and their boundaries. • Region I: Given a state ul ∈ A1 , we first consider the 3-stationary wave that connects ul to u+ ul , ar ) ∈ A1 and then the 1-rarefaction connecting u+ ur belonging 0 = SW ( 0 to a state  to the line h = 0. Since ul , u+ 0 ∈ A1 , both the head and the tail of the 1-rarefaction are positive, so that it can follows the stationary wave. In conclusion, we have found a solution of the form: u+ ur ). (36) W3 (ul , u+ 0 ) ⊕ R1 ( 0 , Figure 5 shows the projection of the intermediate states and the simple waves on the (h, u)plane and a sketch of the free surface at a time t > 0. ul ) ∩ C + . This • Region II: Given a state ul in region II, we first consider the state u− 0 ∈ R1 ( state has to satisfy: ⎧  ⎪ ⎨ u− = gh− 0 0 ,   (37) √ − − ⎪ gh0 − ghl . ⎩ u0 = ul − 2 Some easy computations allow us to solve this system: h− 0

1 = g



√ 2 ul + 2 ghl , 3

u− 0

√ ul + 2 ghl . = 3

(38)

We consider then the 1-rarefaction linking ul to u− 0 , followed of the 3-stationary wave linking − = SW ( u , a ) ∈ A . Finally u+ ur in the this latter state to u+ sup r 1 0 0 0 is connected to a state  ∈ C + , the line h = 0 through a 1-rarefaction: see Figure 6. Since ul ∈ A2 ∪ A3 and u− 0 11

Figure 5: Solution in Region I. Left: projection of the intermediate states and the simple waves on the (h, u)-plane. Right: sketch of the free surface at a time t > 0. The direction of the movement of the simple waves is indicated with a red arrow. head of the first 1-rarefaction is negative and the tail is 0, so that it can be followed by the 3-stationary wave. In turns, since u+ ur ∈ A1 , the second 1-rarefaction has positive 0 ∈ A1 and  head and tail, so that it can follows the 3-stationary wave. We obtain thus a solution whose structure is u− u+ u+ ur ). (39) R1 (ul , u− 0 ) ⊕ W3 ( 0 , 0 ) ⊕ R1 ( 0 , • Region III: A state ul in Region √ III can be linked through a 1-rarefaction to the state ur = (0, ur , al ), where ur = ul + 2 ghl : see Figure 7. The wet-dry front travels at the speed ur < 0, so that at time t > 0 there is vacuum in x > ur t. The structure of the solution is thus R1 (ul , ur ). (40) • Bounday between Region I and Region II: For states ul ∈ C + a solution can be ul , ar ). constructed like in Region I, i.e. (36) with u+ 0 = SWsup ( √ • Boundary between Region II and Region III: For states such that {ul = −2 ghl } a solution can be constructed like in Region III, i.e. (40). In this case, there is vacuum at x > 0 for every t > 0.

3.2

Case 2: initial condition (35)

To find a solution of the Riemann problem, we have to connect the state ur to a state ul belonging to the line {h = 0} of the plane a = al through admissible simple waves. The following result holds: Theorem 3. Suppose that the Riemann problem with initial condition (35) has a self-similar solution u composed by admissible simple waves. Let us denote by u− u+ u to 0 and  0 the limits of  − + − the left and to the right of x = 0. Then, necessarily u0 ∈ A3 ∪ C and u0 ∈ A3 ∪ SWsub (C − , ar ).

12

Figure 6: Solution in Region II. Left: projection of the intermediate states and the simple waves on the (h, u)-plane. Right: sketch of the free surface at a time t > 0. The direction of the movement of the simple waves is indicated with a red arrow. Proof. Similar to the proof of Theorem 2. We are going to show that, if al < ar , the Partial Riemann problem with initial condition (35) may have zero, one, or two solutions. To see it, let us consider the following six regions of the plane a = ar (see Fig 8): √ • Region I: A1 ∩ {(h, u) : u > 2 gh}. √ • Region II: States with u > 0 that are between the curves {(h, u) : u = 2 gh} and W2 (urest ), where urest = (ar − al , 0, ar ) is the state corresponding to a situation of water at rest with ηr = −al . • Region III: States between the curves SWsub (C − , ar ) and W2 (urest ). • Region IV: States with u < 0 that are between the curves W2 (urest ) and SWsup (C − , ar ). • Region V: States between the curves SWsub (C − , ar ), SWsup (C − , ar ), and W2 (urest ). • Region VI: States below the curve SWsup (C − , ar ). Remark 3. The mechanical energy of the states belonging to Regions IV and V is not enough to go up the step (see Remark 2). Let us study the Riemann problem with initial condition (35) for a right state belonging to any of these 6 regions and their boundaries. • Region I: A state ur in Region √ I can be linked through a 2-rarefaction to the state ul = (0, ul , ar ), where ul = ur − 2 ghr : see Figure 9. The wet-dry front travels at the speed

13

Figure 7: Solution in Region III. Left: projection of the intermediate states and the simple waves on the (h, u)-plane. Right: sketch of the free surface at a time t > 0. The direction of the movement of the simple waves is indicated with a red arrow. ul > 0, so that at time t > 0 there is vacuum in x < ul t. The structure of the solution is thus R2 (ul , ur ). (41) • Region III: Let us assume that the curves R2B (ur ) and SWsub (C − , ar ) intersect at one point u+ 0 : see Figure 10. To construct a solution of the Riemann problem, we consider the 2-rarefaction linking ur to u+ 0 followed by the 3-stationary wave linking this latter state to + − = SW ( u , a ). Finally  u ul in h = 0 through a 2-rarefaction. u− l 0 0 0 is connected to a state  + , the head and tail of the first 2-rarefaction are positive, so Since ur ∈ A1 ∪ A2 and u0 ∈ A− 2 − that it can follow the 3-stationary wave. In turns, since u− and ul ∈ A3 , the second 0 ∈ C 2-rarefaction has negative head and the tail is equal to 0, so that it can be followed by the 3-stationary wave. We obtain thus a solution whose structure is R2 (ul , u− u− u+ u+ ur ). 0 ) ⊕ W3 ( 0 , 0 ) ⊕ R2 ( 0 ,

(42)

See Figure 10. To finish, let us check that R2B (ur ) and SWsub (C − , ar ) intersects at one point: Proposition 3. Suppose that al < ar and ur belong to Region III. Then, the intersection of the curves R2B (ur ) and SWsub (C − , ar ) is non-empty and consists of only one state u+ 0. B Proof. A state u+ ur ) ∩ SWsub (C − , ar ) has to satisfy: 0 ∈ R2 (   ⎧ √ ⎪ + + ⎪ u = u + 2 gh − gh ⎨ 0 r r , 0   2 2 (h+ (u+ ⎪ + + 0) 0) ⎪ = a ( u ) = a + − 1 + hmin (u+ a ⎩ l min 0 r 0 ) − h0 , 2 2g hmin (u+ ) 0

14

(43)

Figure 8: Regions of the plane a = ar for the partial Riemann problem with initial conditions (35). where: hmin (u+ 0)

 =

2 + 2 (u+ 0 ) (h0 ) g

 13

.

Therefore h+ 0 has to be a root of the function ψ(h) = ar − al +

3 2



u(h)2 h2 g

where: u(h) = ur + 2(

 13

−

u(h)2 − h, 2g

(44)

gh − ghr ).

Let us study the roots of ψ. To do this, we consider first the state urc ∈ R2B (ur ) ∩ C − . This state has to satisfy:  r √ uc = − ghc , √ (45) urc = ur + 2( ghrc − ghr ). Some easy computations show that: hrc =

1 (ur − 2 ghr )2 , 9g

1 urc = − (2 ghr − ur ). 3

(46)

Next, we consider the state ur0 = (hr0 , 0, ar ) ∈ (R2 (ur ) ∪ R2B (ur )) ∩ {u = 0}. This state has to satisfy:  r u0 = 0,

√ (47) ur0 = ur + 2( ghr0 − ghr ).

15

Figure 9: Solution for states in region I. Left: projection of the intermediate states and the simple waves on the (h, u)-plane. Right: sketch of the free surface at a time t > 0. The direction of the movement of the simple waves is indicated with a red arrow. Some easy computations show that: √ hr0 =

ghr − g

ur 2 2 .

(48)

Observe that, in the interval [hrc , hr0 ], the curve R2B (ur ) travels from C − to u = 0. Therefore, if there is an intersection with SWsub (C − , ar ), the corresponding value of h has to be in this interval. Now, since urc ∈ C − , it verifies:

and, as a consequence:

amin (urc ) = ar .

(49)

ψ(hrc ) = ar − al > 0.

(50)

On the other hand, since ur is in Region III, R2B (ur ) intersects the line u = 0 at the right of the state urest and thus: ψ(hr0 ) = ar − al − hr0 < hr0 − hr0 = 0.

(51)

Therefore, there is at least one root in [hrc , hr0 ]. Let us see that there is only one. First, we rewrite ψ as ψ(h) = f (h) − p(h) with f (h) = ar − al + and: p(h) =

16

3 2



u(h)2 h2 g

u(h)2 + h. 2g

 13

,

(52)

(53)

Figure 10: Solution of a state in Region III. Left: projection of the intermediate states and the simple waves on the (h, u)-plane. Right: sketch of the free surface at a time t > 0. The direction of the movement of the simple waves is indicated with a red arrow. The derivatives of these two functions are: f  (h) =

1 u (h)h + u(h) 1

g3

p (h) = 

1

1

u(h) 3 h 3

u(h)u (h) + 1, g

,

(54)

(55)

g is the derivative of u(h). The denominator of f  is zero when u(h) = 0 h r r r is differentiable and and this only happens when √ h = h0 . Therefore, in (hc , h0 ) the function its derivative vanishes if − gh = u(h), what only happens for h = hrc . So, the sign of f  has to be constant in (hrc , hr0 ). It is easy to see that: 

where u (h) =

lim

h→(hr0 )−

f  (h) = −∞,

(56)

so that f  is negative in (hrc , hr0 ), i.e., f is strictly decreasing in this interval. On the other hand, p (h) only vanishes at hc , so that the sign of p remains constant in (hrc , hr0 ). It is easy to see that: (57) p (hr0 ) = 1 > 0, so that p is positive in (hrc , hr0 ), i.e., p is strictly increasing in the interval. Since f is strictly decreasing and p is strictly increasing in (hrc , hr0 ), there is only a root of ψ in this interval, as we wanted to prove. • Region V: Let us assume that the curves S2B (ur ) and SWsub (C − , ar ) intersect at one point u+ 0 : see Figure 11. To construct the solution of the Riemann problem, we consider the 2u− shock linking ur to u+ 0 followed by the 3-stationary wave linking this latter state to  0 = + − SW (u0 , al ). Finally u0 is connected to a state ul in h = 0 through a 2-rarefaction: see 17

− Figure 11. Since u− and ul ∈ A3 , the 2-rarefaction has negative head and the tail 0 ∈ C is 0, so that it can be followed by the 3-stationary wave. We obtain thus a solution whose structure is u− u+ u+ ur ), (58) R2 (u0 , u− 0 ) ⊕ W3 ( 0 , 0 ) ⊕ S2 ( 0 ,

provided that the speed of the 2-shock linking ur to u+ 0 is positive, so that the shock can follow the 3-stationary wave. This is true for the states of Region V belonging to A2 ∪ C − because of Proposition 1 but it has to be proved for those belonging to A3 . Let us prove first that u+ 0 exists and is unique: Proposition 4. Suppose that al < ar and ur in Region V or Region VI. Then, the intersection of the curves S2B (ur ) and SWsub (C − , ar ) is non-empty and consists of only one state u+ 0. B ur ) ∩ SWsub (C − , ar ) has to satisfy: Proof. A state u+ 0 ∈ S2 (  ⎧  ⎪ g 1 1 ⎪ ⎪ u+ = ur + (h+ − hr ) , ⎪ + + ⎪ 2 0 h ⎨ 0 h0 r

⎪   ⎪ 2 2 ⎪ (h+ (u+ ⎪ + + 0) 0) ⎪ − 1 + hmin (u+ ⎩ al = amin (u0 ) = ar + 0 ) − h0 . 2 2g hmin (u+ 0)

(59)

Then, h+ 0 has to be a root of the function 3 φ(h) = ar − al + 2





where: u(h) = ur +

u(h)2 h2 g

 13

g (h − hr ) 2



−

u(h)2 − h, 2g

(60)

1 1 + . h hr

From this point, the proof is similar to that of Theorem 3, replacing ur0 and urc by us0 ∈ S2B (ur ) ∩ {u = 0} and usc ∈ (S2 (ur ) ∪ S2B (ur )) ∩ C − . The positiveness of the speed of the 2-shock linking ur to u+ 0 for states of Region V belonging to A3 is a consequence of the following result: Proposition 5. Suppose al < ar and that the states ur = (hr , ur , ar ) ∈ A3 and u ˆ0 = ˆ 0, u (h ˆ0 , ar ) ∈ S2B (ur ) can be linked by a stationary shock. Then amin (ur ) < amin (u ˆ0 ). Proof. Since the shock linking u ˆ0 and ur is stationary, we can express u ˆ0 in terms of ur :    2 h 8u r r ˆ0 = −1 + 1 + , (61) h 2 ghr u ˆ0 =

18

ur h r . ˆ0 h

(62)

Figure 11: Solution of a state in Region V. Left: projection of the intermediate states and the simple waves on the (h, u)-plane. Right: sketch of the free surface at a time t > 0. The direction of the movement of the simple waves is indicated with a red arrow. ˆ 0: ˆ0 ) using (61), (62) and hr < h Let’s see that amin (ur ) < amin (u amin (ur ) < amin (u ˆ0 )

⇔ ⇔ ⇔ ⇔ ⇔ ⇔

 1  1 ˆ 0 )2 3 u0 )2 (h 3 u2r h2r 3 u2r 3 (ˆ (ˆ u0 ) 2 ˆ − hr < ar + − h0 ar + − − 2 g 2g 2 g 2g  1  1 3 u2r h2r 3 u2r 3 u2r h2r 3 (ˆ u0 ) 2 ˆ − hr < − h0 − − 2 g 2g 2 g 2g u2 (ˆ u0 )2 ˆ − h0 − r − hr < − 2g 2g u2r ˆ ˆ 0 )2 (h0 + hr ) > g(h 2       8u2r u2r hr 8u2r 8u2r h2r 1+ 1+ 2+ −2 1+ >g 4 ghr 4 ghr ghr       4u2 8u2 8u2 u2r 1+ 1+ r >2 1+ r − 1+ r . ghr ghr ghr ghr

If we define z by z2 = 1 + we must prove:

what is equivalent to

8u2r , ghr

  1 2 1 2 (z − 1)(z + 1) > 2 1 + (z − 1) − z 8 2 (z − 1)(z − 3)2 > 0

and this is true, since z > 1: ur ∈ A3 , so ur = 0. 19

Corollary 1. Suppose that al < ar and that ur ∈ A3 belongs to Region V. Then, the 2-shock linking ur to the state u+ 0 given by Proposition 4 has positive speed. Proof. Observe that, in the half-plane u ≤ 0 of the plane a = ar , the states such that amin (u) = al are those belonging to the curves SWsub (C − , ar ) and SWsup (C − , ar ). On the other hand, since amin (u) = ar > al , ∀u ∈ C − ∩ {a = ar }, in the region of a = ar between SWsub (C − , ar ) and SWsup (C − , ar ) the following inequality holds amin (u) > al , and this is the case, in particular, for ur . Using Theorem 5, we have al < amin (ur ) < amin (u ˆ0 ), where u ˆ0 is the state corresponding to the 2-stationary shock. As a consequence, u ˆ0 has to be in the region between SWsub (C − , ar ) and SWsup (C − , ar ). More specifically, since B u+ ur ) ∩ Ssub (C − , ar ), u ˆ0 has to belong to the arc of S2B (ur ) that links ur to u+ 0 = S2 ( 0 what + ˆ 0 < h . Using Proposition 1 we have that the speed of S2 (u+ , ur ) is positive. implies h 0 0 • Region VI: Given a state ur in Region VI, we first consider the 3-stationary wave that connects ur to u− ur , ar ) and then the 2-rarefaction connecting u− ul 0 = SW ( 0 to a state  belonging to the line h = 0: see Figure 12. Since ur ∈ A3 , by (MC) u− ∈ A and thus 3 0 the head and the tail of the 2-rarefaction are negative, so that it can be followed by the stationary wave. In conclusion, we have found a solution of the form: R2 (ul , u− u− ur ). 0 ) ⊕ W3 ( 0 ,

(63)

For some states of this region we can also construct a solution of the form (58). This solution is possible whenever the 2-shock has positive speed. Region VI can be split into two subregions, Regions VIa and VIb, whose boundary is the curve Γ = {u ∈ A3 : σ2 (u, u+ 0 ) = 0}, composed by the states of A3 such that the 2-shock linking them with the state u+ 0 given u = ˆ (see Figure 13). In Region VIa (including the by Proposition 4 is stationary, i.e. u+ 0 0 boundary) a solution of the form (58) can be constructed. • Regions II and IV: A self-similar solution of the Riemann problem with initial conditions (35) composed by admissible simple waves cannot be constructed in these cases: the states belonging to these regions cannot be linked through a 2-wave with positive or null velocity to a state whose mechanical energy is enough to go up the step (see Remarks 1 and 2). In other words, the step acts in these cases as an obstacle for the fluid. Therefore, in order to construct a physically meaningful solution, we consider x = 0 as a boundary of the problem and impose the boundary condition u = 0: the flow is arrested by the obstacle. To impose

20

Figure 12: Solution of a state in Region VI. Left: projection of the intermediate states and the simple waves on the (h, u)-plane. Right: sketch of the free surface at a time t > 0. The direction of the movement of the simple waves is indicated with a red arrow. properly this boundary condition, following [7], we consider the partial Riemann problem consisting of the homogeneous shallow water system with initial condition:  u(x, 0) = ur , if x > 0, (64) u = 0, if x < 0. To solve this problem, from ur we have to reach a state in the line {u = 0}. In Region II we consider the 2-rarefaction that connects ur to ur0 ∈ R2B (ur ) ∩ {u = 0} given by (48): see Fig 14. Since ur ∈ A1 ∪ A2 and ur0 ∈ A2 , the tail and the head of the 2-rarefaction are positive. The structure of the solution is thus: R2 (ur0 , ur ).

(65)

In Region IV we consider the 2-shock linking ur to us0 ∈ S2B (ur ) ∩ {u = 0}: see Fig 15. This 2-shock has positive speed. The structure of the solution is thus: S2 (us0 , ur ).

(66)

√ • Boundary between Region I and Region II: for states such that {ur = 2 ghr } a solution can be constructed like in Region I, i.e. (41). In this case, there is vacuum at x < 0 for every t > 0. • Boundary between Region II and Region III: for states ur ∈ R2 (urest ) a solution can be constructed consisting of the 2-rarefaction connecting ur to urest followed by the 3stationary wave connecting urest to the vacuum state (0, 0, al ). The structure of this solution is W3 ((0, 0, al ), urest ) ⊕ R2 (urest , ur ). • Boundary between Region II and Region IV: the partial Riemann problem for the homogeneous shallow water system with initial condition (64) in considered. The solution is stationary in this case. 21

Figure 13: The two subregions of Region VI. • Boundary between Region III and Region V: for states ur ∈ SWsub (C − , ar ) we conur , al ) and then the 2-rarefaction sider the 3-stationary wave that connects ur to u− 0 = SW ( ul belonging to the line {h = 0}. Since ur ∈ SWsub (C − , ar ), connecting u− 0 to a state  − and thus the head of the 2-rarefaction has zero speed and its tail is positive. u− 0 ∈ C Thereupon we can construct a solution of the form: u− ur ). R2 (ul , u− 0 ) ⊕ W3 ( 0 ,

(67)

• Boundary between Region IV and Region V: for states ur ∈ S2 (urest ) a solution can be constructed consisting of the 2-shock connecting ur to urest followed by the 3-stationary wave connecting urest to the vacuum state (0, 0, al ). The structure of this solution is W3 ((0, 0, al ), urest ) ⊕ S2 (urest , ur ). • Boundary between Region IV and Region VI: for states ur ∈ SWsup (C − , ar ) that are above the curve W2 (urest ) a solution can be constructed like in Region VI, i.e (63). • Boundary between Region V and Region VI: For states ur ∈ SWsup (C − , ar ) that are below the curve W2 (urest ) two solutions like in Region VIa can be constructed, one of the form (63) and another one of the form (58).

3.3

Summary

Taking into account the order of the eigenvalues at every region, the sign of the velocities of the waves, and Theorems 2 and 3, it can be checked that the solutions constructed above are the only possible self-similar solutions consisting of simple waves linking some intermediate states. We summarize the results in this table: 22

Figure 14: Solution of a state in Region II. Left: projection of the intermediate states and the simple waves on the (h, u)-plane. Right: sketch of the free surface at a time t > 0. The direction of the movement of the simple waves is indicated with a red arrow.

Regions Region I Region II Region III

Number of solutions 1 1 1

Regions Region I Region II Region III Region IV Region V Region VIa Region VIb

Number of solutions 1 0 1 0 1 2 1

Case 1: initial condition (34) Form of the solutions W3 (ul , u+ u+ ur ) 0 ) ⊕ R1 ( 0 , − − + R1 (ul , u0 ) ⊕ W3 (u0 , u0 ) ⊕ R1 (u+ ur ) 0 , R1 (ul , ur ) Case 2: initial condition (35) Form of the solutions R2 (ul , ur ) R2 (ur0 , ur ) R2 (ul , u− u− u+ u+ ur ) 0 ) ⊕ W3 ( 0 , 0 ) ⊕ R2 ( 0 , s S2 (u0 , ur ) R2 (u0 , u− u− u+ u+ ur ) 0 ) ⊕ W3 ( 0 , 0 ) ⊕ S2 ( 0 , − − − R2 (ul , u0 ) ⊕ W3 (u0 , ur ) and R2 (u0 , u0 ) ⊕ W3 (u− u+ u+ ur ) 0 , 0 ) ⊕ S2 ( 0 , − − R2 (ul , u0 ) ⊕ W3 (u0 , ur )

In Regions where the number of solutions is 0, the structure of the solution corresponds to that of the considered partial Riemann problem for the homogeneous shallow water system.

4

Numerical tests

In Section 2 we have seen that in some cases the wet-dry Riemann problem may have zero, one, or two solutions. The goal of this section is to study how some Godunov-type methods behave in these different situations. We consider the following numerical fluxes for the homogeneous shallow water system: • Roe’s flux (ROE): see [15];

23

Figure 15: Solution of a state in Region IV. Left: projection of the intermediate states and the simple waves on the (h, u)-plane. Right: sketch of the free surface at a time t > 0. The direction of the movement of the simple waves is indicated with a red arrow. • Godunov’s flux (GODUNOV): see [21]; and the following numerical treatment of the source term: • Upwind discretization (UPW): see [4]; • Hydrostatic reconstruction (HR): see [2]; • Generalized hydrostatic reconstruction (GHR): see [10]. For the sake of shortness, these numerical techniques are not described here: the interested readers are addressed to the cited references. Let us only mention that, while the three numerical treatments of the source terms lead to methods that preserve water-at-rest stationary solutions (i.e. they satisfy the C-property, according to [4], only GHR leads to schemes that preserve any stationary solution and, in particular, stationary contact discontinuities. The following combination of these techniques will be used: • ROE UPW. • ROE HR. • ROE GHR. • GOD GHR. The numerical treatment of wet-dry fronts proposed in [6] is used in ROE UPW. In all the tests below, uniform meshes of the interval [-2,2], CFL = 0.9, and free boundary conditions are considered.

24

4.1

Tests 1 and 2

The goal of these tests is to study the behavior of the numerical methods in cases where the Riemann problem has no solution, i.e. when the step acts like an obstacle for the fluid. We consider first the initial condition:  ul = (0, 0, 1), if x < 0, u1 (x, 0) = (68) ur = (0.2, −2, 2), if x ≥ 0. It can be easily checked that the problem is of the form (35) and that the state ur is in Region IV: a solution of the form (66) was proposed in this case. In Figure 16 the numerical solutions obtained with a mesh of 400 cells are shown: all of them seem to converge to the proposed solution.

Figure 16: Numerical results of the h component for the initial condition (68). Secondly, we consider the initial condition:  ul = (0, 0, 1), if x < 0, u2 (x, 0) = ur = (2, 5, 2), if x ≥ 0.

(69)

It can be easily checked that in this case the state ur is in Region II: a solution of the form (65) was proposed in this case. In Figure 17 the numerical solutions obtained with a mesh of 400 cells are shown while in Figure 18 we use a mesh of 800 cells: all of them seem to converge again to the proposed solution.

4.2

Test 3

In this test we check the ability of the numerical methods to correctly capture stationary contact discontinuities. We consider the next initial condition:

25

Figure 17: Numerical results with 400 cells of the h component for the initial condition (68).  u3 (x, 0) =

ul = (0, 0, 1), if x < 0, ur = (5, 1, 2), if x ≥ 0.

(70)

Since the right state is in Region III, the exact solution is of the form (42). The numerical results obtained in a mesh of 400 cells are shown in Figure 19 and in Figure 20 we make a zoom on the stationary contact discontinuity. The numerical methods based on UPW or HR, as expected, are not able to correctly capture the stationary contact discontinuities: among them, ROE HR gives the numerical solution farthest to the proposed one. Although GHR leads to schemes that preserve these discontinuities, only GOD GHR seems to converge to the proposed weak solution: ROE GHR captures a contact discontinuity that does not satisfy the (MC) criterion (there is a transition from sub to supercritical through the stationary contact discontinuity).

4.3

Test 4

In [18] a test was shown where the numerical solutions obtained with ROE HR were very far of those produced by other methods. In fact, this initial condition was in Region VIa where two possible solutions have been found: one of the form (58) and the other of the form (63): let us check that ROE HR converges to the former while the other numerical methods converge to the latter. To do this, we consider the initial condition:  ul = (0, 0, 1), if x < 0, (71) u4 (x, 0) = ur = (1, −8, 2), if x ≥ 0. It can be checked that the right state is in Region VIa. Figure 21 shows the result obtained with a mesh of 400 cells. We observe that the only method that converges to the solution of the form (58) is ROE HR, while all the others schemes converge to the solution of the form (63).

26

Figure 18: Numerical results with 800 cells of the h component for the initial condition (68).

5

Conclusions

In this paper the Riemann problem for the shallow water equations corresponding to a wet-dry situation over a step has been studied. Depending on the state at the wet side, zero, one, or two solutions have been found. In the case of non-existence of solution, the problem has been interpreted as a partial Riemann problem for the homogeneous shallow water and a solution has been proposed. Since the problem is nonconservative, all the difficulties concerning the numerical approximation of weak solutions are found. To check this, we have compared the numerical solutions obtained with several standard numerical fluxes and treatment of the source terms. The main conclusions are the following: • In cases where the Riemann problem does not have solutions, the numerical methods converge to the proposed solutions. • In order to capture correctly the stationary contact discontinuities it is necessary to have a numerical method that preserve them. Nevertheless, this is not sufficient to ensure the convergence to the proposed solutions: the numerical solution may converge to weak solutions that do not satisfy the (MC) criterion. • The combination of Godunov’s numerical flux and the GHR technique seems to produce a numerical method that correctly captures all the proposed solutions. • In cases where the Riemann problem has two solutions, the numerical methods may converge to one or to the other.

27

Figure 19: Numerical results of the h component for the initial condition (70).

Figure 20: Zoom of the numerical results of the h component for the initial condition (70).

28

Figure 21: Numerical results of the h component for the initial condition (71).

29

References [1] F. Alcrudo and F. Benkhaldoun. Exact solutions to the riemann problem of the shallow water equations with a bottom step. Computers & Fluids, 30(6):643 – 671, 2001. [2] E. Audusse, F. Bouchut, M.O. Bristeau, R. Klein, and B. Perthame. A fast and stable wellbalanced scheme with hydrostatic reconstruction for shallow water flows. SIAM Journal on Scientific Computing, 25(6):2050–2065, 2004. [3] E. Audusse, C. Chalons, P. Ung. A simple well-balanced and positive numerical scheme for the shallow-water system. Communications in Mathematical Sciences, 13: 1317–1332, 2015. [4] A. Berm´ udez and M.L. V´aquez Cend´ on. Upwind methods for hyperbolic conservation laws with source terms, Computers & Fluids, 23: 1049–1071, 1994. [5] R. Bernetti, V.A. Titarev, and E.F. Toro. Exact solution of the riemann problem for the shallow water equations with discontinuous bottom geometry. Journal of Computational Physics, 227(6):3212 – 3243, 2008. [6] M.J. Castro, A.M. Ferreiro, J.A. Garc´ıa-Rodr´ıguez, J.M., Gonz´alez-Vida, J. Mac´ıas, C. Par´es, and M.E. V´azquez-Cend´on. The numerical treatment of wet/dry fronts in shallow flows: application to one-layer and two-layer systems, 42: 419 – 439, 2005. [7] M.J. Castro, J.M. Gonzalez-Vida, and C. Par´es. Numerical treatment of wet/dry fronts in shallow flows with a modified roe scheme. Mathematical Models and Methods in Applied Sciences, 16(06):897–931, 2006. [8] M.J. Castro, J.A. L´opez-Garc´ıa, and C. Par´es. High order exactly well-balanced numerical methods for shallow water systems. Journal of Computational Physics, 246:242–264, 2013. [9] M.J. Castro, T Morales de Luna, and C. Par´es . Well-balanced schemes and path-conservative numerical methods. In Handbook of Numerical Analysis, Vol. 18, pp. 131-175. Elsevier, 2017. [10] M.J. Castro, A. Pardo, C. Par´es, Well-balanced numerical schemes based on a Generalized Hydrostatic Reconstruction technique, Mathematical Models and Methods in Applied Sciences, 17: 2055–2113, 2007. [11] G. Chen and S. Noelle. A New Hydrostatic Reconstruction Scheme Based on Subcell Reconstructions. SIAM Journal on Numerical Analysis, 55: 758–784, 2017. [12] G. Dal Maso, Ph. LeFloch, and F. Murat. Definition and weak stability of nonconservative products. Journal de math´ematiques pures et appliqu´ees, 74:483–548, 1995. [13] O. Delestre, S. Cordier, F. Darboux, F. James. A limitation of the hydrostatic reconstruction technique for Shallow Water equations. Comptes Rendus Mathematique, 350: 677–681, 2012. [14] F. Dubois. 3.1 partial riemann problem, boundary conditions, and gas dynamics. Absorbing Boundaries and Layers, Domain Decomposition Methods: Applications to Large Scale Computers, page 16, 2001. [15] P. Glaister. Approximate Riemann solutions of the shallow water equations. Journal of Hydraulic Research, 26: 293–306, 1988.

30

[16] P.G. LeFloch and M.D. Thanh. The Riemann problem for the shallow water equations with discontinuous topography. Communications in Mathematical Sciences, 5(4):865–885, 2007. [17] P.G. LeFloch and M.D. Thanh. A godunov-type method for the shallow water equations with discontinuous topography in the resonant regime. Journal of Computational Physics, 230(20):7631–7660, 2011. [18] T. Morales de Luna, M.J. Castro, and C. Par´es. Reliability of first order numerical schemes for solving shallow water system over abrupt topography. Applied Mathematics and Computation, 219(17):9012–9032, 2013. [19] C. Par´es and M.J. Castro. On the well-balance property of roe’s method for nonconservative hyperbolic systems. applications to shallow-water systems. ESAIM: mathematical modelling and numerical analysis, 38(5):821–852, 2004. [20] G. Rosatti and L. Begnudelli. The riemann problem for the one-dimensional, free-surface shallow water equations with a bed step: Theoretical analysis and numerical simulations. Journal of Computational Physics, 229(3):760 – 787, 2010. [21] E. F Toro. Shock-capturing methods for free-surface shallow flows. Wiley and Sons Ltd., 2001. [22] X. Xia, Q. Liang, X. Ming, and J. Hou. An efficient and stable hydrodynamic model with novel source term discretization schemes for overland flow and flood simulations. Water Resources Research, 53: 3730–3759, 2017.

31