Riemann solution for a class of morphodynamic shallow water dam-break problems

This paper investigates a family of dam-break problems over an erodible bed. The hydrodynamics is described by the shallow water equations, and the bed change by a sediment-conservation equation, coupled to the hydrodynamics by a sediment transport (bed-load) law. When the initial states $\boldsymbol{U}_{l}$ and $\boldsymbol{U}_{r}$ are sufficiently close to each other the resulting solutions are consistent with the theory proposed by Lax (Hyperbolic Systems of Conservation Laws and the Mathematical Theory of Shock Waves, 1973, SIAM), that for a Riemann problem of $n$ equations there are $n$ waves associated with the $n$ characteristic families. However, for wet–dry dam-break problems over a mobile bed, there are three governing equations, but only two waves. One wave vanishes because of the presence of the dry bed. When initial left and right bed levels ( $B_{l}$ and $B_{r}$ ) are far apart, it is shown that a semi-characteristic shock may occur, which happens because, unlike in shallow water flow on a fixed bed, the flux function is non-convex. In these circumstances it is shown that it is necessary to reconsider the usual shock conditions. Instead, we propose an implied internal shock structure the concept of which originates from the fact that the stationary shock over a fixed-bed discontinuity can be regarded as a limiting case of flow over a sloping fixed bed. The Needham & Hey (Phil. Trans. R. Soc. Lond. A, vol. 334, 1991, pp. 25–53) approximation for the ambiguous integral term $\int \!h\,\text{d}B$ in the shock condition is improved based on this internal shock structure, such that mathematically valid solutions that incorporate a morphodynamic semi-characteristic shock are arrived at.


Introduction
A Riemann problem consists of an initial value problem composed of a set of conservation equations together with initial piecewise constant data having a single discontinuity.
In nonlinear shallow water flows, piecewise continuous solutions frequently develop. This is because the equations commonly used for describing them admit shocks (discontinuities) as solutions. These are usually interpreted as breaking waves (or bores), and therefore possess a straightforward physical significance, as well as a mathematical structure. These shocks are weak solutions in the sense that they satisfy the integral form of the flow equations. Smooth (differentiable) flow regions may be matched across these shocks by shock conditions, which can be derived by considering mass and momentum conservation across the shock.
Therefore, Riemann problems commonly occur in shallow water flows. Indeed, if one interprets all data in a shallow water numerical model as being piecewise continuous, then a whole series of such problems is solved at each time, and this interpretation forms the basis of a class of numerical shallow water solvers (Toro 2001).
From a physical standpoint, Riemann solutions in shallow water flows are important because they provide us with solutions to idealised problems that can be used as verification cases for numerical solvers. Additionally, these idealised problems serve to highlight fundamental shallow water dynamics. In shallow water flows a variety of Riemann problems are of interest. Some of the simplest of these are dam-break problems, which comprise a Riemann problem with zero initial velocities.
The simplest dam-break problem involves one wet and one dry side (wet-dry), and with the point of discontinuity corresponding to the position of a notional dam wall that at the initial time is instantaneously removed. Solution to this simplest shallow water dam-break problem is given in Stoker (1957). Although one could insist that a dam-break problem only be wet-dry, here we relax this description so as to include so-called wet-wet dam-break problems. These more generalised dam-break problems have a richer structure (see Toro 2001). We further consider dam-break problems in which the initial bed levels also are different (see Bernetti, Titarev & Toro 2008).
The Riemann problems in Stoker (1957), Toro (2001), Bernetti et al. (2008) are those with a fixed bed. If we allow the bed to become erodible, coupling flow velocity to movement of sediment via a bed-load sediment transport relation, then a more complex picture emerges. The solution of the shallow water mobile-bed wet-dry dambreak problem with no bed discontinuity at the dam location dates back to Fraccarollo & Capart (2002), who considered a system with separate layers for fluid and sediment. The equivalent problem without separate layers was considered by Kelly & Dodd (2009), amongst others. These problems are also important from a physical perspective because in real dam-break events considerable scour may result due to the high flow velocities.
In this paper we go further and consider generalised mobile-bed dam-break problems, in which initial bed levels are not, in general, equal across the initial data: see figure 1. Aside from being important in the context of true dam-break events, they also have relevance in the dynamics of waves on a beach. This is because a so-called backwash bore (i.e. a hydraulic jump) is frequently created when water runs back down the beach after a single wave uprush (Hibberd & Peregrine 1979). As the water drains the conditions at the hydraulic jump may be such that flow is minimal and thus subsequent development is predictable as a solution to a mobile-bed dam-break problem (Zhu & Dodd 2015).
Furthermore, we reconsider the usual shock conditions. For non-constant bed levels the bed-slope term in the flux-conservative form of the momentum equation is not integrable, necessitating approximation using conditions on each side of the shock. Here we reinterpret this term based on a new approximation of the internal shock structure. Although the motivation for this comes from the limiting case of flow over a sloping bed (and against which we subsequently verify the method), we instead approximate the internal morphodynamic shock structure as a series of sub-shock problems.
Dam FIGURE 1. Schematic diagram for a dam-break problem.
In the next section we present our governing equations, as well as some of the theory governing the determination of the wave structure across the evolving Riemann solution. We then use this theory to solve this class of Riemann problem in § 3. Finally we present conclusions.

