This paper presents an efficient second-order finite volume method for the simulation of bedload sediment transport which is capable of modeling wet/dry fronts. The governing equations comprise the shallow water equations (SWEs) for the hydrodynamic phase and the Exner equation for the estimation of bedload sediment transport. These sets of equations are then solved using a weakly coupled scheme based on an augmented Riemann solver (WCAR). In this approach, first, the morphodynamic equation is solved, and then, updated bedload changes with the same Riemann structures are used as a source term within the SWEs. The Smart formula is utilized to estimate the bedload sediment discharge within the Exner equation. The proposed numerical model is first used to model a parabolic sediment layer. Then, a bedload hump propagation with an initial subcritical condition is considered. Next, the simulation of dam-break flow over a mobile bed is investigated. Finally, the dam failure due to over-topping is studied and the computed results are compared with available experimental data. Numerical results indicate that the introduced weakly coupled method, developed based upon the augmented Riemann solver, can be effectively used for modeling all investigated flow regimes, including dry-state interfaces.

Accurate understanding and modeling of sediment transport in free-surface flows and the interactions between bedload and flow are crucial for the prediction of morphological changes in rivers, seas and estuaries. Numerical methods employed to simulate morphodynamic systems should be able to evaluate the evolution of sediment transport and predict water surface changes. Generally, these methods comprise the shallow water equations (SWEs) for the modeling of free-surface flow and the Exner equation which approximates the bedload sediment transport (Delis & Papoglou 2008; Castro Díaz et al. 2009) and is considered as an equilibrium sediment transport approach. These two sets of equations usually lead to a nonlinear system of hyperbolic conservation laws and are mainly solved by coupled or splitting approaches (Canestrelli et al. 2010).

The equilibrium sediment transport models for the morphodynamic systems in particular for the rapidly varying transient flows sometimes provide less accurate results and therefore the sediment transport models based on non-equilibrium states have been adopted (Wu et al. 2004a). Non-equilibrium formulations are suitable for strong bedload and free-surface interactions. These models mainly comprise spatial and temporal lag effects between transport capacity and local hydrodynamic conditions (Wu et al. 2018).

Decoupling or splitting techniques solve the SWEs and the Exner equation in an entirely separate matter (Wu et al. 2004b). In this approach, the Exner equation is first solved and the resulting bed profile update is used as a source term for the SWEs. In this method no additional computation is required for the hydrodynamic phase. However, the main weakness of this scheme, as reported in the literature (Saiedi 1997), is the appearance of instabilities caused by the use of different time steps for each equation. In addition, this method is only applicable for weak or mild interactions that take place between bedload sediment transport and free-surface flow (Wu 2008).

Despite being susceptible to numerical instabilities, finite volume-based splitting approaches have been conventionally adopted for simulating morphodynamic systems. Cordier et al. (2011) presented a well-balanced splitting approach and used it to simulate a classical shallow water system. The results of their model have been compared with those computed from a coupled approach and showed that the splitting approach often produces unphysical instabilities which may be avoided in some cases by reducing the Courant–Friedrichs–Lewy (CFL) condition.

Coupled techniques solve the entire morphodynamic system in a single time step which leads to more stability within the solution. Moreover, any strong interactions that may occur in transcritical and supercritical regimes can be readily modeled (Canestrelli et al. 2010; Cordier et al. 2011). However, the method requires a Jacobian matrix for the entire morphodynamic system which consists of complicated sets of eigenvalues and eigenvectors. These eigenvalues and eigenvectors are mainly associated with the choice of the bedload sediment discharge formula with variable interaction parameter and are often difficult to approximate. To avoid this issue some explicit expressions for both coupled eigenvalues and eigenvectors have been developed for the one-dimensional (1D) problems (Castro Díaz et al. 2008; Murillo & García-Navarro 2010).

Different numerical methods have been developed for solving coupled morphodynamic systems, among which the Godunov-type finite-volume methods have been widely employed (Castro Díaz et al. 2008; Serrano-Pacheco et al. 2012). Delis & Papoglou (2008) used an upwind relaxation scheme for the solution of bedload sediment transport based upon the coupled strategy. Wu & Wang (2008) developed a 1D sediment transport coupled model which considered both suspended and bedload sediment discharge in the computations. Goutière et al. (2008) introduced a one-dimensional bedload sediment transport model using the HLL Riemann solver. Soares-Frazão & Zech (2011) defined a coupled system based on the HLLC scheme and chose a different pair of eigenvectors for approximating the flux. In another work, Rosatti & Fraccarollo (2006) employed a well-balanced method with a new strategy for the treatment of non-conservative fluxes. This work was later modified by Murillo & García-Navarro (2010) who defined a novel approximate coupled Jacobian matrix (CJM) method using a triangular mesh. Although the proposed model was accurate, it was proven to be computationally expensive for real morphodynamic problems (Juez et al. 2014). To overcome the computational drawback, Juez et al. (2014) introduced a first-order scheme which mainly solves the morphodynamic and hydrodynamic phases in a weakly coupled form but with less accurate results. Franzini & Soares-Frazão (2018) presented a coupled strategy based upon the augmented Roe solver. In another work, Martínez-Aranda et al. (2019) introduced a 1D coupled Roe scheme adapted to arbitrary cross-section channels. Finally, Mahdizadeh & Sharifi (2019) defined the two-dimensional bedload model based upon a modified flux-wave formula.

