Numerical Simulation of Bed Load and Suspended Load Sediment Transport Using Well-Balanced Numerical Schemes

Noname manuscript No. (will be inserted by the editor) Numerical simulation of bed load and suspended load sediment transport using well-balanced numerical schemes Gonz´alez-Aguirre, J. C. ¨Gonz´alezV´azquez, J. A. ¨Alavez-Ram´ırez, J. ¨ Silva, R. ¨V´azquez-Cend´on, M. E. Received: date / Accepted: date Abstract Sediment transport can be modelled using hydrodynamic models based on shallow water equations coupled with the sediment concentration conservation equation and the bed conservation equation. The complete system of equations is made up of the energy balance law and the Exner equations. The numerical resolution for this complete system is done in a segregated manner. First, the hyperbolic part of the system of balance laws is resolved employing a finite volume scheme which is based on the Q-scheme of van Leer for computing the numerical flux also the numerical flux is computed by using an HLLCS approximate Riemann solver. As well, the hyperbolic system of balance laws is solve taken into account a non conservative form. The discretization of the source terms is carried out according to the numerical flux chosen. In the second stage, the bed conservation equation is resolved by using the approximation computed for the system of balance laws. The numerical schemes have been validate making comparison between the obtained numerical results and the experimental data for some ones physical experiments. The numerical results show in good agree with the experimental data. Divisi´on Ac´ademica de Ciencias B´asicas, Universidad Ju´arez Aut´onoma de Tabasco Carr. Cunduac´an-Jalpa de Mend´ez KM. 1, CP 86690, Cunduac´an, Tabasco, M´exico E-mail: [email protected] Universidad Aut´onoma Metropolitana-Unidad Iztapalapa Av. San Rafael Atlixco no.186 Ciudad de M´exico, Mexico E-mail: [email protected] Divisi´on Ac´ademica de Ciencias B´asicas, Universidad Ju´arez Aut´onoma de Tabasco Carr. Cunduac´an-Jalpa de Mend´ez KM. 1, CP 86690, Cunduac´an, Tabasco, Mexico E-mail: justino.ala[email protected] Instituto de Ingenier´ıa, Universidad Nacional Aut´onoma de M´exico 04510 Ciudad de M´exico, Mexico E-mail: rsilv[email protected] Departamento de Matem´atica Aplicada, Universidade de Santiago de Compostela 15706 Santiago de Compostela, Spain E-mail: elena.v[email protected] 2 Gonz´alez-Aguirre, J. C. et al. Keywords Sediment transport ¨Suspended load ¨Bedload ¨Finite volume method ¨Numerical simulation ¨Well-balanced schemes 1 Introduction The mathematical modelling of sediment transport has become more important subject as of extreme hydrological events intensify and increase in numbers, as a result of climate change, see [8]. According to [23], the mathematical models proposed can be classified as fully coupled models (FCM), partially coupled models (PCM) and decoupled models (DM). The main difference between them is the manner in which they deal with the interaction between water flux, sediment transport and bed evolution. The FCM take into account the hydrodynamic and both the suspended load and the bedload transport. In these models, the bed evolution is modelled using the Exner equation (see [16, 17]) with source term. The PCM consider the hydrodynamic and suspended load sediment transport. In this kind of model, the bedload transport is not taken into account for modelling the bed evolution. The DM compute the sediment transport using the shallow water equations coupled with the bed conservation equation. In this kind of model the suspended load sediment transport is neglected. Another difference between the models is the way as each one them consider the interaction between the flux, sediment transport and the morphological evolution inside the structure of mass conservation equation and momentum conservation equation. The FCM can have up to two additional terms on the right hand side of mass conservation equation. These terms quantify the rate of water entrainment (e.g. rainfall, infiltration, etc.) and the bed deformation. On the right hand side of momentum conservation equation two additional terms are added, which measure the effects of variable sediment concentration and momentum transfer due to sediment exchange between the water and movable bed ([8]). Examples of FCM can be found in [8,23,27,3]. These models consider hydrostatic pressure. Unlike them, in [7] a non-hydrostatic system is used, in this kind of approach, special focus must be put in the appearances of numerical instabilities which are given by the non-hydrostatic distribution. The PCM do not consider these additional terms in the shallow water equations because the sediment conservation and bed conservation are computed with the information provided by the hydrodynamics of clear water. Therefore, the interaction that may exist between the flux and sediment is ignored for the hydrodynamic equations. Examples of this kind of model can be found in [9] and [29]. In order to compute the morphological evolution, the majority of PCM do not take into account the bedload transport, although in works such as [36,29], it is included. The DM use the shallow water equations coupled with the bed conservation equation for modelling the sediment transport flux. The bed conservation, in this kind of model, presents a convective term that quantifies the bedload Title Suppressed Due to Excessive Length 3 transport. Examples of these models have been used in [1,10,31,24]. A summary of these different kinds of models can be seen in Table 1. The first equation to model the bed balance was proposed by Felix Exner in his studies river morphology, [16,17]. In this equation, the sediment flux was approximated through flux velocity. In the currently literature different formulations are used to compute the bedload flux, for instance [30,14,21,40]. The majority of formulations proposed consider the critical Shields parameter, [35], to control the bed movement, nevertheless, the Grass law for bedload transport [21] does not consider this parameter, bed movement starts when the flux velocity is different to zero. The Grass law also allows measure if the interaction between bed and flux is strong or weak. In [10], a formulation where some bedload transport models are rewritten in analogous form to law of Grass, is proposed. In order to solve numerically the sediment transport problem, several numerical strategies have been proposed. The majority of them are based on the finite volume numerical schemes, other strategies are based on cellular automata, although the use of cellular automata for modelling the sediment flux is a field that have been few explored, such as is pointed in [23]. As well, recently, the finite elements framework has been employed to compute the numerical solution of the sediment transport on shallow water equation, for instance [15,2,18] have employed this technique. Up to our knowledge the extension of approximate Riemann solvers is the most suitable option for the numerical resolution of the sediment transport, some one of them are exposed in [38] and [25]. In order to solve the complete system for sediment transport modelling, there are two means found in the literature, regardless of whether the model is FCM, PCM or DM. In the first one, the hyperbolic system of balance laws is solved through an approximate Riemann solver scheme. Then, the bed conservation equation is solved by using the information provided by this scheme. This method has been used in works like [8,23,31,24]. In the work [13] is discussed that this way to solve the sediment transport can fail due to the characteristic velocity related to the bedload is not taken into account for to solve the hydrodynamics of water sediment mixture. The second method computes the solution in coupled form, as in [27,3,28,29,33]. In order to compute the depth of bed together with the other variable the bed conservation equation is rewritten introducing a new conservative variable, which depends on the depth of bed and the suspended load. According to Cao [8], the water flux, sediment transport and morphological evolution of the bed are strongly coupled. The rate of deformation can be significant in comparison with the flux evolution, then the use of a fully coupled model that takes into account all conservation laws is needed. In this paper, a FCM is considered, which takes into account both the suspended load sediment transport and the bedload sediment transport. The structure of this work is as follows: in section 2 the mathematical model is presented and the employed formulations for erosion and deposition are described. In section 3 a general compact formulation will be presented, as well as an alternative way 4 Gonz´alez-Aguirre, J. C. et al. to write the model considering the presence of non conservative product ([10, 11]) and the study of the hyperbolicity of the model. The numerical schemes to be used will be shown in section 4, as well as the steady stationary solution which we want to preserve. Finally, in section 5 we validate the mathematical model and the numerical schemes by carrying out several numerical experiments and making comparison with between the obtained numerical results and the experimental data. 2 Mathematical model In this work, an one-dimensional shallow water flow, over an erodible bottom composed of uniform, non-cohesive sediment, is analysed. This phenomenon is depicted in Figure 1. hpx, tq zpx, tq upx, tq Suspended load Bed load Rolling Sliding Jumping Suspended particle   Epx, tq@ @R Dpx, tq Free surface of water Fig. 1 Sketch of bedload and suspended-load sediment transport. The dynamics of the movable bottom as a result of hydrodynamic behaviour can be modelled by using the system of balance laws made up by the mass conservation equation, the momentum conservation equation, the sediment conservation equation and the bed conservation equation, namely, Bh Btpx, tq` Bhu Bxpx, tq “ Epx, tq´Dpx, tq 1´p,(1) Title Suppressed Due to Excessive Length 5 Bphuq Btpx, tq` B`hu2`1 2gh2˘ Bxpx, tq “ ghpx, tqˆ´Bz Bxpx, tq´Sfpx, tq˙ ´δpρs´ρwqgh2px, tq 2ρpx, tqBc Bxpx, tq ´pρ0´ρpx, tqqpEpx, tq´Dpx, tqqupx, tq ρpx, tqp1´pq, (2) Bphcq Btpx, tq` Bphucq Bxpx, tq “ Epx, tq´Dpx, tq,(3) Bz Btpx, tq` 1 1´pBqb Bxpx, tq “ Dpx, tq´Epx, tq 1´p,(4) where –hpx, tqis the depth of water (m), –upx, tqis the averaged-flux velocity (m/s), –cpx, tqis the averaged-flux volumetric sediment concentration (kg/m3), –zpx, tqis the depth of bed (m), –ρpx, tqis the density of water-sediment mixture (kg/m3), –Sfpx, tqis the friction slope computed by using equation (8)), –gis the gravitational acceleration (m/s2), –pis the bed sediment porosity, –Epx, tqis the sediment entrainment (m/s) (see [8]), –Dpx, tqis the sediment deposition (m/s) (see [8]), –ρwis the water density (kg/m3), –ρsis the sediment density (kg/m3), –ρ0is the saturated bed density (kg/m3), –qbpx, tqis the bed-load discharge (m2/s) (see [21]). –δis defined as δ“"1 if θěθc, 0 if θăθc,(5) where θand θcare the Shields parameter and Critical Shields value, respectively. These values are defined in section 2.1. The right-hand site of the mass conservation law for the water sediment mixture (1), is significant for processes where the sediment transport and the morphological evolution are active (see [8]). The second term on the right-hand side of the momentum conservation equation (2), measures the effects of the variable sediment concentration on the flux direction. In [29] is pointed out that non-uniform sediment concentration over irregular bottom in quiescent water may cause perturbation of water at rest state, in order to deal with this issue, the second therm on the right hand is taken into account from the moment that the flux can move the bottom, it is mean, when the Shields parameter is greater than Critical Shield parameter. The third term represents the 6 Gonz´alez-Aguirre, J. C. et al. momentum transfer due to the exchange between the water and the erodible bed, [8]. Equation (3) is the suspended sediment conservation equation. The morphodynamics are well-modelled by the Exner equation (4), which states that the time variation of the sediment layer, in a certain volume, is due to the net variation of the solid transport through of the boundaries of the volume, [16,17]. The formulation of the bed load discharge qbpx, tqcan be based on empirical laws, some of these most often used were proposed in [30], [21] and [40]. In order to compute the bed load discharge, the law of Grass [21] was chosen for this work qbpx, tq “ Agu3px, tq,(6) where Agis a constant value that takes into account the diameter of the particle and the kinematic viscosity, usually it is obtained experimentally. The values for Agrange from zero to one, so that if Agis close to zero, then the model shows a weak interaction between the sediment bottom and fluid. On the other hand, if Agis close to one, then the interaction between the sediment bottom and fluid is strong, [10]. 2.1 Morphological conditions The density of the water-sediment mixture, ρ, and the density of the saturated bed, ρ0, are computed as ρpx, tq “ ρwp1´cpx, tqq`ρscpx, tqand ρ0“ρwp`ρsp1´pq.(7) The friction term is computed using the following formula Sfpx, tq “ cfupx, tq|upx, tq|,(8) where cfis computed as –cf“η2 hpx, tq4 3 , being ηthe dimensionless Manning coefficient. –cfis the friction coefficient. The sediment exchange between the bed and the water column is determined through the sediment entrainment due to turbulence, Epx, tq, and the sediment deposition due to gravitational action, Dpx, tq, (see [8]). The sediment entrainment is given by Epx, tq “ $ ’ & ’ % ϕpθpx, tq´θcqupx, tq hpx, tqd0.2if θpx, tq ě θc, 0 else, (9) the sediment deposition reads Dpx, tq “ αcpx, tqω0p1´αcpx, tqqm,(10) where: Title Suppressed Due to Excessive Length 7 –ϕis a constant dimensionless value determined by ϕ“φ560p1´pqν0.8 3psgq0.4θc .(11) This expression was deducted following to [8]. –φis a dimensionless value that depends on phenomenon to recreate, –νis the kinematic viscosity of the fluid (m2/s), –s“ρs ρw´1 is the submerged specific gravity of the sediment (kg/m3), –θcis the dimensionless critical Shields parameter which indicates the start of sediment movement. –The dimensionless Shields parameter or mobility parameter is computed using the formula θpx, tq “ u2 ˚px, tq sgd ,(12) where –dis the diameter of the particle (10´3m), –u˚px, tqis the friction velocity (m/s), that is given by the relation u˚px, tq “ dτpx, tq ρ, -τpx, tqis the threshold stress of bottom computed using the following expression τ“gρη2upx, tq|upx, tq| hpx, tq1{3. –The drop velocity ω0(m/s) is computed using the formulation proposed in [42] ω0“dˆ13.95ν d˙2 `1.09sgd ´13.95ν d,(13) –mis an exponent computed by using the Reynolds number of the particle, see [26]. m“4.45R´0.1where, R “ω0d ν.(14) –αis the nonequilibrium adaptation coefficient of suspended load, [41]. Following the work [8], this coefficient is defined by α“min t2,p1´pq{cu. This relation holds that near-bed supend-load concentration αc must be smaller than the bed material concentration 1 ´p, [41]. In Table 1 are shown the different kinds of models which have been discussed in section 1, where the one dimensional shallow water equations are give by 8 Gonz´alez-Aguirre, J. C. et al. Bh Btpx, tq` Bhu Bxpx, tq “ 0,(15) Bhu Btpx, tq` Bhu2`1 2gh2 Bxpx, tq“´gh Bz Bxpx, tq.(16) Table 1 Model summary. Model Equations Remarks FCM (1), (2), (3), (4) Complete model used in this work PCM (15), (16), (3), (4) Bedload discharge is neglected DM (15), (16), (4) RHS of (4) is neglected 3 Problem statement In order to write the mathematical model (1)-(3) in a compact manner the following notation is employed –Conservative variables, –w1“hpx, tqis the depth of water (m), –w2“hpx, tqupx, tqis the unit discharge (m2/s), –w3“hpx, tqcpx, tqis the suspended sediment concentration (kg/m2), – Wpx, tq:“ pw1, w2, w3qtis the vector of conservative variables. –Physical flux, FpWq “ ¨ ˚ ˚ ˚ ˚ ˚ ˝ w2 w2 2 w1`1 2gw2 1 w2w3 w1 ˛ ‹ ‹ ‹ ‹ ‹ ‚ .(17) –Source term, Spx, t, Wq “ ¨ ˚ ˚ ˚ ˚ ˚ ˚ ˝ E´D 1´p ´gw1ˆBz Bx`Sf˙´pρs´ρwqgw2 1 2ρB Bxˆw3 w1˙´ρ0´ρ ρ E´D 1´p w2 w1 E´D ˛ ‹ ‹ ‹ ‹ ‹ ‹ ‚ . (18) The source term, Spx, t, Wq, is split into four parts: bed slope, S1px, t, Wq, friction slope, S2px, t, Wq, sediment entrainment and sediment deposition process, S3px, t, Wq, and volumetric sediment concentration slope, Title Suppressed Due to Excessive Length 9 S4px, t, Wq, so that: Spx, t, Wq “ 4 ÿ K“1 Skpx, t, Wq, where S1px, t, Wq “ ¨ ˚ ˚ ˚ ˚ ˝ 0 ´gw1Bz Bxpx, tq 0 ˛ ‹ ‹ ‹ ‹ ‚ ,(19) S2px, t, Wq “ ¨ ˚ ˚ ˚ ˝ 0 ´gw1Sfpx, tq 0 ˛ ‹ ‹ ‹ ‚,(20) S3px, t, Wq “ ¨ ˚ ˚ ˚ ˚ ˚ ˝ Epx, tq´Dpx, tq 1´p ´ρ0´ρpx, tq ρpx, tq Epx, tq´Dpx, tq 1´p w2 w1 Epx, tq´Dpx, tq ˛ ‹ ‹ ‹ ‹ ‹ ‚ ,(21) and S4px, t, Wq “ ¨ ˚ ˚ ˚ ˚ ˝ 0 ´gw2 1pρs´ρwq 2ρpx, tqB Bxˆw3 w1˙ 0 ˛ ‹ ‹ ‹ ‹ ‚ .(22) Using this notation, the equations (1)-(3) can be rewritten in vectorial form as BW Btpx, tq` BFpWq Bxpx, tq “ Spx, t, Wq,(23) for xP r0,Lsand tP r0, Ts, where Lis the length of channel. –Initial conditions: Wpx, 0q “ W0pxq,and zpx, 0q “ z0pxq, x P r0,Ls. In practice, the initial values for averaged depth of water h0pxq, averaged flux velocity u0pxq, averaged flux volumetric sediment concentration c0pxq and depth of bed z0pxq, are known. –Boundary conditions: Conditions at the left end of the channel are detailed and the conditions of the right end are similar 16 Gonz´alez-Aguirre, J. C. et al. Ci Ci´1Ci´1 F´ i´1 2 F` i´1 2 F´ i`1 2 F` i`1 2 Fig. 3 Numerical flux in the cell Cifor HLLCS scheme. relation for steady shock, of shock speed S, at x“0 and t“0 is followed, namely, F`pWpxi`1, tqq ´F´pWpxi, tqq´ 2 ÿ k“1 ¯ Sk´xi`1 2, t, Wpxi`1 2, tq¯ “S`W`pxi`1, tq´W´pxi, tq˘“0, (42) where ¯ Sk(k“1,2), is the averaged approximated source term at x“0 and t“0, so that 2 ÿ k“1 ¯ Skpxi`1 2, t, Wpxi`1 2, tqq “ » – ¯ S1 ¯ S2 ¯ S3fi fl.(43) Due to the source terms ¯ S1px, t, Wqand ¯ S2px, t, Wqare taken into account inside the numerical flux, the hyperbolic system of balance laws (23) can be solved by using the numerical scheme Wpxi, tn`1q “ Wpxi, tnq´ ∆t ∆x ´FHLLCS´ i`1 2´FHLLCS` i´1 2¯ `∆t 4 ÿ k“3 Skpxi, tn,Wpxi, tnqq, (44) where FHLLCS˘ i˘1 2 are the numerical flows with the source terms inside it. The contributions made by the numerical flux in the cells depends on the sign of the eigenvalues at the boundaries, see Figure 3. In order to define the numerical flux we need to have established values λi, λi`1, and λ˘ ˚, which represent a middle wave related to the source terms ¯ S1 and ¯ S2, see Figure 4. These values are given by expressions (49)-(53). In order to compute the numerical flux, we shall extend for the bottom variations the work done in [32], where the numerical flux is computed as: –If 0 ďλi, then FHLLCS´ i`1 2“FpWpxi, tqq, FHLLCS` i`1 2“FpWpxi, tqq` 2 ÿ k“1 ¯ Skpxi`1 2, t, Wpxi`1 2, tqq. (45) Title Suppressed Due to Excessive Length 17 x t (a) Positive middle wave. xixi`1 ¯ S1`¯ S2 ∆t FiFi`1 WiWi`1 W´ i W` i`1 W`` i`1 λiλ` ˚ λi`1 x t (b) Negative middle wave. xixi`1 ¯ S1`¯ S2 ∆t FiFi`1 WiWi`1 W´´ iW´ iW` i`1 λiλ´ ˚λi`1 Fig. 4 HLLCS numerical flux. –If 0 ěλi`1, then FHLLCS´ i`1 2“FpWpxi`1, tqq´ 2 ÿ k“1 ¯ Skpxi`1 2, t, Wpxi`1 2, tqq, FHLLCS` i`1 2“FpWpxi`1, tqq. (46) –If λiď0ďλi`1, then –if λ` ˚ą0, then FHLLCS´ i`1 2“FpWpxi, tqq`λipW´pxi, tq´Wpxi, tqq, FHLLCS` i`1 2“FHLLCS´ i`1 2` 2 ÿ k“1 ¯ Skpxi`1 2, t, Wpxi`1 2, tqq, (47) where W´pxi, tq “ W`pxi`1, tq´ s H` i`1 2 , s H` i`1 2“ ´ ¯ S2 r λ1r λ2¨ ˝1 0 cpxi, tq˛ ‚,W`pxi`1, tq “ w` 1,i`1¨ ˝1 λ` ˚ cpxi, tq˛ ‚, and w` 1,i`1“˜pw1,iqpui´λiq´λis H` i`1 2,1 λ` ˚´λi¸. –If λ´ ˚ă0, then FHLLCS´ i`1{2“FpWpxi`1, tqq`λi`1pW`pxi`1, tq´Wpxi`1, tqq, FHLLCS` i`1{2“FHLLCS´ i`1{2´ 2 ÿ k“1 ¯ Skpxi`1 2, t, Wpxi`1 2, tqq,(48) 18 Gonz´alez-Aguirre, J. C. et al. where W`pxi`1, tq “ W´pxi, tq` s H´ i`1 2 , s H´ i`1 2“ ´ ¯ S2 r λ1r λ2¨ ˝1 0 cpxi`1, tq˛ ‚,W´pxi, tq “ w´ 1,i ¨ ˝1 λ´ ˚ cpxi`1, tq˛ ‚, and w´ 1,i “˜w1,i`1pui`1´λi`1q´λi`1s H´ i`1 2,1 λ´ ˚´λi`1¸. The values λi,λi`1are given by the following expression: λi“#min !r λ1;ui´?gw1,i;ui`1´?gw1,i`1)if |S2| “ 0, r λ1if |S2| ‰ 0.(49) λi`1“#max !r λ2;ui`?gw1,i;ui`1`?gw1,i`1)if |S2| “ 0, r λ2if |S2| ‰ 0.(50) The average values r λ1and r λ2are the Roe mean [34], with expressions (51) ru“ui?w1,i `ui`1?w1,i`1 ?w1,i `?w1,i`1 (51) and, r λ1“r u´cgw1,i `w1,i`1 2,r λ2“r u`cgw1,i `w1,i`1 2, Finally the values for λ˘ ˚take the form λ` ˚“ λiw1,i`1pui`1´λi`1q´λi`1w1,i pui´λiq`λi`1pλis H` i`1 2,1´S1q w1,i`1pui`1´λi`1q´w1,i pui´λiq`λis H` i`1 2,1 , (52) and λ´ ˚“ λiw1,i`1pui`1´λi`1q´λi`1w1,i pui´λiq`λi`1pλis H´ i`1 2,1´S1q w1,i`1pui`1´λi`1q´w1,i pui´λiq`λis H´ i`1 2,1 . (53) Title Suppressed Due to Excessive Length 19 4.3.2 Bed slope approximation In order to get fully described the discretization of bed slope source term, we shall follow to [32], where the approximation for the bed slope is given by Bz Bxpxi`1 2, tnq « g ∆x ˆwn 1,j ´|δz1| 2˙δz1,(54) where j“"iif δz ě0, i`1 if δz ă0,and δz1“$ & % wn 1,i if δz ě0 and diăzn i`1, wn 1,i`1if δz ă0 and di`1ăzn i, δz else, being δz “zpxi`1, tnq´zpxi, tnqand di“w1pxi, tnq`zpxi, tnq. 4.3.3 Friction term approximation The friction term is approximated by using the following relation, (see [32]). ˇˇˇˇ Sf grw1 ∆xˇˇˇˇi`1 2“min ˆˇˇˇˇ Sf grw1 ∆xˇˇˇˇ,ˇˇˇˇru|umin| 2gˇˇˇˇ˙i`1 2 ,(55) where |umin| “ min p|ui|,|ui`1|q,(56) and Sfis computing using equation (8) by choosing cfas the friction coefficient. 4.3.4 Approximation of sediment concentration slope, sediment erosion and sediment deposition process The source terms S3px, tqand S4px, tqin the numerical scheme (44) are approximated as follows 4 ÿ k“3 Sn k,i “ ¨ ˚ ˚ ˚ ˚ ˚ ˚ ˚ ˝ En i´Dn i 1´p ´δgw2 1,i ρs´ρw 2ρi cn i´cn i´1 ∆x ´ρ0´ρ ρ En i´Dn i 1´pun i En i´Dn i ˛ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‚ .(57) 20 Gonz´alez-Aguirre, J. C. et al. 4.4 Numerical scheme for the non-conservative system In order to deal with the presences of non conservative products, a suitable numerical treatment of them is needed. Following [10] the numerical solution of (26) at cell Ciat time tnis computed using a numerical scheme of the form Wn`1 i“Wn i´∆t ∆x ´A` D´Wn i´1 2¯`Wn i´Wn i´1˘ `A´ D´Wn i`1 2¯`Wn i`1´Wn i˘¯`∆t 3 ÿ k“1 Sn k,i, (58) where A˘ D´Wi`1 2¯“1 2´A´Wi`1 2¯˘ˇˇˇA´Wi`1 2¯ˇˇˇ¯, Taking into account the definition of Agiven by equation (29) the numerical scheme (58) becomes in Wn`1 i“Wn i´∆t ∆x ´Gn i`1 2´Gn i´1 2¯ ´∆t 2∆x ´B´Wn i`1 2¯`Wn i`1´Wn i˘`B´Wn i´1 2¯`Wn i´Wn i´1˘¯ `∆t 3 ÿ k“1 Sn k,i (59) where Gn i`1 2“1 2`FpWn i`1q`FpWn iq˘´1 2ˇˇˇA´Wn i`1 2¯ˇˇˇ`Wn i`1´Wn i˘.(60) The discretization of the bed slope source term S1px, tqis carried out following the methodology explained in section 4.2.2, but in this case the upwind process is done with regard to the eigenvalues of matrix A. The source term, which contains the friction slope S2px, tq, is discretized by using either, (38) or (39) while the source term that takes into account the sediment entrainment and deposition processes S3px, tsq, is discretized by using (40). 4.5 Morphological evolution Exner equation (4) is not actually hyperbolic, however it is possible to write a wave speed estimation r λbassociated to sediment flux (see [24]), such that the Exner equation becomes in Bz Btpx, tq`r λbBz Bxpx, tq “ Dpx, tq´Epx, tq 1´p,(61) Title Suppressed Due to Excessive Length 21 where r λb“1 1´pBqb Bz,(62) this wave speed is the basis to build an upwind strategy to solve the morphological evolution. In order to solve the Exner equation (4), we integrate that equation over a space-time arbitrary rectangle Ciˆrtn, tn`1s, as it was done in section 4.1. zn`1 i“zn i`∆t 1´ppDn i´En iq´1 ∆x żtn`1 tn 1 1´p´qbpxi`1 2, tq´qbpxi´1 2, tq¯dt. (63) Following the finite volume methodology, an approximation qn b,i˘1 2 for the integral of bedload discharge at the common boundary xi˘1 2between two finite volume cells, is needed, namely qn b,i˘1 2«1 ∆t żtn`1 tn qbpxi˘1 2, tqdt. (64) Following [24], this approximation is defined as follows qn b,i`1 2 :“#qbpxi, tnqif r λb,i`1 2ą0, qbpxi`1, tnqif r λb,i`1 2ă0, qn b,i´1 2 :“#qbpxi´1, tnqif r λb,i´1 2ą0, qbpxi, tnqif r λb,i´1 2ă0, (65) where r λb,i`1 2“1 1´p qbpxi`1, tnq´qbpxi, tnq zpxi`1, tnq´zpxi, tnq.(66) Finally, the scheme to solve the Exner equation (4) reads zn`1 i“zn i`∆t 1´pˆDn i´En i´1 ∆x ˆqn b,i`1 2´qn b,i´1 2˙˙.(67) 4.6 Numerical treatment of boundary conditions The numerical methods described earlier are applied to internal cells, Ci, i “ 1, . . . , N ´1. At the boundary cells we cannot use these methods because for boundaries nodes xB“x0or xB“xN, left or right neighbouring nodes, which are necessary to define the numerical flux or upwind the source term, do not exist. So as to compute the values of conservative variables at the boundaries, fictitious neighbouring nodes p xB(see [4]) or ghost cells CB, (see[25]), are introduced, see Figure 4.6. The values of the conservative variables in these cells are defined according to the boundary conditions, which can be either inflow, outflow or wall. For 22 Gonz´alez-Aguirre, J. C. et al. -p xB CB x0 x0´1 2 C0 Wn 0 p xB CB xN xN´1 2 Wn N CN xi xi´1 2 xi`1 2 Wn i Ci Fig. 5 Fictitious neighbouring nodes, pxB, and ghost cells, CB. each case, the values of the conservative variables on fictitious neighbouring nodes are: –Inflow (w0,2ą0 or wN,2ă0) –Left end: WppxB, tnq “ Wpx0, tnq,zppxB, tnq “ zpx0, tnq, –Right end: Wpp xB, tnq “ WpxN, tnq,zppxB, tnq “ zpxN, tnq. –Outflow (w0,2ă0 or wN,2ą0) –Left end: WppxB, tnq “ Wpx0, tnq,zppxB, tnq “ zpx0, tnq, –Right end: WppxB, tnq “ WpxN, tnq,zppxB, tnq “ zpxN, tnq. –Wall (w0,2“0 or wN,2“0) –Left end: w1ppxB, tnq “ w1px0, tnq,w2ppxB, tnq“´w2px0, tnq,w3ppxB, tnq “ w3px0, tnq,zppxB, tnq “ zpx0, tnq, –Right end: w1ppxB, tnq “ w1pxN, tnq,w2ppxB, tnq“´w2pxN, tnq,w3ppxB, tnq “ w3pxN, tnq,zpp xB, tnq “ zpxN, tnq. The approximate solution at the boundaries is computed by using the values of the conservative variables on ghost cells and either of the numerical schemes described earlier. 4.7 Steady stationary solution A numerical scheme is called well balanced if it can preserve a family of stationary solution. For that reason the study of the stationary solution is important. In order to get the steady-stationary solution for the system of balance laws (23) and (4), let us suppose that upx, tq “ 0, in this case θăθc, then the system of balance laws reduces to Bh Btpx, tq “ ´1 1´pDpx, tq,(68) Bh Bxpx, tq` Bz Bxpx, tq “ 0,(69) Bhc Btpx, tq“´Dpx, tq,(70) Title Suppressed Due to Excessive Length 23 Bz Btpx, tq “ 1 1´pDpx, tq.(71) From equations (68) and (71) we can see that both, the depth of water and the depth of the bottom may vary, but the free surface of water ςpx, tq “ hpx, tq`zpx, tqstays constant with respect of time, namely Bς Btpx, tq “ Bh Btpx, tq` Bz Btpx, tq “ 0.(72) Therefore, from (69) and (72) is concluded that the free surface of water is a constant with respect of time and space. From equation (70) and by using (68), we get Bc Btpx, tq “ Dpx, tq hpx, tqˆcpx, tq 1´p´1˙.(73) This equation shows that the variation of the volumentric sediment concentration over time is equal to zero if either Dpx, tq “ 0 then cpx, tq “ 0 or cpx, tq “ 1´p. Moreover if cpx, tq ă 1´pthen cpx, tqis a decreasing function over time. Definition 1 The steady stationary solution of balance laws system (23) and (4) is given by the follows conditions –The water free surface is a constant function, see equations (69) and (72). –Null velocity in the whole domain, u“0. –The variation of suspended sediment concentration over time is equal to the opposite of deposition, see equation (70). Proposition 1 Let us suppose that Wpxi, tnqand zpxi, tnqfulfill the steady stationary solutions conditions, then the numerical scheme defined by equations (41), (34), (35), either (38) or (39), (40) and (67), computes the steady stationary solutions conditions exactly, that is mean –The free surface is constant. –The velocity is null over whole domain. –The suspended sediment concentration variation with respect to time is equal to opposite of deposition. Proof The values for the different variables at interface xi˘1{2will be computed as the arithmetic mean. Under water at rest hypothesis, the different elements of numerical scheme take the following forms: 24 Gonz´alez-Aguirre, J. C. et al. –Convective discretization, 1 ∆xpFn i`1 2´Fn i´1 2q “ ¨ ˚ ˚ ˚ ˚ ˚ ˚ ˚ ˚ ˚ ˝ wn 1,i ´wn 1,i´1 2∆x bg wn 1,i´1 2´wn 1,i`1´wn 1,i 2∆x bg wn 1,i`1 2 g 4∆x `pwn 1,i`1q2´pwn 1,i´1q2˘ cn i´1 2 wn 1,i ´wn 1,i´1 2∆x bg wn 1,i´1 2´cn i`1 2 wn 1,i`1´wn 1,i 2∆x bg wn 1,i`1 2 ˛ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‚ . (74) –Discretization of bed slope source term, Sn 1,i “ ¨ ˚ ˚ ˚ ˚ ˚ ˚ ˚ ˚ ˚ ˚ ˚ ˝ 1 2 zn i`1´zn i ∆x bg wn 1,i`1 2´1 2 zn i´zn i´1 ∆x bg wn 1,i´1 2 ´g 2w1,i`1 2 zn i`1´zn i ∆x ´g 2wn 1,i´1 2 zn i´zn i´1 ∆x 1 2ci`1 2bg wn 1,i`1 2 zn i`1´zn i ∆x ´1 2ci´1 2bg wn 1,i´1 2 zn i´zn i´1 ∆x ˛ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‚ . (75) –Discretization of friction slope source term, Sn 2,i “0. –Discretization of sediment entrainment and sediment deposition process source term, Sn 3,i “¨ ˚ ˚ ˚ ˚ ˝ ´Dn i 1´p 0 ´Dn i ˛ ‹ ‹ ‹ ‹ ‚ .(76) –Discretization of sediment concentration slope source term, Sn 4,i “δ¨ ˚ ˚ ˝ Sn p4,iq,1 Sn p4,iq,2 Sn p4,iq,3 ˛ ‹ ‹ ‚,(77) where Sn p4,iq,1“1 2 cn i`1´cn i ∆x bg wn 1,i`1 2 wn 1,i`1 2 ρs´ρw 2ρn i`1 2 ´1 2 cn i´cn i´1 ∆x bg wn 1,i´1 2 wn 1,i´1 2 ρs´ρw 2ρn i´1 2 , Title Suppressed Due to Excessive Length 25 Sn p4,iq,2“ ´ g 2 cn i`1´cn i ∆x ´wn 1,i`1 2¯2ρs´ρw 2ρn i`1 2 ´g 2 cn i´cn i´1 ∆x ´wn 1,i´1 2¯2ρs´ρw 2ρn i´1 2 , Sn p4,iq,3“1 2 cn i`1´cn i ∆x cn 1`1 2 wn 1,i`1 2bg wn 1,i`1 2 ρs´ρw 2ρn i`1 2 ´1 2 cn i´cn i´1 ∆x cn i´1 2 wn 1,i´1 2bg wn 1,i´1 2 ρs´ρw 2ρn i´1 2 . –Exner equation discretization, zn`1 i“zn i`∆t Dn i 1´p Let ςn`1 ibe the approximation of the free surface of water computed by the numerical scheme (41), therefore ςn`1 i“wn`1 1,i `zn`1 i “wn 1,i ´∆t 2ˆwn 1,i ´wn 1,i´1 ∆x bg wn 1,i´1 2´wn 1,i`1´wn 1,i ∆x bg wn 1,i`1 2˙ `∆t 2ˆzn i`1´zn i ∆x bg wn 1,i`1 2´zn i´zn i´1 ∆x bg wn 1,i´1 2˙ `δ∆t 2˜cn i`1´cn i ∆x wn 1,i`1 2bg wn 1,i`1 2 ρs´ρw 2ρn i`1 2 ´cn i´cn i´1 ∆x wn 1,i´1 2bg wn 1,i´1 2 ρs´ρw 2ρn i´1 2¸´∆t Dn i 1´p`zn i`∆t Dn i 1´p. (78) From (78) and taking into account that for this case δ“0, we get ςn`1 i“ςn i´∆t 2„bg wn 1,i´1 2ˆςn i´ςn i´1 ∆x ˙´bgwn 1,i`1 2ˆςn i`1´ςn i ∆x ˙. (79) Bearing in mind that at time t“tnthe free surface is a constant function with respect of x, then we can conclude that ςn`1 i“ςn i, which is mean that the free surface is a constant function. 32 Gonz´alez-Aguirre, J. C. et al. Fig. 8 Free surface of water and depth of bed at t“60 s. Fig. 9 Error in the compute of the free surface of water at t“60s. Figures 13 and 14 show the results obtained for the velocity at t“20 s and t“120 s, respectively. The results provided by the three numerical schemes at t“20 s are very similar to the results of Cao. However, for t“120 s in the region where the water front is computed, the computed velocity is less than the of [8], but the approximation is correct. Both the expansion wave and shock wave shown by Cao are recreated exactly by the three numerical schemes. Likewise, the velocity and morphological changes are estimated correctly. Thus, the three numerical schemes reproduce accurately the numerical results of [8]. Title Suppressed Due to Excessive Length 33 Fig. 10 Velocity at t“60 s. Fig. 11 Error for the compute of temporal evolution of volumetric sediment concentration at t“60 s. 5.4 Experiment 3 In order to measure the accuracy of our numerical schemes, the experimental dam break flows carried out in Louvain-la-Neuve (Universit`e Catholique de Louvain [20]) were numerically simulated and the results obtained were compared with the experimental data. Figures 15 – 17 show the numerical results of the three numerical schemes against the experimental data. The downstream free surface of water computed by schemes QSVLS1S4 and QSVLNCP is greater than the water free surface given for the experimental data in this region. However the approximation is sufficient. Moreover, the results provided by HLLCS scheme are an 34 Gonz´alez-Aguirre, J. C. et al. Fig. 12 Water free surface and bed the depth computed at t“120 s. Fig. 13 Velocity computed at t“20 s. accurate approximation for the downstream water free surface. It should be noted that the up-stream numerical water free surface obtained is less than the experimental data of the water free surface in this region, but the differences are minor. The hydraulic jump observed in the experimental data is better approximated by the HLLCS scheme. Regarding the water front advance, it can be seen that the three numerical schemes are accurate. Regarding the bed evolution, Figures 15 and 16, show that the up-stream morphological evolution computed by the HLLCS scheme, at times 5 pt0and 7.5 pt0, are less than the morphological changes observed in the experimental data. On the other hand, the upstream morphological evolution computed by both schemes QSVLS1S4 and QSVLNCP is in accordance with the experimental data. Furthermore at time 10 pt0both schemes QSVLS1S4 and QSVLNCP Title Suppressed Due to Excessive Length 35 Fig. 14 Velocity computed at t“120 s. compute more upstream sediment entrainment than that observed one in the experimental data, while the upstream sediment entrainment computed by the HLLCS scheme fit better to the experimental data. The downstream sediment entrainment computed by the three numerical schemes is greater than the experimental data for every computed time, although the computed approximation is good enough. Fig. 15 Water free surface and the bed depth computed at 5 p t0. 36 Gonz´alez-Aguirre, J. C. et al. Fig. 16 Water free surface and the bed depth computed at 7.5 p t0. Fig. 17 Water free surface and the bed depth computed at 10 p t0. 5.5 Experiments 4 and 5 In these numerical experiments the conditions of small scale dam breaks, performed in the Universit`e Catholique of Louvain, were recreated. Both the values for physical parameter and the experimental data were obtained from [37]. The results computed in the fourth numerical experiment are shown in Figure 18. It can be seen that the downstream sediment entrainment computed by the HLLCS scheme is greater than the downstream sediment entrainment provided by the experimental data. However, the differences between them are small. The upstream sediment entrainment computed by the HLLCS scheme is in Title Suppressed Due to Excessive Length 37 accordance with the experimental data. The bed evolution computed by both schemes QSVLS1S4 and QSVLNCP reproduce accurately the bed changes. From Figure 18 it can be seen that the region where the hydraulic jump took place is correct for the three numerical schemes, but the experimental water free surface is not given correctly by the numerical results. The downstream conditions for the water free surface are correctly reproduced by the numerical schemes because the numerical results and experimental data have the same qualitative form. Notice that the downstream free surface of water computed by both schemes QSVLS1S4 and QSVLNCP is greater than the experimental water free surface although the adjustment made is adequate. Moreover, the water free surface computed by the HLLCS scheme fits accurately the experimental water free surface. Regarding the water front advance, it can be seen that the results obtained by the HLLCS scheme are faster than the experimental data while the results provided by the other two methods calculated more accurately. The numerical results of the fifth numerical experiment are shown in Figure 19. We can see that the hydraulic jump is only partially reproduced by the three numerical schemes; there is a gap between the place provided by the experimental data and that calculated by the numerical schemes. The computed water free surface after the hydraulic jump is in good agreement with the experimental data. Notice that the water front advance computed by the three numerical schemes presents a gap in comparison with the water front advance observed in the experimental data. Regarding the morphological changes, the tendency observed in the experimental data is correctly represented by the three numerical schemes. The HLLCS scheme produces more upstream sediment entrainment than that observed one in the experimental data, but the differences are minor. The adjustment of the upstream bed evolution made by the schemes QSVLS1S4 and QSVLNCP is satisfactory. The downstream bed evolution computed by the three numerical schemes is in good agreement with the experimental bed evolution. In these numerical tests small oscillations, are noticeable where the hydraulic jump took place, produced by the HLLCS schemes while the results of QSVLS1S4 and QSVLNCP gives a smooth tendency. 6 Conclusions Three numerical schemes to simulate dam break flows with suspended load and bedload sediment transport have been presented. Two of them implement the Q-scheme of van Leer to compute the numerical flux and the other one implements a HLLC Riemann solver. In one of the presented numerical scheme, the source terms related to the bed slope and the sediment concentration variations are discretised in an upwind way. In a second numerical scheme, the bed slope source term is discretised in an upwind way but the source term related to the variation of sediment concentration is treated as a non 38 Gonz´alez-Aguirre, J. C. et al. Fig. 18 Free surface of water and bed depth computed at t“1.25 s. Fig. 19 Free surface of water and bed depth computed at t“1.25 s. conservative product. The last one numerical scheme, the HLLCS Riemann solver, take into account the bed slope source term inside of the numerical flux whiles the sediment concentration variation is computed in a explicit way. We have proved that the numerical schemes are well-balanced, by the numerical results of the steady stationary solution test. The contrast made between the numerical results and the experimental data shows that the numerical schemes accurately compute the water front advance. It can be seen that the computed free surface of water recreates in successfully the free surface of water of the experimental data. Regarding to bed evolution, we have seen that the three numerical scheme compute a successful approximation for the observed bed evolution in the experimental data. Due to that the Grass parameter is close to zero, the interaction between the flux and the bottom is Title Suppressed Due to Excessive Length 39 weak for the physical experiments simulated in this work. From the different numerical tests carried out, we can conclude that the computed water front advance depends on the wet-dry parameter εW D, the threshold value for sediment entrainment εED or the interval pH0,H1q. We have also seen that the constant value φ, in the sediment entrainment formula, helps to improve the computation of the water front advance. Therefore an interesting future work is try to developed an uncertainty quantification study about these parameters. Acknowledgements The authors are grateful to the Consejo Nacional de Ciencia y Tecnolog´ıa for the scholarship granted to carry out this project. This work was partially supported by the Spanish MICINN project MTM2013-43745-R and MTM201786459-R and by the Xunta de Galicia, the FEDER under research project ED431C 2017/60 -014 and was partially supported by PRODEP project UAMPTC-669. References 1. D. D. Apsley and P. K Stansby. Bed-load sediment transport on large slopes: model formulation and implementation within a RANS solver. J. Hydraul. Eng., 134(10):1440– 1451, 2008. 2. Faranak Behzadi and James C Newman III. An exact source-term balancing scheme on the finite element solution of shallow water equations. Comput. Methods. Appl. Mech. Eng., 359:112662, 2020. 3. F. Benkhaldoun, S. Sari, and M Seaid. A flux-limiter method for dam-break flows over erodible sediment beds. Applied Mathematical Modelling, 36(10):4847 – 4861, 2012. 4. A. Berm´udez, X. L´opez, and M. E V´azquez-Cend´on. Numerical solution of nonisothermal non-adiabatic flow of real gases in pipelines. J. Comput. Phys., 323:126 – 148, 2016. 5. A. Berm´udez and M. E V´azquez-Cend´on. Upwind methods for hyperbolic conservation laws with source term. Comput. Fluids, 23(8):1049–1071, 1994. 6. P. Brufau, P. Garc´ıa Navarro, and M. E V´aquez-Cend´on. Zero mass error using unsteady wetting-drying conditions in shallow water flows over dry irregular topography. Int. J. Numer. Methods Fluids, 45(10):1047–1082, 2004. 7. F. N. Cantero-Chinchilla, O. Castro-Orgaz, S. Dey, and J. L Ayuso-Mu˜noz. Nonhydrostatic dam break flows. II: one-dimensional depth-averaged modeling for movable bed flows. J. Hydraul. Eng., 142(12):04016069, 2016. 8. Z. Cao, G. Pender, S. Wallis, and P Carlling. Computational dam-break hydraulics over erodible sediment bed. J. Hydraul. Eng., 130(7):689–703, 2004. 9. H. Capart and D Young. Formation a jump by the dam-break wave over a granular bed. J. Fluid Mech., 372:165–187, 1998. 10. M. J. Castro-D´ıaz, E. D. Fern´andez-Nieto, and A. M Ferreiro. Sediment transport models in shallow water equations and numerical approach by high order finite volume methods. Comput. Fluids, 37(3):299 – 316, 2008. 11. M. J. Castro-D´ıaz, E. D. Fern´andez-Nieto, T. Morales de Luna, G. Narbona-Reina, and C Par´es. A HLLC scheme for nonconservative hyperbolic problems. application to turbidity currents with sediment transport. Esaim Math. Model Numer. Anal., 47:1–32, 2013. 12. L Cea. An unstructured finite volume model for unsteady turbulent shallow water flow with wet-dry fronts: numerical solver and experimental validation. PhD thesis, Universidade da Coru˜na, 2005. 40 Gonz´alez-Aguirre, J. C. et al. 13. S. Cordier, M.H. Le, and T Morales de Luna. Bedload transport in shallow water models: Why splitting (may) fail, how hyperbolicity (can) help. Adv. Water. Resour., 34(8):980 – 989, 2011. 14. H. A Einstein. The bed load function for sediment transportation in open channel flows. Bulletin 1026, Department of Agriculture, Soil Conservation Service, Washington, D.C., 1950. 15. Fateme Ebrahimi Erami and Ali Rahmani Firoozjaee. Numerical solution of bed load transport equations using discrete least squares meshless (dlsm) method. Applied Mathematical Modelling, 77:1095–1109, 2020. 16. F. M Exner. Zur physik der d¨unen. Akad. Wiss. Wien Math. Naturwiss, Klasse, 129:929–952, 1920. 17. F. M Exner. ¨ Uber die wechselwirkung zwischen wasser und geschiebe in fl¨usen. Akad. Wiss. Wien Math. Naturwiss, Klasse, 134:165–2014, 1925. 18. Ali Rahmani Firoozjaee and Mostafa Sahebdel. Element-free galerkin method for numerical simulation of sediment transport equations on regular and irregular distribution of nodes. Eng. Anal. Bound. Elem., 84:108–116, 2017. 19. L. Fraccarollo and E. F Toro. Experimental and numerical assessment of the shallow water model for two-dimensional dam-break type problems. J. Hydraul. Res., 33(6):843– 864, 1995. 20. L. Fraccarrollo and H Capart. Riemann wave description of erosional dam-break flows. J. Fluid Mech., 461:183–238, 2002. 21. A. J Grass. Sediment transport by waves and currents. Technical Report FL29, SERC. London Cent. Mar, 1981. 22. A Harten. On class of high resolution total-variation-stable finite-difference schemes. SIAM J. Numer. Anal., 21(1):1–23, 1984. 23. P. Hu and Z Cao. Fully coupled mathematical modeling of turbidity currents over erodible bed. Adv. Water. Resour., 32(1):1 – 15, 2009. 24. C. Juez, J. Murillo, and P Garc´ıa-Navarro. A 2D weakly-coupled and efficient numerical model for transient shallow flow and movable bed. Adv. Water. Resour., 71:93 – 109, 2014. 25. R. J LeVeque. Finite-Volume methods for hyperbolic problems. Cambridge University Press., 2002. 26. J. Li, Z. Cao, G. Pender, and L Qingquan. A double layer-averaged model for dam-break flows over mobile bed. J. Hydraul. Res., 51(5):518–534, 2013. 27. S. Li and C. J Duffy. Fully coupled approach to modeling shallow water flow, sediment transport, and bed evolution in rivers. Water Resour. Res., 47(3):W03508, 2011. 28. X. Liu, J. A. Infante Sedano, and A Mohammadian. A robust coupled 2-D model for rapidly varying flows over erodible bed using central-upwind method with wetting and drying. Can. J. Civ. Eng., 42(8):530–543, 2015. 29. X. Liu, A. Mohammadian, A. Kurganov, and J. A Infante Sedano. Well-balanced central-upwind scheme for a fully coupled shallow water system modeling flows over erodible bed. J. Comput. Phys., 300:202 – 218, 2015. 30. E. Meyer-Peter and R M¨uller. Formulas for bed-load transport. Technical report, 2nd Meet. Int. Assoc. Hydraul. Struct. Res., Stockholm, 1948. 39-64. 31. J. Murillo and P Garc´ıa-Navarro. An exner-based coupled model for two-dimensional transient flow over erodible bed. J. Comput. Phys., 229(23):8704 – 8732, 2010. 32. J. Murillo and P Garc´ıa-Navarro. Augmented versions of the HLL and HLLC riemann solvers including source terms in one and two dimensions for shallow flow applications. J. Comput. Phys., 231(20):6861–6906, 2012. 33. Khawar Rehman and Yong-Sik Cho. A novel well-balanced scheme for spatial and temporal bed evolution in rapidly varying flow. J. Hydro-environ. Res., 27:87–101, 2019. 34. P. L Roe. A basis for upwind differencing of the two-dimensional unsteady Euler equations. Numerical Methods for Fluid Dynamics II, pages 55–80, 1986. 35. A Shields. Anwendung der Aehnlichkeitsmechanik und der Turbulenzforschung auf die Geschiebebewegung. PhD thesis, University Berlin, 1936. 36. S. Soares-Fraz˜ao and Y Zech. HLLC scheme with novel wave-speed estimators appropriate for two-dimensional shallow-water flow on erodible bed. Int. J. Numer. Methods Fluids, 66(9):1019–1036, 2011. Title Suppressed Due to Excessive Length 41 37. B. Spinewine and Y Zech. Small-scale laboratory dam-break waves on movable beds. J. Hydraul. Res., 45(sup1):73–86, 2007. 38. E. F Toro. Shock-capturing methosd for free-surface shallow flows. Wiley & Sons, Ltd, 2001. 39. E. F. Toro, M. Spruce, and W Speares. Restoration of the contact surface in the HLLRiemann solver. Shock waves, 4(1):25–34, 1994. 40. L. C Van Rijn. Sediment transport, part I: Bed load transport. J. Hydraul. Eng., 110(10):1431–1456, 1984. 41. W. Wu and S. S Wang. One-dimensional modeling of dam-break flow over movable beds. J. Hydraul. Eng., 133(1):48–58, 2007. 42. R. Zhang and J Xie. Sedimentation research in China: Systematic selection. China and Water and Power Press, Beijing, 1993.