Governing equations
A one-dimensional (1-D) idealised configuration for the initial set-up of a generalised dam-break problem is shown in figure 1. As mentioned, the nonlinear shallow water equations (NSWEs) have often been used for describing one-or two-dimensional dambreak flows (Ritter 1892;Stoker 1957;Toro 2001;Bernetti et al. 2008). For dambreak problems over a mobile bed, and if only bed load is considered (see Soulsby 1997), the governing equations are the NSWEs and a sediment-conservation equation: ht +ûĥx +ĥûx = 0, (2.1) wherex represents horizontal distance (m),t is time (s),ĥ represents water depth (m),û is a depth-averaged horizontal velocity (m s −1 ),B is the bed level (m),q is sediment flux due to bed load (m 2 s −1 ), ξ = 1/(1 − p) with p being bed porosity and g is acceleration due to gravity (m s −2 ). In order to reveal the shock dynamics by solving a strictly hyperbolic system, we do not include the downslope diffusion effect in our model, although morphodynamic shocks are considered where vertical bed steps occur.
In general,q is strongly dependent onû and a weak function ofĥ. Here, a simple but commonly used formulaq = Aû 3 (see Grass 1981) is employed for the bed load (see e.g. Kelly & Dodd 2010;Zhu, Dodd & Briganti 2012), with A being the bed mobility parameter (s 2 m −1 ). Note that this formulation is an over-simplification of the complex process of bed-load transport (see e.g. Pritchard & Hogg 2005) but that the purpose here is to construct the mathematical solution so that the basic dynamics can be understood, and to provide a mathematical test case for numerical models.
Therefore, (2.3) becomesBˆt whereĥ 0 is a length scale, andû 0 = (gĥ 0 ) 1/2 . Substituting (2.5) into the governing equations (2.1), (2.2) and (2.4) gives The vector form of these three non-dimensional governing equations is The eigenvalues of A are the roots of the polynomial equation (2.11) Equation (2.11) has three roots, denoted λ 1 , λ 2 and λ 3 such that λ 1 λ 3 λ 2 . For the solution of λ 1 , λ 2 and λ 3 we refer to Kelly & Dodd (2009, 2010. 2.3. Generalised simple wave theory When there are two equations or fewer in a Riemann problem, we can derive a Riemann invariant along each characteristic (Stoker 1957). Jeffrey (1976) shows, however, that a direct extension of the concept of a Riemann invariant is not possible when there are more than two equations and dependent variables in a Riemann problem. Therefore, the generalised simple wave theory and the generalised Riemann invariants are introduced to solve such Riemann problems.
If we consider, for the moment, a general quasilinear hyperbolic system where U is a vector of n dependent variables, given by then the assumption that a simple wave region exists in a Riemann problem (or indeed generally) is such that U = U(u 1 ) holds across a wave, where u 1 is chosen without loss of generality. This means that there is a functional dependence between u i and u 1 of the form u i = f i (u 1 ), and the wave is called a generalised simple wave (Jeffrey 1976). For a simple wave, equation (2.12) can be written as (2.14) This system can have a non-trivial solution for dU/du 1 only if This implies that u 1 , and therefore u i , are constant along characteristics (dx/dt) = λ i , which are themselves straight lines. A simple wave can also be defined as the wave, across which one family of characteristics are all straight lines. There are n families of characteristics, and the simple wave associated with the ith characteristic family is called the λ i simple wave. From (2.14), when λ = λ i the vector dU/du 1 must be proportional to the right eigenvector R (i) of A, which gives (Jeffrey 1976) where r (i) j is the jth component of R (i) . Integrating (2.17) yields the jth generalised Riemann invariant K j associated with the λ i simple wave: (2.18) Further details of the simple wave theory can be found in Jeffrey (1976).
Returning to the present system we have For a fixed t, where h varies continuously, equations (2.20) and (2.21) may be integrated across the λ i simple wave, in the (h, u, B) phase space, to yield the structure. However, variables do not always vary continuously across the Riemann structure. Therefore, we must understand the structure before we can solve it. 2.4. Wave structure for a Riemann problem In general the Riemann problem associated with (2.12), is composed of n waves associated with the n characteristic families, which are separated by n − 1 newly formed constant regions provided that the values in the initial piecewise constant states are sufficiently close (Lax 1973;Fraccarollo & Capart 2002;Toro 2009): see figure 2. Note that a wave could be a rarefaction fan, a shock or a contact wave.
As we integrate (2.20) and (2.21) across a λ i simple wave, λ i varies such that dλ i = (∂λ i /∂u j ) du j , and in view of (2.17): where we have assumed, without loss of generality, that r (i) 1 = 1. Equation (2.22) therefore represents the rate of change of characteristic velocity across the λ i wave. If ∇ U λ i · R (i) = 0 for all U, (2.23) then the characteristic velocity increases or decreases continuously across the λ i wave, thus implying the existence of either an expansion (rarefaction) wave, or a compressive wave (which becomes a shock ultimately). Then the λ i wave field is said to be genuinely nonlinear (or convex) (Sharma 2010). Note that if the wave field is said to be genuinely linear. If we return again to the present system we find that , and it is therefore not clear that (2.25) will be strictly > or < 0. In this case the problem is said to be non-convex (Sharma 2010). Note that this is in contrast to the classical (fixed-bed) shallow water system, wherein This implies that the solution to the dam-break problem resulting from the system (2.6)-(2.8) may possess semi-characteristic shocks (Sharma 2010).
. Schematic diagrams depicting possible structures for the combination of a rarefaction wave and a semi-characteristic shock within a simple wave.

Semi-characteristic shock
In general, ∇ U λ i · R (i) = f (u 1 ) across a simple wave. If, as seems possible for (2.25), f (u 1 ) passes through 0, this implies that at a rarefaction fan, across which characteristics must diverge, a point is reached at which this can no longer happen because (dλ i /du 1 ) = f (u 1 ) = 0. The only way then of accommodating this behaviour without the solution becoming multivalued is for a shock to be contiguous with the fan: the semi-characteristic shock, with the characteristic at one edge of the fan coinciding with this semi-characteristic shock. The possible three wave structures are shown in figure 3. Note that only (a) represents a physical structure in general. This is because for figure 3(b,c) the structures will, in general, be overdetermined. For further details and exceptions see Sharma (2010).

Shock conditions
We require shock conditions to be satisfied across shocks and semi-characteristic shocks. For derivations of the shock conditions we refer to Kelly & Dodd (2010), Zhu et al. (2012); the shock conditions are where L and R represent variables on the left and right side of a shock, and W is shock velocity.
The term is not well defined along the face of bed step (Kelly 2009). Needham & Hey (1991) performed the integration by approximation: This expression is adopted by Kelly & Dodd (2010), Zhu et al. (2012), Zhu & Dodd (2013. Bernetti et al. (2008) equated it to the hydrostatic pressure force exerted on the bed step by the water: (2.31) There are four cases depending on from which side the water is exerting a force on the bed step, and whether the free surface on the side of lower bed level is above the top of the bed step. Note that when free surface elevations on both sides are equal, equations (2.30) and (2.31) are identical.
Different interpretations are possible. This can be seen by assuming that in the vicinity of a sudden change in bed level (a bed step) is the Dirac delta function, and x = x 0 is the location of the bed step. Using (2.32) Therefore, the remaining ambiguity is in how to define h(x 0 ). From this perspective the expression of Needham & Hey (1991) (hereinafter NH91) is a simple average depth across the discontinuity; similarly, the middle two of the expressions of Bernetti et al. (2008) also correspond to composite depths, although the first and last cannot obviously be interpreted in this way. More generally, there is an implied variation of depth across the shock, which, it turns out, is important for the shocks we consider here.
Subcritical Supercritical FIGURE 4. Bed level geometry for this case. F =û/ gĥ is the Froude number. In these scenarios the whole flow is either sub-or supercritical.

Investigation of
x R x L hB x dx If we know the internal structure of a shock, i.e. h = h(B) across the shock, we can calculate the ambiguous integral numerically or analytically straightforwardly: Needham & Hey (1991) approximation, i.e. (2.30). We refer to this approach as the n-step approach; it is based on an implied internal structure of a shock. To see how this internal structure might be arrived at from a physical standpoint we now consider flow over a fixed-bed slope.

Interpretation of a stationary shock across a fixed-bed discontinuity
We consider a steady state flow across a fixed, linear slope (see figure 4), and therefore only the hydrodynamic behaviour. From shallow water theory, the flow over the slope is either entirely sub-or supercritical (see appendix D). We return here to dimensional variables, and then use an alternative non-dimensionalisation.
From continuity we haveĥ For a steady state, the flux-conservative form of (2.2) is (2.36) Now we introduce a different set of non-dimensional variablesh,ũ,x andB on the sloping section withĥ =ĥ Lh ,û =û Lũ ,x =ĥ Lx / tan α andB =ĥ LB . This gives Note that the slope tan α is now absent, and the only free parameter is the inflow Froude number, F L . Straightforwardly, we then obtaiñ for the variation ofh across the slope. If we consider an abrupt change in bed level to be the limiting case as α → π/2 of this linear slope variation, and, moreover, that this variation is independent of slope (tan α), we may then assume that this variation may be used across a fixed-bed step as the implied internal shock structure.
It should be noted that (2.40) can also be directly derived from an energyconservation law for the shallow water equations, because the flow down the slope is continuous. There is some debate about whether energy conservation or the momentum balance equation, in which B R B Lh dB has to be approximated, should be used for the shock across a fixed-bed step (Cozzolino et al. 2011;Valiani & Caleffi 2017). The energy-conserving approach has been widely used to study the stationary shock across a fixed-bed discontinuity (Karelskii & Petrosyan 2006;Valiani & Caleffi 2017). Valiani & Caleffi (2017) uses both energy conservation and momentum balance to derive a depth at the bed step, i.e. h(x 0 ) in (2.33). However, the energy loss across a morphodynamic shock is a priori unknown. Therefore in this work, we utilise the momentum equation to solve the stationary shock across a fixed-bed step and also for morphodynamic shocks. Accordingly, we now focus on the approximations for for a stationary shock across a fixed-bed discontinuity, in whichh is calculated using (2.40). The exact solution (2.41) together with (2.40) allows us to see how well (2.30), (2.31) or (2.34) describe this usually ambiguous integral. The performances of these approximations are presented in appendix A, from which we can see that (2.34) yields greater accuracy than (2.30) or (2.31).

Application of n-step approach to morphodynamic shocks
In this section we consider whether the n-step approach is valid for morphodynamic shocks. Hereafter we return to the non-dimensionalisation introduced in § 2.2.
In the shallow water morphodynamical system that we consider here, there is in general one characteristic speed much smaller than the other two. This can be seen Note that the characteristic polynomial for λ depends only on F and σ : The characteristic speed that is generally much smaller than the other two is associated with bed wave movement. This property pertains everywhere except for transcritical flows, as also indicated in figure 5. Here we define a morphodynamic shock as a shock formed by the convergence of two characteristics of one family, at least one of which is a bed characteristic, and for which the shock speed |W| 1. Now, note that σ = ξ Ag 1 (see (2.8) and (2.29)). This is becauseq = O(10 −3 ) m 3 s −1 m −1 or less, whereasû 0 = O(1) m s −1 in our original nondimensionalisation (2.5). This implies that at the hydrodynamical time scalê t 0 = ĥ 0 /g, equation (2.8) becomes B t ≈ 0, implying no bed change at this time scale. However, at the morphodynamical time scale,t m =t 0 /σ , ⇒ both sides of (2.29) are of comparable magnitude. This is consistent with W 1 for a morphodynamic shock, so that in these circumstances (2.27)-(2.29) become  .27) and (2.28) to obtain the internal shock structures (h i , u i and B i ) and apply the n-step approach (2.34) for (2.27)-(2.29). In the remainder of this paper we use this approach, which we refer to as the n-step approach, to construct the Riemann solution. The performances of these other approximations for morphodynamic shocks are examined in appendix C. 2.7.3. Implementation of n-step approach for morphodynamic shocks To use (2.34) across a morphodynamic shock with (2.27)-(2.29), the procedure is as follows: (i) Obtain initial estimates for U R = U (1) R by solving (2.27)-(2.29), with x R x L hB x dx approximated using (2.30). x L hB x dx using the n-step approach (2.34). (iv) We then solve (2.27)-(2.29) using the calculated

Solution of dam-break problems over an initially piecewise flat mobile bed
The dam-break problem system consists of three equations, and according to Lax's theorem (Lax 1973) there are at most four constant states separated by three elementary waves associated with the three characteristic families. Note that for wet-dry dam-break problems over a mobile bed there are two waves separated by one newly formed constant region (Kelly & Dodd 2009). One wave vanishes because of the presence of the dry bed.

Initial conditions
The dimensional initial conditions for a generalised dam-break problem are shown in figure 1. Withĥ 0 =ĥ l , the non-dimensional initial conditions of the left side are h(x 0) = h l = 1, u(x 0) = u l = 0, B(x 0) = B l = 0. For wet-dry dam-break problems, h(x 0) = h r = 0, and u(x 0) = u r = 0. For wet-wet dam-break problems, we set h(x 0) = h r = 0.1, and u(x 0) = u r = 0. In this paper, we consider conditions of both B r = 0 and B r = 0 to investigate the wave structures in these more generalised dam-break problems.

Wet-dry dam-break problem
We first assume that the wet-dry dam-break solution over an erodible bed for various B r consists only of elementary waves, i.e. rarefaction waves or shocks. We introduce the semi-characteristic shock when the assumption no longer applies. The obtained wave structures at t = 1 are shown in figure 6.

B r B l
Figures 6(a) and 6(b) show, respectively, water surface elevation and bed level, and velocity. The wave structure is that of a λ 1 rarefaction wave, a constant region, U * , and a λ 3 rarefaction wave, which is consistent with that presented by Kelly & Dodd (2009), who considered only B r = B l = 0. As B r → B l + h l , B * + h * → B l + h l and the extents of the λ 1 and λ 3 rarefaction waves decrease, and the solution (at t = 1) resembles the initial conditions more. Note that the volume of water set in motion at time t is h l |λ 1 (h l )|t, and is independent of B r because the left edge of the λ 1 fan (λ 1 (h l )) is unaffected by changes in B r . As the downstream elevation B r increases, velocities across the Riemann solution decrease, as, therefore, does the sediment movement. As B r increases, the flow in the constant region changes from supercritical (e.g. when B r = 0) to subcritical flow (e.g. when B r = 0.4 or 0.8), and the λ 3 wave close to x = 0 changes from a hydrodynamic into a bed wave. In the λ 3 rarefaction fan, where (2.11) has been used. Thus, the large bed change that occurs near the dam location for B r = 0.4 or 0.8 is connected by the λ 3 simple wave with λ 3 → 0 in (3.1). The lip of the initial bed discontinuity is eroded by the flow, and the initial discontinuity in bed level is transformed into a steep continuous variation. The small bed step (Kelly & Dodd 2010) at the flow tip (x = x s (t)) remains a feature of the solutions with decreasing height as flow velocity there decreases as B r increases. For B r = h l + B l no flow ensues, and there is no erosion. The wave development is closely related to Froude number (F), and Froude number acts as a proxy for position x. In figure 7( f ), we show the relationship between F h are plotted as a function of Froude number (F). Again, note that the λ curves are invariant for all dam-break solutions with σ = 0.01 unless there is a discontinuity in F when a shock develops. In contrast, the black line superimposed on parts of these curves indicates the variation of λ 1 across the λ 1 fan, and the variation of λ 3 across the λ 3 fan, with the jump from one to the other also depicted, for B r = 0 and 0.4. If we follow the (black) λ i values along the λ i curves in figure 7(a,b) from F = 0 we see that λ 1 increases as F increases for B r = 0.4 and 0. For B r = 0 the jump from λ 1 to λ 3 for B r 0, which corresponds to the constant region, occurs for a larger F value than that for B r = 0.4. Both jumps occur prior to the point at which dλ 1 /dF = 0. Therefore, λ increases monotonically across both fans. (dλ 1 /dF) = 0, indicating a convergence in λ 1 , occurs at F ≈ 1.6 for all dam-break solutions with σ = 0.01.

B r B l
In figure 6(c,d) we can see the dam-break structure in this case, which is similar to the preceding one in that two rarefaction fans (λ 1 and λ 3 ) form, separated by a constant region.
Before commenting further on the structure of these solutions it is instructive first to consider their representation in (λ , F) space. In figure 7(c-e) we do this. Figure 7(c,d) illustrates the behaviour for B r = −0.4. For figure 7(c) we see the elementary wave solution. Note, however, that the jump via the constant region, from λ 1 to λ 3 curve, occurs when (dλ 1 /dF) < 0. This implies that at some point within the λ 1 fan the λ 1 derived characteristics start decreasing as F increases. In contrast figure 7(d) shows the behaviour with a semi-characteristic shock included.
Representing the solution in (λ , F) space is appealing because these curves are invariant with the continuous Riemann solution (i.e. size of bed step). However, it is the convergence of λ 1 characteristics (not λ 1 ) in the (x, t)-plane that determines whether or not a semi-characteristic shock must be fitted. So, to determine this point of change in the Riemann solution it is appropriate to examine variation in λ 1 . The multivalued solution of λ 1 starts at B r −0.175, at the location at which F ≈ 1.87 (and h ≈ 0.301). This point is illustrated in figure 8. λ 1 increases as F increases, but at F ≈ 1.87 the characteristic velocity (λ 1 ) starts to decrease at the leading edge of the λ 1 rarefaction fan. The constant region, corresponding to the jump from λ 1 to λ 3 curve, thus occurs when (dλ 1 /dF) < 0. This indicates the convergence of λ 1 characteristics within the λ 1 fan. The solution is multivalued. The part of the λ 1 curve for which F > 1 behaves like a characteristic associated with a bed wave. Because λ 1 < 0 this behaviour results in a bed wave propagating against the flow. Similar behaviour can be observed in the propagation of anti-dunes in supercritical open channel flow on a mobile bed, which propagate against the flow (see Kennedy 1963).
To obtain a valid mathematical structure here, the λ 1 fan must terminate prior to the point at which (dλ 1 /dh) = 0 at a semi-characteristic shock with λ 1L = W; downstream, in the constant region we must have λ 1R < W, for a valid shock structure, which is possible because F increases across the shock and therefore λ 1 decreases.
Physically, the main difference between this case and that for B r > B l is the larger velocities induced by the lower downstream elevation, which drives the early supercritical flow development and therefore the formation of the semi-characteristic shock. Now, the constant region is shifted mostly to x > 0, with a large decrease in h * . The flow is supercritical at the dam location, and the λ 1 wave has become a bed wave close to x = 0, where we can also see the large bed change. The large bed decrease occurs at the shock, which helps to connect widely separated values of B l and B r .
We assume an implied internal shock structure for all the semi-characteristic shocks, and n = 2 is adopted here. For the semi-characteristic shocks for all negative B r values, we have λ 1L = W > λ 1R . The convergence of characteristics implies that the shocks are physical. Note that if we use the approximation (2.30) we do not arrive at physical solutions for B r = −0.8, and if we use the approximation (2.31) we do not get solutions for any semi-characteristic shocks with the examined B r values. It is therefore critical that we introduce this more accurate approach. In appendix C we show the equivalent solutions and examine dependence on n.

Varying upstream bed level
We also examine the effect of varying the upstream bed level only, while keeping the upstream surface elevation and downstream bed level fixed (h l + B l = 1 and B r = 0). The Riemann solutions (B and B + h only) are shown in figure 9. The structures are similar to those already observed, but now with particularly large variations between solutions for x < 0, because h l varies and so therefore does the driving force. The wave solutions of these dam-break problems are similar to those of fixed h l and B l values (h l = 1, B l = 0) but varying B r values. Actually, these two kinds of dambreak problems can be converted to each other through scaling and transformation. It is the water depth on the left and bed difference that determines the wave structure. Two dam-break problems with the same ratios of water depths and bed differences, i.e. h l1 /h l2 = (B l1 − B r1 )/(B l2 − B r2 ) are similar. Therefore the wave structures after dam collapse are similar.

Wet-wet dam-break problem
We now turn to the wet-wet dam-break solution over an initially piecewise flat erodible bed for various B r values, which consists of elementary waves. The solutions at t = 1 are shown in figure 10. 3.3.1. B r B l Figure 10(a,b) shows, respectively, water surface level and bed level, and velocity for this case. The wave structure is of a λ 1 rarefaction wave, a λ 3 rarefaction wave and a λ 2 shock. There are two newly formed constant regions (U l * and U r * ) separating the three waves. For B r = 0, the λ 2 shock corresponds to the λ + shock in the equivalent fixed-bed wet-wet dam-break problem (e.g. Toro 2001), and the λ 1 and λ 3 rarefaction waves correspond to the λ − rarefaction wave in the fixed-bed case. When the bed mobility σ → 0, the λ 1 and λ 3 rarefaction waves tend to combine into one wave.
As B r increases, the λ 3 wave becomes more confined to the original bed step position, and less water flows into x > 0 region. When h r + B r → 1 flow ceases.

B r B l
In figure 10(c,d), we can see the dam-break structure in this case. There is still a λ 1 rarefaction wave, and a λ 2 shock, as for B r > B l . However, for the investigated negative B r values, the λ 3 wave changes from a rarefaction wave into a shock. There are also two newly formed constant regions. For B r < 0 as B r decreases further (see figure 10c), h l * decreases very quickly, and when h l * < h r * the λ 3 fan observable for B r = 0 becomes a shock as the characteristics converge.
When B r −0.207, we get multivalued solutions within the λ 1 fan, (see figure 10c). The B r value at which this occurs will depend on h r , which here is 0.1. We again assume the existence of a λ 1 semi-characteristic shock, and the corresponding wave structure is shown in figure 10(c,d). Once again, we assume an implied internal shock structure, expressed through (2.34) in (2.27)-(2.29), in order to obtain a physical shock (see appendix C).

Varying downstream water depth
In figure 11 we look at the effect that varying h r has on the structure of these problems. As h r decreases, we expect this wet-wet problem to start to resemble previously examined wet-dry problems (see figure 6). Accordingly, the λ 3 shock diminishes such that between h r = 0.03 and 0.015 it becomes a rarefaction fan. As h r decreases further the λ 3 fan extends towards the λ 2 shock such that in the limit h r → 0 the leading edge of this λ 3 fan becomes the wet-dry boundary (with zero depth) and a semi-characteristic shock, and the λ 2 shock disappears. The λ 1 wave is a combination of a fan and a semi-characteristic shock, which is consistent with the equivalent wet-dry dam-break solution.

Conclusion
Generalised wet-dry and wet-wet dam-break problems over an erodible, initially piecewise flat bed with water initially at rest and discontinuous bed levels are investigated and solved based on the Riemann theory, and quasi-exact solutions are presented. The solutions are consistent with the theory proposed by Lax (1973) that for a Riemann problem of n equations there are at most n + 1 constant states separated by n elementary waves associated with the n characteristic families. However, for wet-dry dam-break problems, one wave vanishes because of the presence of the dry bed. For the examined wet-dry dam-break problems, there are two waves separated by one newly formed constant region, which is in agreement with that in Kelly & Dodd (2009).
For some dam-break problems with negative bed steps, in which the initial states (B l and B r ) are not sufficiently close, there are multivalued solutions when applying Lax's theorem. The semi-characteristic shock is introduced to describe the flow and physical wave structures are obtained. This is consistent with the solution for a Riemann problem of one single equation of non-convex flux function (Sharma 2010).
The ambiguous integral x R x L hB x dx in shock conditions, which is usually approximated by the Needham & Hey (1991) approach, is reconsidered. An implied internal shock structure is proposed initially by considering the limiting case of flow down a linear slope over fixed bed. Based on the internal shock structure, the integral x R x L hB x dx is discretised into many steps and each step is approximated by the Needham & Hey (1991) approach. This is to reduce the effects of curvature between h and B, which is source of inaccuracy in the Needham & Hey (1991) approximation. This strategy is then extended to morphodynamic shocks which are here by assumption slow moving. However, because of the more general implied internal shock structure approach we ultimately adopt it would appear that this approach is more generally applicable. The resulting semi-characteristic shocks are physical, in that the characteristics converge across them. Here, we takeB L = 0 andB R = −x, withx > 0 being a variable. The approximation (2.30) gives and it is exact if the relationship betweenh andB is linear. The approximation (2.31) is We also use the n-step approach (2.34) with n = 5 for the approximation withB i = −ix/5, in which h i is calculated by (2.40). The results of (A 1), (A 2), (A 3) and (2.41) are illustrated in figure 12. Comparison shows that the approximation (A 1), i.e. (2.30), is generally quite accurate. Nearer to critical conditions the free surface curvature introduces significant discrepancies. Also, asB R −B L increases, accuracy diminishes. The n-step approach (A 3), i.e. (2.34), greatly increases the accuracy. The approximation (A 2), i.e. (2.31) (Bernetti et al. 2008), is indistinguishable from the exact value when it is subcritical flow. It slightly overestimates the exact value when it is closer to critical flow. However, when the flow is supercritical, the discrepancies become significant.

Appendix B. Comparison between wet-dry dam-break Riemann solution of fixed and nearly fixed-bed case
In this section, we test the Riemann solver by comparing the nearly fixed-bed wetdry dam-break solution (σ = 1 × 10 −5 ) against the fixed-bed solution.
The dam-break problem over a fixed bed with a bed step will lead to a stationary shock at the bed step (Bernetti et al. 2008). For the positive B r values examined, the wave structure over a fixed bed from the left to right is a λ − rarefaction wave, a constant region, a stationary shock and a λ − rarefaction wave. That for the negative B r values examined is a λ − rarefaction wave, a stationary shock, a constant region and a λ − rarefaction wave ( figure 13). This stationary shock corresponds to a semicharacteristic shock if we interpret it in a morphodynamic context. It is a λ 3 semicharacteristic shock with λ 3L = W = λ 3R = 0 because σ = 0. On a nearly fixed bed, the wave structure for the non-negative B r values examined, is a λ 1 rarefaction wave, a constant region and a λ 3 rarefaction wave. The stationary shock in the fixed-bed case disappears because of the bed mobility, and the left part of the λ 3 rarefaction wave becomes a steep but smooth part, which is to some extent similar to the stationary shock in the fixed-bed case. The nearly fixed-bed solutions are in good agreement with the fixed-bed solutions (figure 13a).
The Riemann solutions with the NH91 condition for the examined negative B r values over a nearly fixed bed are similar to those presented in § 3.2.2. The wave structures are a λ 1 rarefaction, a λ 1 semi-characteristic shock and a λ 3 rarefaction. The λ 1 rarefaction fan corresponds to the stationary shock on the fixed bed. The results compare favourably with the corresponding fixed-bed results (figure 13b). In this section, we investigate the approximation of B R B L h dB by comparing the solutions of wet-dry and wet-wet dam-break problems using different approximation methods (see appendix A and § 2.7.2). In summary, the approximation methods include (a): (2.41) together with (2.40) and (2.29), i.e. h R and u R are directly calculated for a known B R as for the stationary shock across a fixed-bed discontinuity;  method (c) and (d) are close. When the bed step height increases, the difference grows. However, the semi-characteristic shocks predicted by method (a) and (b) for B r = −0.4 (and B r = −0.2, not shown) and that by method (d) for B r = −0.8 (and B r = −0.6, not shown) are non-physical because λ 1L W λ 1R cannot be satisfied. In contrast, method (c), used throughout this paper, results in a physical semi-characteristic shock for all these B r values. It is therefore necessary to retain the W terms in (2.27)-(2.29).
The effects of how many steps the integral B R B L h dB is discretised into are also investigated. The comparison for both wet-dry and wet-wet dam-break problems using method (c) (n-step approach) and (d) (NH91 approximation) is shown in figure 15. We can see that the results with the implied internal shock structure are initially very close to those directly using the NH91 approximation. As the size of bed step increases, the difference between solutions with the NH91 approximation and the new approach also increases slightly.
In the wet-dry dam-break problems with B r = −0.2, −0.4, the semi-characteristic shocks have λ 1L = W > λ 1R when solved directly with the NH91 approximation. However, as previously mentioned, for B r = −0.6, −0.8, we have λ 1L = W < λ 1R , which indicates a non-physical shock. It is further noted that in the B r = −0.4, −0.6, −0.8 cases, the water on the left and right sides of the λ 1 semi-characteristic shock are separated by the high bed step ( figure 16). Similarly, in the wet-wet dam-break problems with B r = −0.6, −0.8, the semi-characteristic shocks are also non-physical. However, when we use the implied internal shock structure, the shocks become physical.
In order to further investigate this, we take the wet-dry dam-break problem with B r = −0.6 as an example to illustrate the characteristics across the semi-characteristic shocks; see figure 17. We can see the characteristics diverging when using the NH91 approximation, and characteristics converging when using implied internal shock structure. This indicates that the NH91 approximation becomes less accurate when the bed step becomes large in which case the curvature of h with B becomes enhanced. In this case, the importance of considering the internal shock structure becomes obvious.
In figure 18 we see how the multivalued problem in a wet-dry dam-break problem is rationalised by introducing a semi-characteristic shock. According to Whitham (1974), the areas A 1 = A 2 . The results obtained here are consistent with this . Structure of the wave solution (NH91 approximation) for a wet-dry dambreak problem with B r = −0.6 which shows that the water on the two sides of the semicharacteristic shock is separated by the high bed step (t = 1). law, and when n increases, | A 1 − A 2 | decreases. This also demonstrates that the n-step approach applied for (2.27)-(2.29) gives more accuracy. Note that when the NH91 condition is used, the jump at the semi-characteristic shock occurs outside the multivalued region. As a result, the shock becomes non-physical, because λ 1 * > λ 1L = W.
Appendix D. Impossibility of a smooth flow from sub-to supercriticality down a slope For smooth transition from sub-to supercriticality down a slope we need a situation like that depicted in figure 19. This flow is described by (2.39). Note that the right-hand side of (2.39) > 0 for all x; therefore the left-hand side > 0 as well.
If we differentiate the left-hand side of (2.39) with respect to x we get  However, for the flow in figure 19,hx < 0. Therefore, we must have Now, Therefore, if flow is subcritical the flow cannot exist. Therefore, the flow in figure 19 cannot exist. The authors could not find an example in the literature of this analysis being presented, hence its inclusion here.