The aim of this paper is to approximate morphodynamic equations using a Godunov-type finite volume method, known as the Godunov wave propagation algorithm, following a weakly coupled scheme and based on an augmented Riemann solver firstly introduced by George (2008) for wet/dry front propagation. The weakly coupled method is a form of a splitting method that holds the features of coupled methods. Unlike other splitting methods that are dependent on the Courant number and flow regime, and also fail to model strong interactions between water flow and sediments, a weakly coupled scheme developed based on an augmented Riemann solver, with new speed approximations for Riemann waves and wave decomposition, is able to accurately simulate sediment transport over dry interfaces. To the best of the author's knowledge, no extension of the augment Riemann solver has been used for bedload sediment transport modeling. To obtain a basic understanding the problem is considered one-dimensional.

The remainder of the paper is organized as follows. First, the background to the governing equations of one-dimensional morphodynamic systems and different empirical formulas for sediment transport is presented. Then, a brief explanation of the wave propagation algorithm and the augmented Riemann solver for a weakly coupled system is given. Finally, in order to validate the proposed method, several test cases with wet/dry states are simulated, and numerical results are compared with analytical and reference solutions.

### Governing equations

The general governing equations for approximating sediment transport include the SWEs and the Exner equation (Exner 1925):
(1)
which can be also re-rewritten in the conservation law form as
(2)
where
(3)
In the above equations, g is the gravitational acceleration (m/s2), h is the fluid depth (m), u is the depth-averaged velocity in the x-direction (m/s), B is the height of the sediment layer (m), is the bedload sediment discharge in the x-direction (m2/s), where represents the porosity of the bedload materials, which is non-dimensional, with, is the bed shear stress calculated by
(4)
where is the water density and is the bed friction coefficient which can be evaluated by Manning's equation:
(5)
So far, many analytical laws have been developed for determining sediment discharge among which the Grass formula follows a simple power law given as
(6)
where is the Grass coefficient that is usually determined from different empirical formulations.
The sediment transport rate, , is often expressed by the following dimensionless parameter:
(7)
where is the ratio between sediment material () and water densities, and is the median diameter of the sediment material (mm). The sediment transport rate is a function of the dimensionless bed shear stress or Shields parameter which is usually used in a dimensionless form:
(8)

Different empirical formulae to evaluate are presented in Table 1.

Table 1

Various empirical equation for determining parameter

In the equations listed in Table 1, the threshold is called the critical Shields parameter. It depends on the physical properties of the sediment layer and is usually computed empirically. It should be noted that sediment transport occurs when exceeds.

By substituting (7) in (6), the following expression for the Shields parameter is obtained:
(9)
Using Equations (4) and (7), the bedload transport formula can be expressed as
(10)
with and varying in each case as shown in Table 3 (Juez et al. 2013).

In Table 2, d30 and d90 are the grain diameters for which 30% and 90% of the weight of a non-uniform sample is finer, respectively, and S0 is the bed slope.

Table 2

Sediment transport formulas in the form of and

Formula
Meyer-Peter & Müller (1984)   0.0470
Camenen & Larson (2005)   0.0400
Fernandez Luque & Van Beek (1976)   0.037–0.0455
Smart (1984)   0.0470
Formula
Meyer-Peter & Müller (1984)   0.0470
Camenen & Larson (2005)   0.0400
Fernandez Luque & Van Beek (1976)   0.037–0.0455
Smart (1984)   0.0470
Table 3

Comparison between the RMSE computed between first- and the second-order WCAR methods and the approximation solution

Cells
B50100200400800
RMSE First order 0.1064 0.0812 0.0548 0.0331 0.0181
0.0818 0.0436 0.0223 0.0128 0.0104
Second order 1 0.0437 0.0232 0.0107 0.0146 0.0178
5 0.0494 0.0311 0.0239 0.0172 0.0047
Cells
B50100200400800
RMSE First order 0.1064 0.0812 0.0548 0.0331 0.0181
0.0818 0.0436 0.0223 0.0128 0.0104
Second order 1 0.0437 0.0232 0.0107 0.0146 0.0178
5 0.0494 0.0311 0.0239 0.0172 0.0047

### Wave propagation algorithm

The wave propagation algorithm described by LeVeque (1998) implements a Godunov-type finite volume method for the solution of Riemann problems by calculating updated versions of the unknowns for the next calculation step and has been successfully applied for many hyperbolic conservation laws problems such as gas dynamics (Bale et al. 2002), SWEs (George 2008; Mahdizadeh et al. 2011, 2012), and water-hammer equations (Mahdizadeh 2019; Mahdizadeh et al. 2018). The Godunov-type wave propagation algorithm used to solve the morphodynamic system can be written as
(11)
where is the vector of unknowns, and are the time and spatial steps, indicates the updated version for the vector of unknowns, represent the second-order correction fluxes required to achieve the high-order solution. If,, then the first-order Godunov method is obtained. and are left- and right-going fluctuations, respectively, and can be obtained by solving the Riemann problems at each cell interface and can be defined for the cell interfaces:
(12)
where is called the flux-wave propagating with a wave speed described in a later section.

### Second-order correction terms

