Bed Evolution under Rapidly Varying Flows by a New Method for Wave Speed Estimation

Abstract: This paper proposes a sediment-transport model based on coupled Saint-Venant and Exner equations. A finite volume method of Godunov type with predictor-corrector steps is used to solve a set of coupled equations. An efficient combination of approximate Riemann solvers is proposed to compute fluxes associated with sediment-laden flow. In addition, a new method is proposed for computing the water depth and velocity values along the shear wave. This method ensures smooth solutions, even for flows with high discontinuities, and on domains with highly distorted grids. The numerical model is tested for channel aggradation on a sloping bottom, dam-break cases at flume-scale and reach-scale with flat bottom configurations and varying downstream water depths. The proposed model is tested for predicting the position of hydraulic jump, wave front propagation, and for predicting magnitude of bed erosion. The comparison between results based on the proposed scheme and analytical, experimental, and published numerical results shows good agreement. Sensitivity analysis shows that the model is computationally efficient and virtually independent of mesh refinement.


Introduction
Numerical techniques for modeling morphodynamics are preferable over field and laboratory observations at large scales, due to reduced computational costs and time [1][2][3][4]. Mathematical modeling of morphodynamics can be performed using depth-integrated shallow-water equations (continuity and momentum equations) representing the hydrodynamics, and an Exner equation (continuity of sediment) representing the morphological evolution. These equations can either be modeled using a coupled approach in which the hydrodynamic and morphodynamics flows are solved simultaneously, or using an uncoupled (quasi two-phase) approach, in which flow hydrodynamics are modeled first and then the obtained velocity is used for modeling morphological evolution [5]. According to [6,7], the coupled approach is best suited for analysis of rapidly varying flow conditions, where flow unsteadiness strongly affects sediment evolution; while the uncoupled approach is mostly used for slowly varying flows. The applicability of the uncoupled approach is limited due to the assumption of constant sediment concentration in the lower layer [8].
Coupled shallow-water and Exner equations have been widely used to model sediment transport in rivers, bed evolution in dam-break flows, siltation of reservoirs, and sediment flushing in sewers [4,[9][10][11][12][13][14]. In this study, the authors used a coupled model based on approaches adopted by Goutière, Soares-Frazão et al. (2008) and Murillo and García-Navarro (2010) to simulate a bed load sediment-transport initiated by water flow. Different equations are used in the literature to model solid transport sediment flux. Commonly used sediment-transport formulas include the Grass formula [15], Meyer-Peter and Muller's (MPM) equation [16], Van Rijn's equation [17][18][19], and Nielsen's formula [20]. The performance of various transport formulas is reported in [21]. For simplicity, the authors have adopted Grass's model in this study to define a solid transport discharge. This is a simple sediment formula that combines a global case dependent calibration and an easily differentiable power law of velocity [22].
We use the volume method of Godunov-type to solve coupled shallow-water and Exner equations models. Two approaches are usually employed for the computation of interface fluxes in coupled models, the Roe approximate Riemann solver [5,10,[22][23][24][25][26][27], and the Harten-Lax-van Leer-Contact (HLLC) approximate Riemann solver [7,14]. A predictor-corrector step is adopted for evolving the solution to get second-order temporal accuracy. In the predictor step, the solution is updated by half a time step using a non-conservative method; while in the corrector step, the conservation is recovered over a full time step by solving local Riemann problems using the data obtained from the predictor step. The fluxes at the interface are computed using the HLLC approximate Riemann solver for the hydrodynamic part, and by the Harten-Lax-van Leer (HLL) scheme for the morphological changes.
Since this study employs the Riemann problem, the general solution for the hydrodynamic part of the solution consists of three waves, a left and right moving wave and a middle wave (a shear wave). In such a problem, the most important feature of flux computation with the HLLC scheme is the choice of relations for approximating the water depth and velocity along the shear wave. The usual practices of computing these variables are either based on the assumption that the two non-linear waves are either shocks or rarefactions, or by using depth positivity condition (readers are referred to [28] for more details). However, in cases of coupled shallow-water and Exner equations, these equations are not applicable due to differences in wave structure.
Wave speed estimation in HLLC solver requires the values of Roe-averaged states (water depth, velocity, energy states) [29]. These states are obtained from Roe approximations (see details in [30]). Given that the Roe-averaged states are not uniquely defined [31], different approached are adopted for estimating Roe-averaged states. Most of the numerical models for sediment transport studies [6,13,14,22,23] follow Glaister's approach [32] for computing Roe averages for a coupled shallow-water Exner equations system. However, Hu et al. [33] argued that artificially introduced intermediate states, like that of Glaister's, may not be admissible for some cases. Therefore, they proposed an extension of the Roe averages for estimation of wave speeds by avoiding artificial states and relying only on intermediate states for state averages.
The HLLC studies [7,14] used the average of left and right Riemann states for computation of variables along the shear wave, but this method makes the solution too dependent on reconstructed values, which can be erroneous for cases in which the grids are skewed or when the problem presents strong discontinuities, and when the reconstruction procedure is not accurate. We propose a new method for approximating the values along the shear wave by Roe averages in which the variable values do not exclusively depend upon reconstructed values. The basic concept is to predict average states directly from the intermediate states by utilizing vertex values at the interface along with reconstructed variables. To make the dependence more even, the variable values at vertex are included in finding the Roe averages. This method avoids potential problems of inaccuracy and instability, which can be caused by meshes with skewed elements, problems with strong discontinuities, and may address the shortcomings of the reconstruction procedure to some extent.
The aim of this work is to develop a simple, but robust, coupled shallow-water flow model that can accurately handle the discontinuities and capture shocks associated with dam-break flows on erodible beds. The objective is to develop a model that can help us study the bed evolution under the sudden release of large volume of water and to accurately capture the wave structure associated with such flows.

Governing Equations
The governing equations describe the Saint-Venant system commonly known as shallow-water equations coupled to the Exner equation for sediment-transport and forming a non-linear system of equations.

Hydrodynamic Component
Flow evolution can be represented by the shallow-water equations. For an inviscid and incompressible flow, these equations can be expressed as: where t is the time, h is the water depth, u is the velocity component, g is the gravitational acceleration, and S and τ are the bed slope and bottom frictional terms, respectively, estimated as: where b is bed elevation, and n is Manning's roughness coefficient.

Morphodynamic Component
The sediment dynamics can be expressed by a simple sediment continuity or Exner equation. Mathematically, this equation is given as [12]: where ξ " 1{1´p, p is the sediment porosity and q s is sediment-transport discharge. Physically, Equation (3) implies that time variation of a sediment layer in a control volume is due to the net variation of solid transport through the boundaries of the control volume [22]. In this study, the Exner equation is used in its reduced form, as the main focus is bed evolution and flow dynamics. A sediment discharge transport formula is used to close Equation (3). To check the performance of the proposed model, we compute bed load discharge using two approaches, the bed load equation of Grass [15] and the MPM sediment transport formula.

Grass Formula
Sediment transport flux in the Grass formula follows a power law, given as: where A is the intensity of interaction between flow and sediment particles, the value of which ranges between 0 and 1 with increasing interaction strength, and u is the velocity. This equation does not take critical shear stress into account, and it relates the sediment discharge to the flow velocity as in [12]. For the cases modeled in this study, A is either adopted from previous studies, or calculated using Equation (5) for cases with detailed sediment parameters; A " dˆ0.005´d h¯0 where d is the local sediment medium diameter; ρ s and ρ are the sediment and water density, respectively; and the particle diameter D * is estimated as where υ is the kinematic viscosity of water.

MPM Formula
The MPM formula is given as [16]: where s is the specific gravity of sediment, θ is the Shield's parameter and θ c is threshold Shield's parameter, set equal to 0.047 in our study, following [14]. Inserting the expression of Shield's parameter into Equation (7) forms

Numerical Scheme
In conservative form, Equations (1) and (3) can be written as: where U is the vector of conserved variables, F is a flux vector functions, and S is a vector of source terms. In vector form, the coupled equations can be written as [13]: Integration of Equation (9) over the cell gives where Ω is the volume of a control volume. Equation (11) can be expressed as: where n denotes the unit outward normal to the boundary of the cell, and F(U)¨n is the flux vector normal to the interface, expressed as: Discretization of Equation (12) takes the form: where ∆t and ∆x are step sizes in space and time, respectively, and x i´1{2 and x i`1{2 are the boundaries of the cell. The explicit nature of this scheme demands that the time step is limited by the CFL (Courant-Friedrichs-Lewy coefficient) condition.
where ν is the CFL number, set equal to 0.6 for all cases modeled in this study. Equation (14) is evolved in two steps: the predictor step and the corrector step. In the predictor step, the solution evolves by half a time step employing reconstructed cell center variables. This step is given as: where F i`1{2 and F i´1{2 are fluxes at the right and left sides of the cell interface, respectively, n indicates the initial values of variables, and U n i is the vector of conserved variables at the cell center at initial time n. Fluxes F i˘1{2 are computed by the HLLC and HLL schemes for the hydrodynamic and morphodynamic parts, respectively. The results from Equation (17) are advanced to a full time step in the second stage by the corrector scheme given as:

Flux Computation
Flux calculation using coupled shallow-water Exner equations is complicated due to the additional wave, related to sediment information in the wave structure. Approximate Riemann solvers (especially Roe and HLLC schemes) can accurately estimate fluxes for shallow-water flows, as these solvers have high resolution and shock capturing capability. However, for cases with strong interface interactions, like dam-break flows over erodible beds, complications arise due to the additional wave (associated with the sediment layer) issued from an initial discontinuity in the coupled equations.
We use HLLC and HLL schemes in this study to compute hydrodynamic and morphological fluxes, respectively. The wave speeds for the coupled system of equations, Equation (7), are estimated from the approximate eigenvalues derived in [13] by analysis of the Jacobian matrix for the coupled system of equations.
where c " a gh is the wave celerity, and d is a the term involving the normal derivative of sediment transport rate. Differentiating the sediment discharge, given in Equations (4) and (8), gives us: The largest eigenvalue λ 3 accounts for hydrodynamic effects and λ 1,2 give eigenvalues carrying sediment characteristic information. The HLLC scheme used for computing hydrodynamic mass and momentum fluxes is given as [34]: The speed of the middle wave is λ˚, and is approximated as where K = L, R; S L " min`λ 1L , λ 1R , u˚´agh˚˘, S R " max`λ 3L , λ 3R , u˚`agh˚˘, h˚and u˚are variable values along the shear wave estimated by Equations (29) and (30), respectively. The eigenvalue analysis in [13] shows that sediment flow information (λ 1,2 ) is not affected by the fast moving wave (λ 3 ), since the speed with which sediment information is transmitted is very slow compared to the evolution of other hydrodynamic information. We choose the HLL scheme for computing sediment fluxes because of the constant intermediate state between λ 1 and λ 2 , which is the assumption of HLL scheme. Application of the HLL scheme is simple, as intermediate waves between the slowest and fastest moving waves are not present in this approach. The speed of these waves is estimated by eigenvalues of the coupled equation; the speed of the slowest and fastest travelling waves is estimated by Equation (21). The HLL flux is calculated as [34]: The important step in arriving at the solution is to approximate the water depth and velocity across the shear wave. Their constant values, in Equations (29) and (30), are computed by a new method in which the contribution of vertex values is also accounted for u˚" where v represents the two vertices of the interface through which the flux is being computed.
The vertex values can be computed by the pseudo-Laplacian relation proposed by [35]. This new method for approximation of variable values in the star region also gives weight to the variable values at the vertices. This method is effective for highly skewed control volumes, as well as for structured grids, as the variable approximation does not strongly depend upon the reconstructed values at cell interfaces. In previous studies, the estimation of shear waves is based solely on Roe averages that are obtained by solving locally linearized Riemann problem at each cell interface [12,30,36]: This method might make the solution susceptible to error for cases in which the mesh is composed of highly distorted grids, which lead to cell-centered variables that are not reconstructed exactly at the mid-point of an interface. We show the performance of the proposed method in Section 3 by modeling various test cases on different meshes.

Source Term Treatment
Accuracy of the numerical solution is greatly influenced by the source term treatment. In the present study, source terms include variable bed slope and shear stress terms. The scheme is well-balanced if its flux gradient is balanced by the slope source terms for a "lake at rest" condition where η is water surface elevation. We extend and modify the source term treatment, the revised surface gradient method (RSGM), proposed by [37], for shallow water flows on non-erodible beds. We extend and modify RSGM for erodible beds by introducing reconstructed values of water surface elevation to estimate flow depth average and bed level difference.
S " The proposed scheme satisfies well-balanced conditions and gives stable and accurate solutions as shown in application of the model. Friction terms were treated explicitly.

Wetting and Drying
Inaccurate treatment of wetting and drying fronts results in significant errors in solutions, including negative water depths and unphysical high velocities. A simple treatment is adopted in this study, which involves setting a tolerance depth h min (0.0001 m in our model), below which cells are considered dry, and consequently the momentum flux across the cell interface becomes zero. This treatment of wet-dry fronts is similar to applying solid boundary conditions on the interface between wet and dry cells.

Boundary Conditions
Wall and open boundary conditions are imposed in the proposed model. Unlike the pure hydrodynamic case, there are three characteristic waves entering and leaving the domain at the inlet and outlet. The additional wave accounts for the presence of sediment. Solid wall conditions are imposed on side walls making the velocity normal to the walls equal to zero. Physical boundary conditions are required at the inlet. For cases in which some physical conditions are not available, the theory of characteristics is used. The inlet boundary conditions are implemented by the same procedure as literature references, including [23,38]. A free outflow condition is imposed at the outlets to enable waves to pass through the boundary without reflection.

Dam-Break over a Rigid Bed
The proposed scheme is first tested for a benchmark hydrodynamic case of dam-break on a rigid non-erodible bed widely employed to verify model's performance on dry and wet beds. The dam-break is instantaneous on a flat bottom with negligible friction. Based on the water level in the right half of the domain, there are two cases for this problem; if water level is over zero the case is said to be dam-break on a wet domain, and if it is zero the case is dry bed. The initial conditions for this problem constitute a Riemann problem. For initially wet bed the conditions are given as The initial conditions for dry bed are where L is the total length of the channel taken equal to 200 m, and x o " 100 m is the position of dam.
The results are compared with analytical solution [39] at t = 1, 12, 21 and 27 s for flow over wet and bed in Figure 1. The numerical results show excellent agreement with the analytical solution; the transition of wet/dry front is accurately captured showing that the scheme is well-balanced. Proposed solution shows excellent agreement with analytical solution where the results are almost indistinguishable over both wet and dry beds in Figure 1a,b, whereas the scheme without modification also gives very satisfactory predictions with slight difference from the analytical solution when flow changes from subcritical to supercritical; the difference, however, is indistinguishable on the plotted scale. We observe that the water level estimated by proposed scheme gives slightly better results than the scheme without modification. conditions are required at the inlet. For cases in which some physical conditions are not available, the theory of characteristics is used. The inlet boundary conditions are implemented by the same procedure as literature references, including [23,38]. A free outflow condition is imposed at the outlets to enable waves to pass through the boundary without reflection.

Dam-Break over a Rigid Bed
The proposed scheme is first tested for a benchmark hydrodynamic case of dam-break on a rigid non-erodible bed widely employed to verify model's performance on dry and wet beds. The dam-break is instantaneous on a flat bottom with negligible friction. Based on the water level in the right half of the domain, there are two cases for this problem; if water level is over zero the case is said to be dam-break on a wet domain, and if it is zero the case is dry bed. The initial conditions for this problem constitute a Riemann problem. For initially wet bed the conditions are given as The initial conditions for dry bed are where L is the total length of the channel taken equal to 200 m, and  Figure 1a,b, whereas the scheme without modification also gives very satisfactory predictions with slight difference from the analytical solution when flow changes from subcritical to supercritical; the difference, however, is indistinguishable on the plotted scale. We observe that the water level estimated by proposed scheme gives slightly better results than the scheme without modification. (a)

Bed-Slope Aggradation
This test case models a sediment-laden flow over a sloping bed that emulates flows in mountainous terrains triggered by torrential rains. A uniform supercritical flow with constant sediment supply enters the domain from the left boundary. Bed evolution is observed until the equilibrium state is reached by gradual aggradation of the bed. Boundary conditions at x = 0 and L are set equal to s q and 0, respectively. In order to compare our results with [13]-based on similar governing equations to this study-we used the MPM formula to calculate sediment discharge with a mean sediment diameter of 1.65 mm and porosity (p) of 0.42. The analytical solution of this problem was presented by [40], given as: where ∞ denotes the equilibrium state; An is the Fourier Series coefficient, estimated by Equation (39); and 0 T is the response time for the problem given as: The channel characteristics are: length (L) = 6.90 m; width = 0.5 m; bed roughness = 0.0165 s·m −1/3 ; and slope S0 = 2.4%. The flow parameters are fixed as: water discharge q = 7.98 m 3 /s and sediment discharge upstream qs = 0.147 m 3 /s. These parameters give a response time of T0 = 1670 s and a final equilibrium slope, S, equal to 3.03%. The results are compared with the analytical and numerical solution in Figure 2. After start of the flow at t = 0, bed evolution goes through a transient state (observed at t = 900 and 1800 s) until the final equilibrium state is reached at t = 3600 s. Good agreement is observed between the computed and analytical solutions at all time instants. Our scheme gives slightly better approximations for aggradation than the results of [13]. This application shows the applicability of our scheme for modeling sloping rivers. This also demonstrates the capability of the model to compute the equilibrium slope based on aggradation and degradation of channels, which is useful for limiting sedimentation of water bodies.

Bed-Slope Aggradation
This test case models a sediment-laden flow over a sloping bed that emulates flows in mountainous terrains triggered by torrential rains. A uniform supercritical flow with constant sediment supply enters the domain from the left boundary. Bed evolution is observed until the equilibrium state is reached by gradual aggradation of the bed. Boundary conditions at x = 0 and L are set equal to q s and 0, respectively. In order to compare our results with [13]-based on similar governing equations to this study-we used the MPM formula to calculate sediment discharge with a mean sediment diameter of 1.65 mm and porosity (p) of 0.42. The analytical solution of this problem was presented by [40], given as: where 8 denotes the equilibrium state; A n is the Fourier Series coefficient, estimated by Equation (39); and T 0 is the response time for the problem given as: The channel characteristics are: length (L) = 6.90 m; width = 0.5 m; bed roughness = 0.0165 s¨m´1 /3 ; and slope S 0 = 2.4%. The flow parameters are fixed as: water discharge q = 7.98 m 3 /s and sediment discharge upstream q s = 0.147 m 3 /s. These parameters give a response time of T 0 = 1670 s and a final equilibrium slope, S, equal to 3.03%. The results are compared with the analytical and numerical solution in Figure 2. After start of the flow at t = 0, bed evolution goes through a transient state (observed at t = 900 and 1800 s) until the final equilibrium state is reached at t = 3600 s. Good agreement is observed between the computed and analytical solutions at all time instants. Our scheme gives slightly better approximations for aggradation than the results of [13]. This application shows the applicability of our scheme for modeling sloping rivers. This also demonstrates the capability of the model to compute the equilibrium slope based on aggradation and degradation of channels, which is useful for limiting sedimentation of water bodies.

Dam-Break in an Erodible Channel
We simulate a dam-break case on a mobile bed to assess the performance of proposed scheme for high speed flows associated with rapid bed deformation and high speed wave propagation. This case was previously also tested by [41,42] to verify their models which were also based on coupled shallow-water and Exner equations.
The computational domain consists of a frictionless erodible channel bed of 50 m length with a dam at the middle of the channel. Following [41,42] the initial water depth of 1 m is fixed upstream of the dam and downstream is considered dry. Computations of sediment discharge by Grass formula and MPM formula have no significant effects on the numerical results for bed evolution, water depths, and velocity [41,42]. We used Grass formula to compute sediment discharge with sediment porosity (p) of 0.4 and A is taken equal to 0.004. At time t = 0, the dam is instantaneously removed; the consequent water level, bed deformation, and velocity is compared with exact Riemann solution of the problem in Figure 3 (readers are referred to [42] for details on exact Riemann solution).
For comparison with the exact solution, the results for water level and bed deformation are non-dimensionalized with respect to initial water depth in the channel, and velocity is non-dimensionalised by wave celerity. The results are plotted at 2 s interval up to 5 s. The numerical results show the correct propagation of wave-front, bed evolution, and velocity. Shocks introduced by the bed discontinuity are well captured at all instants and the numerical solution is free of oscillations. Numerical solution (proposed and without modification) show lag to exact solution at all time instants, which was also observed by [41]. In Figure 3b, we observe that numerical solution slightly overestimates bed deformation with time advancement, which is a characteristic feature of computing sediment discharge by Grass formula due to neglecting water depth in flux estimation. Overall, flow variables show very good agreement with exact solution with the proposed scheme giving better results in all cases.

Dam-Break in an Erodible Channel
We simulate a dam-break case on a mobile bed to assess the performance of proposed scheme for high speed flows associated with rapid bed deformation and high speed wave propagation. This case was previously also tested by [41,42] to verify their models which were also based on coupled shallow-water and Exner equations.
The computational domain consists of a frictionless erodible channel bed of 50 m length with a dam at the middle of the channel. Following [41,42] the initial water depth of 1 m is fixed upstream of the dam and downstream is considered dry. Computations of sediment discharge by Grass formula and MPM formula have no significant effects on the numerical results for bed evolution, water depths, and velocity [41,42]. We used Grass formula to compute sediment discharge with sediment porosity (p) of 0.4 and A is taken equal to 0.004. At time t = 0, the dam is instantaneously removed; the consequent water level, bed deformation, and velocity is compared with exact Riemann solution of the problem in Figure 3 (readers are referred to [42] for details on exact Riemann solution).
For comparison with the exact solution, the results for water level and bed deformation are non-dimensionalized with respect to initial water depth in the channel, and velocity is non-dimensionalised by wave celerity. The results are plotted at 2 s interval up to 5 s. The numerical results show the correct propagation of wave-front, bed evolution, and velocity. Shocks introduced by the bed discontinuity are well captured at all instants and the numerical solution is free of oscillations. Numerical solution (proposed and without modification) show lag to exact solution at all time instants, which was also observed by [41]. In Figure 3b, we observe that numerical solution slightly overestimates bed deformation with time advancement, which is a characteristic feature of computing sediment discharge by Grass formula due to neglecting water depth in flux estimation. Overall, flow variables show very good agreement with exact solution with the proposed scheme giving better results in all cases.

Dam-Break over Sand Bed
Experimental studies for sand bed evolution under dam-break waves were conducted at the Civil Engineering Department of Universite Catholique de Louvain (UCL), Belgium [43]. The experimental results were later validated by numerical simulations in references [8,22,44,45]

Dam-Break over Sand Bed
Experimental studies for sand bed evolution under dam-break waves were conducted at the Civil Engineering Department of Universite Catholique de Louvain (UCL), Belgium [43]. The experimental results were later validated by numerical simulations in references [8,22,44,45] as well as by others.
The tests were conducted for idealized dam-break cases over an erodible bed. The experimental set-up included a flume 6 m long, 0.25 m wide and 0.70 m deep. A thin steel gate, fixed in the center of flume, was used as a dam in the experiments. In the original experiments [43], the erosion was observed for two bed materials, uniform coarse sand and Polyvinyl Chloride (PVC) pellets. For brevity purposes, this study simulates the experiment with a sand bed, extending on both sides of the dam. The thickness of the sand bed varied for different configurations to account for bed discontinuity at dam location, and these configuration were modeled in the study.
The following [43] four different cases are simulated in this study, which primarily differ in upstream bed layer thickness and initial water depth. Following previous studies [22], the values of the dimensional constant (A), sediment porosity (p), and Manning's coefficient (n) are set equal to 0.01, 0.4, and 0.0165, respectively, for all configurations. Figure 4 shows the initial configurations for which the dam-break flow is simulated using the proposed scheme. The domain was discretized into 400 cells. The model was run using a refined mesh, but no significant improvement in the results was observed (discussed further in Section 4.5). A CFL number of 0.6 was used due to the explicit nature of the scheme. The dam was removed instantaneously at t = 0, and the results were obtained for the four water surface and bed profiles. Model results were compared to the experimental measurements of Spinewine and Zech (2007) and to the numerical scheme without the modification introduced in Equations (31) and (32). pellets. For brevity purposes, this study simulates the experiment with a sand bed, extending on both sides of the dam. The thickness of the sand bed varied for different configurations to account for bed discontinuity at dam location, and these configuration were modeled in the study.
The following [43] four different cases are simulated in this study, which primarily differ in upstream bed layer thickness and initial water depth. Following previous studies [22], the values of the dimensional constant (A), sediment porosity (p), and Manning's coefficient (n) are set equal to 0.01, 0.4, and 0.0165, respectively, for all configurations. Figure 4 shows the initial configurations for which the dam-break flow is simulated using the proposed scheme. The domain was discretized into 400 cells. The model was run using a refined mesh, but no significant improvement in the results was observed (discussed further in Section 4.5). A CFL number of 0.6 was used due to the explicit nature of the scheme. The dam was removed instantaneously at t = 0, and the results were obtained for the four water surface and bed profiles. Model results were compared to the experimental measurements of Spinewine and Zech (2007) and to the numerical scheme without the modification introduced in Equations (31) and (32).  We observe that the model without modification overestimates bed erosion particularly immediately after the dam-break, and is not able to follow the experimental wave-front, whereas the proposed scheme gives very good agreement with the erosion and the wave-fronts observed in the experiment. The position of wave-fronts given by the proposed scheme at t = 0.25 s and 0.75 s are in good agreement with previous results, while at t = 1.25 s it is slightly above the measured, showing slow movement of the wave, due to bed friction. As water propagates downstream, more particles come in contact with the bed surface compared with the initial time intervals resulting in slower propagation of the wave-front. However, considering human error in experimental measurements, and considering the limitations of numerical modeling, we can conclude that, overall, the agreement between numerical results and experimental results using the proposed scheme is very satisfactory.  We observe that the model without modification overestimates bed erosion particularly immediately after the dam-break, and is not able to follow the experimental wave-front, whereas the proposed scheme gives very good agreement with the erosion and the wave-fronts observed in the experiment. The position of wave-fronts given by the proposed scheme at t = 0.25 s and 0.75 s are in good agreement with previous results, while at t = 1.25 s it is slightly above the measured, showing slow movement of the wave, due to bed friction. As water propagates downstream, more particles come in contact with the bed surface compared with the initial time intervals resulting in slower propagation of the wave-front. However, considering human error in experimental measurements, and considering the limitations of numerical modeling, we can conclude that, overall, the agreement between numerical results and experimental results using the proposed scheme is very satisfactory.  Figure 6 shows the comparison of water surface and bed profiles for configuration 2 at t = 0.5 s. Although the bed profile is in fair agreement with experimental results, the wave-front in the present study slightly leads the measured. The leading wave-front may be due to the excessive erosion of the bed step, which, according to Evangelista et al. (2013), occurs due to a greater angle of internal friction than a vertical slope of the bed step. The angle of internal friction is not considered in the Grass formula, which may be the reason for the slightly lower bed scour estimates and the leading wave fronts.  Figure 6 shows the comparison of water surface and bed profiles for configuration 2 at t = 0.5 s. Although the bed profile is in fair agreement with experimental results, the wave-front in the present study slightly leads the measured. The leading wave-front may be due to the excessive erosion of the bed step, which, according to Evangelista et al. (2013), occurs due to a greater angle of internal friction than a vertical slope of the bed step. The angle of internal friction is not considered in the Grass formula, which may be the reason for the slightly lower bed scour estimates and the leading wave fronts.  Figure 7. In Figure 7a, the numerical studies under-predict the bed scour; this is due to the sudden removal of the gate in experiments, which causes excessive soil collapse at the start of dam-break. However, at later instants, when the bed layer smoothens, the results of the proposed scheme are in good agreement with the experimental results as compared to the scheme without modification, as shown in Figure 7b,c. The differences in position of wave fronts likely occur for the same reason discussed for configuration 1.   Figure 7. In Figure 7a, the numerical studies under-predict the bed scour; this is due to the sudden removal of the gate in experiments, which causes excessive soil collapse at the start of dam-break. However, at later instants, when the bed layer smoothens, the results of the proposed scheme are in good agreement with the experimental results as compared to the scheme without modification, as shown in Figure 7b,c. The differences in position of wave fronts likely occur for the same reason discussed for configuration 1.  Figure 7. In Figure 7a, the numerical studies under-predict the bed scour; this is due to the sudden removal of the gate in experiments, which causes excessive soil collapse at the start of dam-break. However, at later instants, when the bed layer smoothens, the results of the proposed scheme are in good agreement with the experimental results as compared to the scheme without modification, as shown in Figure 7b,c. The differences in position of wave fronts likely occur for the same reason discussed for configuration 1.

Sensitivity Analysis
A sensitivity analysis is performed to study the effects of mesh resolution on computational time, flow variables (discharge), and bed elevation. The model is simulated for five consecutively refined mesh systems; the number of elements from the coarsest to finest mesh is 400, 1000, 2000, 3500, and 5000. The analysis was performed on an Intel ® Core TM i7-3770 (8 M Cache, 3.4 GHz) computer with 8 GB DDR3 RAM. Figure 9 shows the increase in simulation time with mesh refinement for configurations 1 and 4. It can be observed in Figure 9 that there are no significant changes in water depth and velocity with mesh refinement. This shows that the proposed scheme satisfies the primary requirement of robust solvers to be computationally efficient and independent of the mesh resolution. Sensitivity of the sediment parameters, including Manning's roughness (n) and porosity (p) are also analyzed for the laboratory scale test discussed in Section 4.4. We consider configurations 1 and 4 for comparison purposes. Figure 10 shows the variation in results as each parameter varies with the rest kept constant. In Figure 10a,b, no significant changes are observed in bed erosion with changes in sediment parameters, which again can be attributed to the presence of downstream water. Figure 10b shows that with changes in bed roughness, the hydraulic jump extends beyond the measured value (shown in Figure 8b), which is due to the lagging wave front causing a delay in the formation of the jump.

Sensitivity Analysis
A sensitivity analysis is performed to study the effects of mesh resolution on computational time, flow variables (discharge), and bed elevation. The model is simulated for five consecutively refined mesh systems; the number of elements from the coarsest to finest mesh is 400, 1000, 2000, 3500, and 5000. The analysis was performed on an Intel ® Core TM i7-3770 (8 M Cache, 3.4 GHz) computer with 8 GB DDR3 RAM. Figure 9 shows the increase in simulation time with mesh refinement for configurations 1 and 4. It can be observed in Figure 9 that there are no significant changes in water depth and velocity with mesh refinement. This shows that the proposed scheme satisfies the primary requirement of robust solvers to be computationally efficient and independent of the mesh resolution.

Sensitivity Analysis
A sensitivity analysis is performed to study the effects of mesh resolution on computational time, flow variables (discharge), and bed elevation. The model is simulated for five consecutively refined mesh systems; the number of elements from the coarsest to finest mesh is 400, 1000, 2000, 3500, and 5000. The analysis was performed on an Intel ® Core TM i7-3770 (8 M Cache, 3.4 GHz) computer with 8 GB DDR3 RAM. Figure 9 shows the increase in simulation time with mesh refinement for configurations 1 and 4. It can be observed in Figure 9 that there are no significant changes in water depth and velocity with mesh refinement. This shows that the proposed scheme satisfies the primary requirement of robust solvers to be computationally efficient and independent of the mesh resolution. Sensitivity of the sediment parameters, including Manning's roughness (n) and porosity (p) are also analyzed for the laboratory scale test discussed in Section 4.4. We consider configurations 1 and 4 for comparison purposes. Figure 10 shows the variation in results as each parameter varies with the rest kept constant. In Figure 10a,b, no significant changes are observed in bed erosion with changes in sediment parameters, which again can be attributed to the presence of downstream water. Figure 10b shows that with changes in bed roughness, the hydraulic jump extends beyond the measured value (shown in Figure 8b), which is due to the lagging wave front causing a delay in the formation of the jump.   Sensitivity of the sediment parameters, including Manning's roughness (n) and porosity (p) are also analyzed for the laboratory scale test discussed in Section 4.4. We consider configurations 1 and 4 for comparison purposes. Figure 10 shows the variation in results as each parameter varies with the rest kept constant. In Figure 10a,b, no significant changes are observed in bed erosion with changes in sediment parameters, which again can be attributed to the presence of downstream water. Figure 10b shows that with changes in bed roughness, the hydraulic jump extends beyond the measured value (shown in Figure 8b), which is due to the lagging wave front causing a delay in the formation of the jump. For configuration 1, the variation in porosity (0.35-0.38) and Manning's roughness (0.03-0.04) show significant increases in bed erosion, as seen in Figure 11a,b, respectively. The increase in values of both parameters result from increased interaction of sediment and water related to direct contact of water particles with sediment, due to absence of tail-water. This scenario also results in a higher resistance offered to the flow in comparison to configuration 4. We conclude that the absence of tail-water makes the erodible bed more sensitive to changes in sediment parameters, and can significantly affect the bed evolution under wave-front propagation, as observed in Figure 11.  For configuration 1, the variation in porosity (0.35-0.38) and Manning's roughness (0.03-0.04) show significant increases in bed erosion, as seen in Figure 11a,b, respectively. The increase in values of both parameters result from increased interaction of sediment and water related to direct contact of water particles with sediment, due to absence of tail-water. This scenario also results in a higher resistance offered to the flow in comparison to configuration 4. We conclude that the absence of tail-water makes the erodible bed more sensitive to changes in sediment parameters, and can significantly affect the bed evolution under wave-front propagation, as observed in Figure 11. For configuration 1, the variation in porosity (0.35-0.38) and Manning's roughness (0.03-0.04) show significant increases in bed erosion, as seen in Figure 11a,b, respectively. The increase in values of both parameters result from increased interaction of sediment and water related to direct contact of water particles with sediment, due to absence of tail-water. This scenario also results in a higher resistance offered to the flow in comparison to configuration 4. We conclude that the absence of tail-water makes the erodible bed more sensitive to changes in sediment parameters, and can significantly affect the bed evolution under wave-front propagation, as observed in Figure 11.

Conclusions
This paper introduces a simple and robust model for bed evolution and its effects on water surface profile under rapidly varying flows using a formulation comprised of coupled Saint-Venant and Exner equations. A finite volume method of Godunov type is adopted to obtain solutions for the coupled system. The proposed method for the calculation of variables in the star region has proven to give better results than the previous numerical studies, and is able to reproduce the experimental results very satisfactorily. The sensitivity analysis shows that the proposed scheme is computationally very efficient, which is important for schemes modeling sediment flow since the wave structure involved is usually complicated. It was observed for all cases that the values of bed friction, porosity, and dimensional constant (A), which represents the intensity of interaction between water and sediment particles, significantly influence the stability and accuracy of the model, and the position of the hydraulic jump. Analytical solutions, original experimental data, and numerical results are adopted for comparisons. The accuracy of the model is tested for a number of cases, highlighting its efficiency and adaptability under various conditions.

Conclusions
This paper introduces a simple and robust model for bed evolution and its effects on water surface profile under rapidly varying flows using a formulation comprised of coupled Saint-Venant and Exner equations. A finite volume method of Godunov type is adopted to obtain solutions for the coupled system. The proposed method for the calculation of variables in the star region has proven to give better results than the previous numerical studies, and is able to reproduce the experimental results very satisfactorily. The sensitivity analysis shows that the proposed scheme is computationally very efficient, which is important for schemes modeling sediment flow since the wave structure involved is usually complicated. It was observed for all cases that the values of bed friction, porosity, and dimensional constant (A), which represents the intensity of interaction between water and sediment particles, significantly influence the stability and accuracy of the model, and the position of the hydraulic jump. Analytical solutions, original experimental data, and numerical results are adopted for comparisons. The accuracy of the model is tested for a number of cases, highlighting its efficiency and adaptability under various conditions.