Numerical simulation of bedload sediment transport with the ability to model wet/dry interfaces using an augmented Riemann solver

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.

Different empirical formulae to evaluate are presented in Table 1.

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⁠.

In Table 2, d₃₀ and d₉₀ are the grain diameters for which 30% and 90% of the weight of a non-uniform sample is finer, respectively, and S₀ is the bed slope.

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

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.

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.

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.

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.

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.

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 ⁠.

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 m²/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/m³ 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.

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.

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.

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/m³ and Manning's coefficient of n = 0.0165. Table 6 summarizes the initial conditions for three different dam-break test cases.

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.

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.

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.

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 d_m = 1.13, d₃₀ = 0.52, d₅₀ = 0.86 and d₉₀ = 3.8. The input discharge is set to 1.23 ls⁻¹.

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.

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.