To improve the order of accuracy for the wave propagation algorithm the second-order correction terms can be used. The second-order correction terms are mainly determined by the flux-waves at the cell interface given by
(13)
where is a limited version of flux-wave, , evaluated as where implies the limiter function and . The index I is utilized to represent the upwind side at the cell interface :
(14)
Then is computed using a choice of limiter (LeVeque 2002) with the monotonized centered (MC) limiter used (LeVeque 2002; Mahdizadeh 2010) in this work:
(15)

### Wave decomposition and flux-wave method

For the hyperbolic conservation laws, LeVeque & Pelanti (2001) introduced a wave decomposition and flux-wave formula for the wave propagation algorithm in the following form:
(16)
where and are the fluxes at the right and left side of the cell interface , and and denote the vector of unknowns at the same locations, m is the number of waves propagating from each Riemann interface, is the flux-wave and denotes the wave (LeVeque 2002). For the wave propagation algorithm, each wave can be related to the flux-wave through the flowing equation (LeVeque 1998, 2002):
(17)
where again represents the wave speed. It should be stated that the first component of Equation (13) gives the wave decomposing formula for the wave propagation algorithm and the second equation provides the flux-wave approach which can be also described in the following form:
(18)

Numerical approximation of bedload sediment transport by a weakly coupled approach.

In order to predict the evolution of the bed profile by a weakly coupled approach developed based upon the augmented Riemann solver, shallow water and Exner equations are solved separately. For this purpose, the sediment transport formula is first solved for a specific time step. Then, the resulting bed profile is inserted as a fixed source term into the momentum equation for the shallow water system. Finally, SWEs are solved using the flux-wave approach. In this research, for the first time, the flux-wave method is used in order to approximate the sediment transport equation. Thus, Equation (16) is rewritten in the form of (19) for the sediment transport formula:
(19)
where and are the amount of sediment discharge for i and i − 1 cells computed by the Grass-type equations. In cases where the differences between the bed profile of adjacent cells tends to zero, i.e. , it can be presumed that the interaction between the sediment and water surface is rather weak and, the depth-average velocity, u, in Equation (6) should be substituted by where and q is the constant water discharge (steady-state condition). Thus the Grass-type equation can be rewritten as
(20)
and the wave speed can be computed as
(21)
To evaluate waves within each computational finite volume cell interface, only the differences between adjacent bed profile cells are required. Therefore, the flux wave for sediment transport formula is obtained by
(22)
According to the flux wave method, can be simply obtained and the bed variation update for the next time step can be computed as
(23)

### Augmented shallow water system

In this section, the augmented Riemann solver developed by George (2008) is introduced for solving the SWEs over mobile topography. The advantage of the augmented Riemann solver is that it incorporates desirable characteristics of different Riemann solvers such as Roe, HLLE and exact Riemann solvers. Therefore, this approach produces depth positivity and well-balanced results in the Riemann solution for wet/dry surfaces. The conservation law shallow water system given in (3) can be augmented by multiplying with the Jacobian matrix, such as
(24)
which can be rewritten as the following equation by replacing :
(25)
For the momentum equation, the above equation becomes:
(26)
By considering Equation (26), the SWEs expressed in (3) excluding the friction term is equivalent to the following system:
(27)
where the new vector of unknowns and corresponding Jacobian matrix are:
(28)
The new conservation quantities vector contains depth, , momentum,, momentum flux,, and bed level, B. The Jacobian matrix has the following eigenvectors and eigenvalues (George 2008):
(29)
Therefore, the corresponding wave decomposition provided in (16) for the system of (26) for two neighboring finite volume cells becomes:
(30)
where m = 4. In the case that the friction terms included within the solution, the wave decomposition becomes:
(31)

As can be observed in Equations (30) and (31) the differences between the bathymetry deviations for the neighboring finite-volume cells have been considered as a fourth wave for the augmented Riemann solver. This provides an outstanding advantage over the flux-wave approach and other conventional Riemann solvers in particular in terms of wet/dry front propagation with the existence of the mobile bed.

### The choice of the wave speed for the flux-wave approach

The Roe (1981) linearization method was proposed for obtaining parameter vectors for Euler equations. This scheme can be applied for deriving a parameter average in each cell for the SWEs. The Roe-average approximation of fluid depth can be easily obtained by taking the arithmetic average over and:
(32)
and the Roe-average for velocity can be calculated as
(33)
Consequently, the Roe speeds become
(34)
For dry and nearly dry states, Roe speeds give nonphysical solutions to the SWEs. This problem can be partially rectified by using Einfeldt speeds (Einfeldt 1988):
(35)
where is the kth eigenvalue for the Jacobian matrix and implies the kth eigenvalue of the Roe matrix. Using Einfeldt speeds with the flux-wave approach can provide depth non-negativity, but cannot be applied to shallow water problems when strong transonic rarefaction waves arise in the solution and in particular where strong rarefaction waves collide with dry states, Einfeldt speeds cannot provide the exact solution. This shortcoming can be overcome by using another choice of speed for the flux-wave approach which is called the Generalized Einfeldt speed (George 2008). This new wave speed contains a combination of the Einfeldt and exact Riemann speeds:
(36)

### Approximate augmented Riemann solver

A Riemann problem for the augmented system provided in (28) with piecewise constant initial data is given as
(37)
where
(38)
An approximate solution to the Riemann problem (37) is expressed by decomposing the Riemann waves:
(39)
where represent the pth Riemann wave propagating in each cell interface and are local approximations to the eigenvectors of
(40)
The values are the Generalized Einfeldt speeds defined in the previous section. Also and are the averages defined according to the following:
(41)
By subtracting the fourth eigenvector of relation (40) from the vector quantities Equation (36), the Riemann problem reduces to a three-equation system. Thus, the Riemann waves are presented by solving a linear system for calculating the coefficients, and:
(42)

### Stability condition

In order to achieve numerical stability, the CFL condition (CFL) should be applied. For the one-dimensional wave propagation algorithm the following CFL condition should be fulfilled:
(43)
with .

### Numerical results

This section validates the accuracy of the numerical scheme proposed and described in the previous section. Numerical results have been compared with exact solutions or experimental data. First, the method is used to simulate parabolic sediment layer transportation over a flat bed. Then, sediment hump propagation is modeled with the initial subcritical flow. Next, one-dimensional dam break on moving bed including three scenarios is considered. Finally, dam failure due to over-topping flow is simulated. For all cases, the morphodynamic model was solved using an in-house FORTRAN code on an Intel Core i7-5500U 2.4 GHz with 8GB of RAM.

### Parabolic sediment layer transportation

The purpose of this problem is to compare the numerical solution obtained by the introduced augmented Riemann solver and Grass model, and the analytical solution provided by (Hudson & Sweby 2005). In this test case, water depth and flow rates are set at constant values of 10 m and 10 m3/s, respectively. The initial condition for sediment layer profile is defined as (Hudson & Sweby 2005):
(44)

By setting the values , , and B = 1, the calculated time t = 238,079 s for the approximate solution is obtained (Hudson 2001). Figure 1 shows the initial bed profile and compares the numerical results obtained by the weakly coupled scheme using augmented Riemann solver (WCAR) and the exact solution. These results have been evaluated at t = 238,079 s and 5,339 s (which correspond to B= 1 and 5, respectively) with 200 computational cells and CFL condition is set equal to 0.95.

Figure 1

Comparison of results obtained by the weakly coupled scheme using augmented Riemann solver (WCAR) and the exact solution for the parabolic sediment propagation test case at (a) t = 238,079 s and B = 1 and (b) t = 5,339 s and B = 5.

Figure 1

Comparison of results obtained by the weakly coupled scheme using augmented Riemann solver (WCAR) and the exact solution for the parabolic sediment propagation test case at (a) t = 238,079 s and B = 1 and (b) t = 5,339 s and B = 5.

Close modal

According to Figure 1(a) and 1(b), a very good sediment layer localization is obtained for both values of B, especially in the upper part of the moving hump where a large difference is observed in the predictions of numerical methods developed based on non-conservative fluxes (Castro Díaz et al. 2008). Moreover, Figure 2 represent the difference between the first- and the second-order versions WCAR with B = 1 and B = 5, respectively.

Figure 2

Comparison between first- and second-order weakly coupled scheme using augmented Riemann solver with (a) B = 1 and (b) B = 5.

Figure 2

Comparison between first- and second-order weakly coupled scheme using augmented Riemann solver with (a) B = 1 and (b) B = 5.

Close modal
Figure 3

Initial condition for the sediment hump propagation with initial subcritical flow test case.

Figure 3

Initial condition for the sediment hump propagation with initial subcritical flow test case.

Close modal

As shown in Figure 2, first-order results are more diffusive than the second order scheme. As expected higher order better approximates the sediment layer thickness. Additionally, Table 3 demonstrates the root mean square error (RMSE) calculated between the first- and second-order WCAR methods and the approximate solutions with different computational cells. As can be observed the second-order scheme produces a smaller amount of RMSE compared to the first-order method.

### Sediment hump propagation with initial subcritical flow

This case study is presented by Cordier et al. (2011) and evaluates the ability of the proposed numerical solver in modeling the propagation of the sediment hump with different values of which creates different interactions between sediment and water surface. The initial condition for the discharge, bed profile and depth-average velocity were selected as (Figure 3):
(45)
Figure 4

Sediment hump propagation in a subcritical regime with Ag = 0.005 at t = 10 s.

Figure 4

Sediment hump propagation in a subcritical regime with Ag = 0.005 at t = 10 s.

Close modal

To create the middle interaction between the sediment layer and fluid, first, the Grass formula with was considered. The numerical results obtained at t = 10 s, using 200 computational cells and CFL = 0.95, are illustrated in Figure 4.

According to Figure 4, the proposed model is well able to simulate the shock waves caused by the movement of the sediment layer. Figure 5 shows the simulation results by the WCAR model with CFL = 0.95, t = 2.1 s and a different value of .

Figure 5

Sediment hump propagation in a subcritical regime with Ag = 0.07 at t = 2.1 s.

Figure 5

Sediment hump propagation in a subcritical regime with Ag = 0.07 at t = 2.1 s.

Close modal

It should be noted that the value of leads to a strong interaction between bedload sediment and free-surface flow, which creates nonphysical results (Cordier et al. 2011), but the proposed method does not create this problem and provides a very good approximation of the solution.

The Meyer-Peter and Müller sediment discharge formula is used to evaluate the efficiency of the proposed numerical model. For this purpose, the initial discharge was set equal to hu = 0.8 m2/s and the values for Manning's roughness, sediment density and the CFL number were chosen to be n = 0.095, ρs = 2,612 kg/m3 and CFL = 0.95, respectively. The value of is no longer constant, and it varies continually with the water depth according to Table 2. Figure 6 demonstrates numerical results at t = 2.1 s using WCAR.

Figure 6

Sediment hump propagation in a subcritical regime using Meyer-Peter formula t = 2.1 s.

Figure 6

Sediment hump propagation in a subcritical regime using Meyer-Peter formula t = 2.1 s.

Close modal

As can be seen, WCAR is able to simulate the shock waves created for bedload, as well as free-surface flow without any numerical oscillations.

Furthermore, a comparison was performed between the results computed by the proposed weakly coupled scheme using augmented Riemann solver (WCAR) and a coupled flux-wave method (CFW), proposed by Mahdizadeh & Sharifi (2019). Table 4 demonstrates the Euclidean norm between the solutions computed using WCAR and CFW models and the reference solutions obtained with 10,000 computational cells with at t = 10 s. Also, in Table 5, a similar comparison is done with at t = 2.1 s.

Table 4

Comparison between the Euclidean norms of the difference between WCAR and CFW simulations and the reference solution achieved with 10,000 numerical cells at time t = 10 s with Ag = 0.005

Cells
50100200400800
WCAR h 0.08554 0.05217 0.03101 0.01486 0.00448
hu 0.01227 0.00792 0.00518 0.00277 0.00137
B 0.05818 0.03536 0.02209 0.01107 0.00327
CFW h 0.05312 0.04556 0.02519 0.01049 0.00297
hu 0.01153 0.00505 0.00191 0.00113 0.0092
B 0.02036 0.01478 0.0157 0.00286 0.00182
Cells
50100200400800
WCAR h 0.08554 0.05217 0.03101 0.01486 0.00448
hu 0.01227 0.00792 0.00518 0.00277 0.00137
B 0.05818 0.03536 0.02209 0.01107 0.00327
CFW h 0.05312 0.04556 0.02519 0.01049 0.00297
hu 0.01153 0.00505 0.00191 0.00113 0.0092
B 0.02036 0.01478 0.0157 0.00286 0.00182
Table 5

Comparison between the Euclidean norms of the difference between WCAR and CFW simulations and the reference solution achieved with 10,000 numerical cells at time t = 2.1 s and Ag = 0.07

Cells
50100200400800
WCAR h 0.09034 0.0593 0.03977 0.02474 0.0139
hu 0.08826 0.05809 0.03989 0.0256 0.01547
B 0.06482 0.03999 0.03367 0.02578 0.01636
CFW h 0.053 0.01593 0.01127 0.01597 0.00915
hu 0.02559 0.02789 0.02096 0.01082 0.01092
B 0.03185 0.01365 0.0223 0.01283 0.01006
Cells
50100200400800
WCAR h 0.09034 0.0593 0.03977 0.02474 0.0139
hu 0.08826 0.05809 0.03989 0.0256 0.01547
B 0.06482 0.03999 0.03367 0.02578 0.01636
CFW h 0.053 0.01593 0.01127 0.01597 0.00915
hu 0.02559 0.02789 0.02096 0.01082 0.01092
B 0.03185 0.01365 0.0223 0.01283 0.01006

As can be observed, the WCAR approach with a simpler set of equations gives rather comparable error norm with the CFW method which has been achieved based upon the coupled strategy with more complex eigenvector structures in particular for the refined meshes. In terms of CPU time with a strong interaction parameter, , the WCAR method takes only 0.0126 while the running time for the CFW method is 0.168 s for the same CFL number equal to 0.95.

### 1D dam-break test case

This test case is based on the experiments performed on a flume designed by Spinewine & Zech (2007). The bed material was gravel sand with dm = 1.82 mm, density ρs = 2,683 kg/m3 and Manning's coefficient of n = 0.0165. Table 6 summarizes the initial conditions for three different dam-break test cases.

Table 6

Initial conditions for 1D dam-break test case

Test
0.25 0.1 0.1 0.0
0.35 0.0 0.0 0.0
0.25 0.0 0.1 0.0
Test
0.25 0.1 0.1 0.0
0.35 0.0 0.0 0.0
0.25 0.0 0.1 0.0

This first test case (A) investigates various types of waves that are produced by a dam break in the presence of wet bed, and assesses the performance of the numerical method using Meyer-Peter and Müller and Smart sediment transport formulations. The numerical results are obtained taking Δx = 0.015 and CFL = 0.95. Figure 7 compares the experimental data and bedload and free-surface simulation results computed by the CJM method (Murillo & García-Navarro 2010) and WCAR approaches with Meyer-Peter and Müller and Smart sediment discharge formulae. Over time, the process of changing the bed profile leads to rarefaction and shock waves along with an intermediate region with a constant level. As can be seen in, the amount of erosion material simulated by WCAR based on the Smart sediment transport equation is very close to experimental values, and shock and rarefaction waves are accurately captured by the numerical model. Furthermore, the computed water level data by the Smart equation is compatible with experimental data. Negligible differences observed in the shock wave are related to the fast transient energy changes associated with the hydraulic jump.

Figure 7

Comparisons between the experimental data, free-surface and bed variations results computed by WCAR and CJM approaches with Meyer-Peter and Smart sediment discharge formulae for the test case (A) at times (a) t = 0.25 s, (b) t = 0.5 s, (c) t = 0.75 s, (d) t = 1 s, (e) t = 1.25 s and (f) t = 1.5 s.

Figure 7

Comparisons between the experimental data, free-surface and bed variations results computed by WCAR and CJM approaches with Meyer-Peter and Smart sediment discharge formulae for the test case (A) at times (a) t = 0.25 s, (b) t = 0.5 s, (c) t = 0.75 s, (d) t = 1 s, (e) t = 1.25 s and (f) t = 1.5 s.

Close modal

The test cases (B) represent a situation where morphological variations have occurred over a dry-state with the bedload height equal to zero. This test case is used to examine the numerical suitability of the defined WCAR in dealing with propagation over dry-state mobile topology using both Smart and Meyer-Peter and Müller sediment discharge formulations. Figure 8 depicts the dam-break results along with the bed sediment variations obtained based upon the WCAR against the experimental data and the weakly coupled method (WCM) presented in Juez et al. (2014). As can be observed, only a rarefaction wave moves over the dry area which contains a small drop at the location of the dam. Additionally, for the dry states test cases both Smart and Meyer-Peter and Müller sediment discharge formulations give rather identical results which might be due to the existence of the rarefaction wave over dry-state. The results are also in very good agreement with the experimental data as well as the WCM method in particular for the later times, confirming that the WCAR method is capable of simulating the rarefaction propagation over a dry state.

Figure 8

Comparisons between the experimental data, free-surface and bed variations results computed by WCAR and WCM approaches with Meyer-Peter and Smart sediment discharge formulae for the test case (B) at times (a) t = 0.25 s, (b) t = 0.5 s, (c) t = 0.75 s, (d) t = 1 s, (e) t = 1.25 s and (f) t = 1.5 s.

Figure 8

Comparisons between the experimental data, free-surface and bed variations results computed by WCAR and WCM approaches with Meyer-Peter and Smart sediment discharge formulae for the test case (B) at times (a) t = 0.25 s, (b) t = 0.5 s, (c) t = 0.75 s, (d) t = 1 s, (e) t = 1.25 s and (f) t = 1.5 s.

Close modal

The final test case (C) describes again the propagation of dam break over a dry state where the bedload itself comprised a Riemann problem at the place of dam failure and follows the same initial condition as Test A. Figure 9 illustrates the comparison between the WCAR results with the Smart and Meyer-Peter and Müller sediment discharge formulae against the experimental data and coupled Roe solver results (CRS) given in Martínez-Aranda et al. (2019). As is evident for all the calculated times, a small shock propagates to the right state, whilst the refection wave moved to the left. Moreover, the WCAR results with the Smart sediment discharge formula are in qualitative agreement with both experimental data and the CRS approach in particular for later times. In terms of the Meyer-Peter and Müller sediment discharge, although it is capable to model small front shock and left-going rarefaction, it fails to follow experimental measurements at the initial location of the dam face.

Figure 9

Comparisons between the experimental data, free-surface and bed variations results computed by WCAR and CRS approaches with Meyer-Peter and Smart sediment discharge formulae for the test case (c) at times (a) t = 0.5 s, (b) t = 1 s, (c) t = 1.5 s.

Figure 9

Comparisons between the experimental data, free-surface and bed variations results computed by WCAR and CRS approaches with Meyer-Peter and Smart sediment discharge formulae for the test case (c) at times (a) t = 0.5 s, (b) t = 1 s, (c) t = 1.5 s.

Close modal

### Dam failure due to over-topping flow

The final test case investigated here is the problem of dam failure over dry bed caused by over-topping flow which is adopted from the study by Tingsanchali & Chinnarasri (2001). The height and width of the dam crest are fixed at 0.8 and 0.3 m, respectively, and the upstream and downstream slopes of the dam are set to 1:3 and 1:2.5, respectively. The dam consists of sand with dm = 1.13, d30 = 0.52, d50 = 0.86 and d90 = 3.8. The input discharge is set to 1.23 ls−1.

The measured experimental data indicate severe erosion on the crest of the dam which is due to the rapid flow and strong reactions between the water flow and dam materials. As time passes, dam materials quickly accumulate by the flow downward of the dam. In this case study, the Meyer-Peter and Smart sediment transport formulas are used to determine the trend of dam erosion caused by over-topping flow. Computed bed levels by the WCAR approach using Mayer-Peter and Müller and Smart sediment discharge formulae at t = 30 and 60 seconds are displayed alongside experimental data in Figure 10. The computational domain is discretized to 0.05 m cells.

Figure 10

Dam failure due to over-topping flow, comparisons between the experimental data and the bed profile simulation results computed by the WCAR approach with Meyer-Peter and Smart sediment discharge formula at (a) t = 30 s and (b) t = 60 s. (c) Evolution in time of the measured water level and computed water level using Meyer-Peter and Smart sediment discharge formulae. (d) Measured over-topping discharge in time and computed over-topping discharge with Meyer-Peter and Smart sediment discharge formulae.

Figure 10

Dam failure due to over-topping flow, comparisons between the experimental data and the bed profile simulation results computed by the WCAR approach with Meyer-Peter and Smart sediment discharge formula at (a) t = 30 s and (b) t = 60 s. (c) Evolution in time of the measured water level and computed water level using Meyer-Peter and Smart sediment discharge formulae. (d) Measured over-topping discharge in time and computed over-topping discharge with Meyer-Peter and Smart sediment discharge formulae.

Close modal

As can be observed, the use of Mayer-Peter and Müller sediment transport formula is not suitable in this test case, since the calculated erosion rate is lower than measured erosion. This is due to the lack of a parameter related to the slope effect in the Meyer-Peter and Müller equation. However, as can be seen, the proposed WCAR approach with the Smart sediment discharge formula provides results comparable to the experimental data. In Figure 10(c) the relevant water reservoir level was depicted which demonstrates a good agreement with the measured data in particular for the Smart formula. Figure 10(d) exhibits the over-topping discharge computed by the WACA method along with the measured data after 120 s just upstream of the breach. As is evident, the prediction of discharge at the location of dam crest using the WCAR method with the choice of Smart formula is far from the experimental measurement. This inaccurate discharge can be rectified using non-equilibrium solid transport formulations.

In this paper, a second-order finite-volume method was proposed for solving a one-dimensional morphodynamic system using a weakly coupled scheme based on an augmented Riemann solver (WCAR). The method used an augmented Riemann solver for the hydrodynamic phase and utilized the wave decomposition and flux-wave formula for the morphodynamic estimation. The main advantages of the method is that it employed the bed variations as a fourth wave propagating from each Riemann cell interface for the hydrodynamic phase using an augmented Riemann solver. This bed variation wave was also updated through the wave decomposition and flux-wave formula for the Exner equation. Additionally, a wave decomposition and flux-wave formula significantly improved the effectiveness of the method for modeling strong interaction regimes. Validation of the proposed method was performed by comparing numerical results with the exact or analytical solution for four different test cases.

First, parabolic sediment layer transportation was simulated. The results in Figure 1 show very good agreement with the exact solution in the entire sediment layer. It was also observed that applying second-order correction expressions to the corresponding equations increased the simulation accuracy (Figure 2). Then, sediment hump propagation was modeled with the initial subcritical flow. It was observed that the WCAR method is able to simulate sediment transport accurately even with strong interaction between water flow and sediment. Next, one-dimensional dam break on moving bed using Meyer-Peter and Müller and Smart sediment discharge formulae was investigated for both wet and dry states. According to the obtained results, the numerical data produced by the WCAR method have very good agreement with the experimental results. Additionally, this method is able to accurately simulate any type of waves generated in the case of a dam break in particular for the dry state. Finally, a dam failure due to over-topping was studied. The results indicate high precision modeling in cases of the strong interaction between flow and sediment over the dry state.

Bale
D. S.
,
Leveque
R. J.
,
Mitran
S.
&
Rossmanith
J. A.
2002
A wave propagation method for conservation laws and balance laws with spatially varying flux functions
.
SIAM Journal on Scientific Computing
24
(
3
),
955
978
.
Camenen
B.
&
Larson
M.
2005
A general formula for non-cohesive bed load sediment transport
.
Journal of Estuarine, Coastal and Shelf Science
63
(
1–2
),
249
260
.
Canestrelli
A.
,
Dumbser
M.
,
Siviglia
A.
&
Toro
E. F.
2010
Well-balanced high-order centered schemes on unstructured meshes for shallow water equations with fixed and mobile bed
.
Advances in Water Resources
33
(
3
),
291
303
.
Castro Díaz
M. J.
,
Fernández-Nieto
E. D.
&
Ferreiro
A. M.
2008
Sediment transport models in Shallow Water equations and numerical approach by high order finite volume methods
.
Computers & Fluids
37
(
3
),
299
316
.
Castro Díaz
M. J.
,
Fernández-Nieto
E. D.
,
Ferreiro
A. M.
&
Parés
C.
2009
Two-dimensional sediment transport models in shallow water equations. A second order finite volume approach on unstructured meshes
.
Computer Methods in Applied Mechanics and Engineering
198
(
33–36
),
2520
2538
.
Cordier
S.
,
Le
M. H.
&
Morales de Luna
T.
2011
Bedload transport in shallow water models: why splitting (may) fail, how hyperbolicity (can) help
.
Advances in Water Resources
34
(
8
),
980
989
.
Delis
A. I.
&
Papoglou
I.
2008
Relaxation approximation to bed-load sediment transport
.
Journal of Computational and Applied Mathematics
213
(
2
),
521
546
.
Einfeldt
B.
1988
On Godunov-type methods for gas dynamics
.
SIAM Journal on Numerical Analysis
25
(
2
),
294
318
.
Exner
F. M.
1925
Über die Wechselwirkung zwischen Wasser und Geschiebe in Flüssen
.
Akademie der Wissenschaften
134
(
2a
),
165
204
.
Fernandez Luque
R.
&
Van Beek
R.
1976
Erosion and transport of bed-load sediment
.
Journal of Hydraulic Research
14
(
2
),
127
144
.
Franzini
F.
&
Soares-Frazão
S.
2018
Coupled finite-volume scheme with adapted Augmented Roe scheme for simulating morphological evolution in arbitrary cross-sections
.
Journal of Hydroinformatics
20
(
5
),
1111
1130
.
Goutière
L.
,
Soares-Frazão
S.
,
Savary
C.
,
Laraichi
T.
&
Zech
Y.
2008
One-dimensional model for transient flows involving bed-load sediment transport and changes in flow regimes
.
Journal of Hydraulic Engineering
134
(
6
),
726
735
.
Hudson
J.
2001
Numerical Techniques for Morphodynamic Modelling
.
PhD Thesis
,
Department of Mathematics, University of Reading
.
Hudson
J.
&
Sweby
P. K.
2005
A high-resolution scheme for the equations governing 2D bed-load sediment transport
.
International Journal for Numerical Methods in Fluids
47
(
10–11
),
1085
1091
.
Juez
C.
,
Murillo
J.
&
García-Navarro
P.
2013
Numerical assessment of bed-load discharge formulations for transient flow in 1D and 2D situations
.
Journal of Hydroinformatics
15
(
4
),
1234
1257
.
Juez
C.
,
Murillo
J.
&
García-Navarro
P.
2014
A 2D weakly-coupled and efficient numerical model for transient shallow flow and movable bed
.
Advances in Water Resources
71
,
93
109
.
LeVeque
R. J.
2002
Finite Volume Methods for Hyperbolic Problems
.
Cambridge University Press
,
Cambridge
.
LeVeque
R. J.
&
Pelanti
M.
2001
A class of approximate Riemann solvers and their relation to relaxation schemes
.
Journal of Computational Physics
172
(
2
),
572
591
.
Mahdizadeh
H.
2010
Modelling of Flood Waves Based on Wave Propa-Gation Algorithms with bed Efflux and Influx Including A Coupled Pipenetwork Solver
.
University of Manchester
,
Manchester
,
UK
.
Mahdizadeh
H.
&
Sharifi
S.
2019
A fully-coupled bedload sediment transport model based on a two-dimensional modified wave propagation algorithm
.
Journal of Hydraulic Research
.
doi:10.1080/00221686.2018.1562998
.
Mahdizadeh
H.
,
Stansby
P. K.
&
Rogers
B. D.
2011
On the approximation of local efflux/influx bed discharge in the shallow water equations based on a wave propagation algorithm
.
International Journal for Numerical Methods in Fluids
66
,
1295
1314
.
Mahdizadeh
H.
,
Stansby
P. K.
&
Rogers
B. D.
2012
Flood wave modeling based on a two-dimensional modified wave propagation algorithm coupled to a full-pipe network solver
.
Journal of Hydraulic Engineering
138
(
3
),
247
259
.
Mahdizadeh
H.
,
Sharifi
S.
&
Omidvar
P.
2018
On the approximation of two-dimensional transient pipe flow using a modified wave propagation algorithm
.
Journal of Fluids Engineering
140
(
7
),
071402
(7 pages)
.
Martínez-Aranda
S.
,
Murillo
J.
&
García-Navarro
P.
2019
A 1D numerical model for the simulation of unsteady and highly erosive flows in rivers
.
Computers & Fluids
181
,
8
34
.
Meyer-Peter
E.
&
Müller
R.
1984
Formulas for bed-load transport
.
Report on Second Meeting of IAHR
,
Stockholm
.
Murillo
J.
&
García-Navarro
P.
2010
An Exner-based coupled model for two-dimensional transient flow over erodible bed
.
Journal of Computational Physics
229
(
23
),
8704
8732
.
Nielsen
P.
1992
Coastal Bottom Boundary Layers and Sediment Transport
.
World Scientific
,
Singapore
.
Roe
P. L.
1981
Approximate Riemann solvers, parameter vectors, and difference-schemes
.
Journal of Computational Physics
43
(
2
),
357
372
.
Rosatti
G.
&
Fraccarollo
L.
2006
A well-balanced approach for flows over mobile-bed with high sediment-transport
.
Journal of Computational Physics
220
(
1
),
312
338
.
Saiedi
S.
1997
Coupled modeling of alluvial flows
.
Journal of Hydraulic Engineering
123
(
5
),
440
446
.
Serrano-Pacheco
A.
,
Murillo
J.
&
Garca-Navarro
P.
2012
Finite volumes for 2D shallow-water flow with bed-load transport on unstructured grids
.
Journal of Hydraulic Research
50
(
2
),
154
163
.
Smart
G. M.
1984
Sediment transport formula for steep channels
.
Journal of Hydraulic Engineering
110
(
3
),
267
276
.
Soares-Frazão
S.
&
Zech
Y.
2011
HLLC scheme with novel wave-speed estimators appropriate for two-dimensional shallow-water flow on erodible bed
.
International Journal for Numerical Methods in Fluids
66
(
8
),
1019
1036
.
Spinewine
B.
&
Zech
Y.
2007
Small-scale laboratory dam-break waves on movable beds
.
Journal of Hydraulic Research
45
(
Suppl. 1
),
73
86
.
Tingsanchali
T.
&
Chinnarasri
C.
2001
Numerical modelling of dam failure due to flow overtopping
.
Hydrological Sciences Journal
46
(
1
),
113
130
.
Wu
W.
2008
Computational River Dynamics
.
Taylor & Francis
,
London
.
Wu
W.
&
Wang
S. S. Y.
2008
One-dimensional explicit finite-volume model for sediment transport
.
Journal of Hydraulic Research
46
(
1
),
87
98
.
Wu
W.
,
Vieira
D. A.
&
Wang
S. S.
2004a
One-dimensional numerical model for nonuniform sediment transport under unsteady flows in channel networks
.
Journal of Hydraulic Engineering
130
(
9
),
914
923
.
Wu
W.
,
Vieira
D.
&
Wang
S.
2004b
One-dimensional numerical model for nonuniform sediment transport under unsteady flows in channel networks
.
Journal of Hydraulic Engineering
130
(
9
),
914
923
.