## Abstract

Catastrophic flood events worldwide have become increasingly more frequent and their dynamics mechanism has attracted much interest as how to predict them by numerical method. As the most common phenomenon occurs in the flowing process, entrainment and deposition can significantly influence flow mobility by increasing in mass and changing in flow character. In this paper, a two-dimension coupling model is presented to simulate water flow, sediment transport and bed evolution based on the shallow water assumption, depth-averaged integration as well as morphological evolution. A new term accounting for the sediment effect on the momentum conservation of water–sediment mixture is added to the model equations by assuming that the flow and the fixed bed is connected by an infinitesimally thin boundary layer, in which the erodible material gains the necessary velocity to enter the flow above. Comparison of numerical results and experimental data indicate the presented model can adequately describe the complex dynamic process, sediment transport and bed evolution. The velocity profile of flow can influence the momentum transfer between the water column and the erodible bottom boundary due to sediment exchange, further influencing flow mobility. Moreover, the velocity profile of flow changes with variations of sediment concentration, bed surface and friction resistance.

## INTRODUCTION

Water ﬂow induces sediment transport and changes in the surface morphology, which in turn, modifies the ﬂow. The dynamical mechanism of this process is an issue that has been studied for a long time (Fagherazzi & Sun 2003; Simpson & Castelltort 2006; Peng & Cao 2009; Schippa & Pavan 2009). As one of the most effective ways to study the dynamics of water flow, the mathematical model of coupled sediment transport and water flow was developed rapidly (Cao *et al.* 2004, 2006; Benkhaldoun *et al.* 2011; Li & Duffy 2011). Cao *et al.* (2004) presented a one-dimensional (1D) dam-break model to study mobile bed hydraulics of dam-break ﬂow and the induced sediment transport and morphological evolution, and provided a detailed description of dam-break ﬂuvial process. Wu & Wang (2007) proposed a similar 1D model to simulate dam-break ﬂows over mobile beds for investigating the mechanisms of morphodynamic processes. On the basis of previous study, a two-dimensional mathematical model coupling water ﬂow and sediment transport dynamics was proposed by Simpson & Castelltort (2006). This model further considered the influences of rainfall and turbulence, and could calculate the changing surface morphology through time and space. Iverson (2012) considered mass and momentum exchange between three continuous layers: the upper layer is a flow with a free surface, the middle layer is an erodible bed-sediment layer, and the lower layer represents a substrate that cannot be entrained. Hillebrand *et al.* (2017) provided a test where detailed field measurements of velocities and bed elevation changes over a three-month period from a prototype reservoir are compared with simulation results. It shows that the bed elevation changes were sensitive to the fall velocities of the finest cohesive sediment particles and sediment cohesion.

In the past decades, it has been generally accepted that reliable numerical methods for predicting the water surface, sediment transport and bed evolution can promote the study of dynamical mechanism of water flow effectively. Initially, several numerical methods were developed to solve the shallow water equations related with mass movement over a fixed bed, including finite element method (Jia & Wang 1999; Schwanenberg & Harms 2004), finite volume method (Fraccarollo & Capart 2002; Murillo & García-Navarr 2010; Benkhaldoun *et al.* 2012) and finite difference method (Jha *et al.* 1995; Wang *et al.* 2000). However, in order to exhibit the geomorphological evolution, it is necessary to couple the shallow water equations describing the flow continuity and the momentum conservation equation, the sediment continuity equation and the bathymetry evolution equation. Thus, increasing attention is being paid to methods of solving the complete coupled equations.

In this work, a new mathematical model is presented for predicting the change of water surface morphology and bed morphology in time and space through water flow and sediment transport dynamics. This model includes the shallow water equations for water flow, conservation of sediment concentration, substrate erosion and deposition as well as morphological changes. The model assumes that a thin layer which connects the flow and the fixed bed, and the sediment material is scoured and instantly mixed with the flow and transferred to the downstream as the boundary is fast evolving (Fraccarollo & Capart 2002; Papa *et al.* 2004; Lê & Pitman 2009; Iverson 2012; Iverson & Ouyang 2015). On this basis, a new term decided by entrainment and deposition is introduced to the momentum equation of the model. In order to solve the model equations more accurately, a coupled numerical method based on the finite volume methods has been applied.

This paper is organized as follows. The governing equations for coupled model of water flow, sediment transport and bed evolution are presented in the section below. The computational scheme is presented in the next section. This is followed by a section detailing the numerical results and examples, and a final section summarizes the concluding remarks.

## METHODS

### Model equations

The study of complex processes of flow over a mobile bed is mainly focused on three aspects, which are the motion of flow, the distribution of eroded sediment in flow and the evolution of the bed surface (Xia *et al.* 2010; Benkhaldoun *et al.* 2011). In order to model these processes more accurately, the shallow water equations for the mass and momentum conservations of water–sediment mixture, the mass conservation of the sediment and bed material are used in the current model. Similar coupled equations have been presented and improved by many researchers (Fagherazzi & Sun 2003; Cao *et al.* 2004; Benkhaldoun *et al.* 2012).

*et al.*2003; Papa

*et al.*2004; Iverson & Ouyang 2015) is used. Moreover, the kinetic theory of granular material is used to describe this boundary layer (James 1987; Savage 1998; Iverson 2012). During the transport process of erodible material, its velocity accelerates from a limited velocity (as the depth-averaged velocity within the boundary layer) to reach the ambient flow velocity (as the depth-averaged velocity within the moving flow) or vice versa. However, how to determine the value of the limited velocity is not fully studied. Based on the kinetic theory originally proposed by Jenkins & Savage (1983), Lê & Pitman (2009) derived the relation between limited velocity and flow velocity as

**u**

*= 0.5*

_{b}**v**when the thickness of boundary layer tends to zero, where

**u**

*is the depth-averaged velocity of boundary layer and*

_{b}**v**is the velocity at the top of the boundary layer. In their research, the value of

**v**is the same as the depth-averaged velocity of moving flow. In this paper, the value of

**v**is further studied to reflect the transfer process of erodible material more realistically. We consider that the value of

**v**relates to the velocity profile of moving flow, and further establish the relation between limited velocity and depth-averaged flow velocity as follows: where

**u**

*is the depth-averaged velocity of moving flow;*

_{ag}*σ*is a scale factor,

*σ*=

**v**/

**u**

*, and its value ranges from 0 to 1 which depends on the velocity profile of the moving flow. Some common velocity profiles investigated by Johnson*

_{ag}*et al.*(2012) are also shown in Figure 1. Considering that the erodible material attains the velocity and mixes with the flow immediately as the mass flows by, the extra momentum for the moving flow is created which will make the flow move faster (Pitman

*et al.*2003). With this assumption, the term being included to account for the sediment effect on the momentum conservation of water–sediment mixture is added to the model equations.

*et al.*2004; Simpson & Castelltort 2006; Yue

*et al.*2008; Zhang & Duan 2011) can be expressed in detail as follows: where

*t*= time;

*h*= water depth;

*u*and

*v*= depth-averaged velocity components in the

*x*and

*y*directions, respectively;

*g*= gravitational acceleration;

*c*= depth-averaged concentration of the volumetric sediment;

*ρ*= water–sediment mixture density,

*ρ*=

*ρ*(1 –

_{w}*c*) +

*ρ*, in which

_{s}c*ρ*= sediment density and

_{s}*ρ*= water density;

_{w}*ρ*= density of saturated bed,

_{b}*ρ*=

_{b}*ρ*(1 −

_{w}*α*) +

_{s}*ρ*, where

_{s}α_{s}*α*= bed sediment porosity;

_{s}*z*= bed elevation; the friction slope terms

*S*and

_{fx}*S*are written as

_{fy}*S*=

_{fx}*n*

^{2}

*u*(

*u*

^{2}+

*v*

^{2})

^{1/2}/

*h*

^{4/3}and

*S*=

_{fy}*n*

^{2}

*v*(

*u*

^{2}+

*v*

^{2})

^{1/2}/

*h*

^{4/3}in the

*x*and

*y*directions, respectively;

*n*= Manning's roughness coefﬁcient;

*D*and

*E*= deposition and entrainment rates at the interface between the water column and bed, respectively;

*u*and

_{b}*v*= limited velocity components of boundary layer. As

_{b}*σ*= 0, the system equations reduce to the common flow model. However, in this study, the suspended and bed loads are treated as a single mode following Cao

*et al.*(2004), and which provides that the sediment transport takes place in a continuum between pure suspended load and pure bed load. Equation (2) represents mass conservation for the water–sediment mixture with the term of mass exchange between the flow and the erodible bed on the right-hand side. Equations (3) and (4) are equations for momentum conservation in the

*x*and

*y*directions, respectively. The first term represents the extra momentum that erodible material has within the boundary layer. The second three terms on the right-hand side of the momentum equations represent, respectively, momentum transfer due to sediment exchange between the water column and the erodible bottom boundary, and spatial variations in sediment concentration. The last term on the right-hand side of the momentum equations represents the effects of bed topography and friction loss. Equation (5) represents sediment mass conservation, and Equation (6) links the local variation in bed elevation to sediment transfer.

*φ*= calibrated constant;

*θ*= Shields parameter,

*θ*=

*u*

_{*}

^{2}/

*sgd*;

*θ*= critical Shields parameter for initiation of sediment movement;

_{c}*d*= grain diameter of bed sediment;

*s*= submerged specific gravity of sediment,

*s*=

*ρ*/

_{s}*ρ*–1;

_{w}*u*

_{*}= friction velocity,

*u*

_{*}= (

*gh*(

*S*

_{fx}^{2}+

*S*

_{fy}^{2}))

^{1/2}. For the deposition of non-cohesive sediment the following relation is used (Cao 1999): where

*C*= near-bed volumetric sediment concentration,

_{a}*C*=

_{a}*ca*, in which

_{c}*a*= coefficient larger than unity, and can be obtained by

_{c}*a*=

_{c}*min*(2, (1 –

*α*)/

_{s}*c*);

*m*= exponent indicating the effects of hindered settling due to high sediment concentrations, and here

*m*= 2 (Cao

*et al.*2004);

*ω*= setting velocity of a single particle in tranquil water, which is expressed as (Li & Duffy 2011): where

_{o}*υ*= water kinematic viscosity.

### Numerical scheme

**U**is a vector of the conserved variables;

**F**and

**G**are the convective ﬂux vectors of the ﬂow in the

*x*and

*y*directions, respectively;

**S**and

**T**is the source term vectors of the flow in the

*x*and

*y*directions, respectively. In this paper, the hyperbolic system of conservation laws given by Equation (10) is discretized using finite volumes on regular rectangular grids (shown in Figure 2). With this grid system, two separate one-dimension problems in

*x*and

*y*directions can be obtained by dividing Equation (10) based on the operator-splitting technique (Liang

*et al.*2006) and expressed as follows:

*L*and

_{x}*L*represent the computational procedures based on predictor-corrector method and are used for solving Equations (11) and (12), respectively. Each computational procedure is both composed of Roe type scheme (Toro 2009), Monotonic Upstream-Centered Scheme for Conservation Laws (MUSCL) approach which is used for reconstructing the interface flux in order to obtain second order accuracy solution in space and avoid spurious oscillations (Zoppou & Roberts 2000) and McCormack scheme which is used for obtaining second order time accuracy (Ouyang

_{y}*et al.*2013; Liu

*et al.*2016). The solution at next time step can be obtained by operating the procedures each twice with the time step calculated by the stability criterion (Brufau

*et al.*2004). For

*L*, the explicit time-marching conservative finite volume formula of Equation (12) is: where superscript

_{x}*n*is the time level; superscript

*i*is the cell index; Δ

*x*is the distance between the centroids of the cell; Δ

*t*is the time step;

**F**

*and*

_{e}**F**

*are the fluxes through the west and east cell interfaces. Here, the internal flux, for example*

_{w}**F**

*evaluated by Roe scheme at the interface, is given as: where*

_{e}**F**

*and*

_{l}**F**

*are the interface fluxes on both sides of the east cell interface, and calculated by the left and right Riemann states*

_{r}**U**

*and*

_{l}**U**

*;*

_{r}**is the diagonal matrix;**

*Λ***R**= (

*e*

_{1},

*e*

_{2}, …,

*e*) and

_{k}*e*is the

_{k}*k*h right eigenvectors of

**J**, where

**J =**

*∂*

**F**

*/*

_{e}*∂*

**U**is the Jacobian matrix of the flux

**F**

*. Furthermore, coupling the MUSCL approach with the Roe scheme allows us to reconstruct the interface functions*

_{e}**U**

*and*

_{l}**U**

*for obtaining accuracy and avoiding spurious oscillations. The*

_{r}**U**

*and*

_{l}**U**

*functions can be expressed as: where The function M is a min-mod flux limiter and can be written as: In order to achieve second order accuracy in time, the McCormack scheme is also adopted as: where the superscripts*

_{r}*pr*and

*cr*indicate the predictor and corrector steps, respectively. The stability criterion can be expressed as: where

*cfl*= Courant number and less than 1.

## STUDY AREA AND DATA

In this study, we focus on the dynamic process of the water flow by considering the effects of sediment transport and bed evolution. For this purpose, we ran the hydrodynamic model and numerical scheme with some cases of dam break which have different experimental conditions. Numerical results are compared with some experimental data in order to verify the applicability of the presented model. The experimental data of water depth from CADAM concerted action (Hiver 2000; Liang & Marche 2009) is used to test the numerical scheme. Experimental measurements of water depth and bed evolution from the Louvain experiment (Fraccarollo & Capart 2002), the Taipei experiment (Capart & Young 1998) and the Civil and Environmental Engineering Department laboratory (Palumbo *et al.* 2008; Goutiere *et al.* 2011) are also used to test the accuracy of the model.

## RESULTS AND DISCUSSION

The presented model and the corresponding numerical scheme are used to simulate several benchmark tests and results are compared with laboratory measurements when available. For all tests, *g* = 9.8 m/s and *ρ _{w}* = 1,000 kg/m

^{3}.

### Simulation of 1D dam-break over triangular hump

In order to verify the accuracy of the current numerical scheme, an experiment of 1D dam-break over triangular hump without erosion is reproduced, which was originally performed with the CADAM concerted action (Hiver 2000; Liang & Marche 2009). The detailed setting of this experiment is shown in Figure 3, where the horizontal domain is 38 m long and the gate located 15.5 m away from the upstream end. The still water surface elevation of the reservoir is 0.75 m, and a symmetric triangular hump of 0.4 m high and 6 m long is 13 m away from the gate at the downstream. A constant Manning coefficient *n* = 0.0125 is used for clear water throughout the domain. The upstream and downstream ends of the domain are assumed to be solid wall and free outlet, respectively. Seven gauges are located at 2 m, 4 m, 8 m, 10 m, 11 m, 13 m and 20 m downstream of the gate for recording the time history of water depth to reflect the complicated process of water wave propagation. It is assumed that the dam suddenly fails and the initial still water in the reservoir rushes onto the downstream ﬂoodplain. The comparisons between simulation results and laboratory measurements of curves of the flow height versus the time are illustrated in Figure 4. It is shown that the numerical water levels agree well with the measured values except for the last gauge furthest away from the gate. A simulation with a small grid of Δ*x* = 0.01 m was carried out in order to assess the less satisfactory agreement with the experimental measurement in Gauge 7, and the numerical result is also presented in the same panel in Figure 4. However, when compared with the previous numerical result, no significant difference can be found. It should be considered that the complex and unstable situation of water flow after the hump is hard to predict with the presented model equations as a possible explanation, and more factors suited to this situation should be included.

### Comparison of computational results with values of *σ* and sediment particle diameter

As discussed above, previous models which used non-sliding boundary condition and neglected the entrainment and deposition effect on the momentum equation underestimated the dam-break flow erosion ability and mobility. Thus comparisons of the model with *σ* = 0, 0.3, 0.7 and 1 were done with the problem of 1D dam-break over an erodible bed. In this problem, the channel was set to 200 m long. The initial upstream water depth was 10 m and the downstream water depth 2 m. The channel bed was composed of non-cohesive uniform sediment and the initial height *z* = 0 m. According to Cao *et al.* (2004), the relevant parameters are shown in Table 1. Figure 5 shows the water depth variations along the channel for a dam-break flow over an erodible bed at time *t* = 4, 6 and 8 s, using different values of *σ*. At *σ* = 0, the numerical result is the same as the results calculated by Cao *et al.* (2004). It can be expected that the erosion depth, volume fraction of particle in the model increases with the *σ* increase, especially in the original dam site, which indicates that the type of velocity profile can influence the erosion ability of flow, further influencing the free surface profiles. It should be noted that the propagation distance of flow is enhanced with a larger value of *σ*, which means more extra momentum is absorbed from the boundary layer. It makes the flow move faster further than it otherwise would. A sensitivity test of sediment particle diameter was performed. As shown in Figure 6, it can be seen that the computed bed elevation changes are sensitive to the sediment particles. The finer the bed material, the greater the deviation from the original bed level. This feature is in agreement with results reported by Benkhaldoun *et al.* (2012).

Coefficients . | n
. | m
. | υ
. | α
. _{s} | φ
. | θ
. _{c} |
---|---|---|---|---|---|---|

Values | 0.03 | 2 | 1.2 × 10^{−6} | 0.4 | 0.015 | 0.045 |

Coefficients . | n
. | m
. | υ
. | α
. _{s} | φ
. | θ
. _{c} |
---|---|---|---|---|---|---|

Values | 0.03 | 2 | 1.2 × 10^{−6} | 0.4 | 0.015 | 0.045 |

### Comparison of computational results with Louvain and Taipei experiments

In this section, model flow dynamics with fully coupled sediment transport is implemented for the Louvain experiment (Fraccarollo & Capart 2002) and the Taipei experiment (Capart & Young 1998) as two laboratory-scale dam-break experiments which involve similar devices and procedures. In the Louvain experiment, the equivalent spherical diameter of sediment particles is 3.5 mm and the speciﬁc gravity is 1.54. In the Taipei experiment, sediment particles have a uniform size diameter of 6.1 mm and specific gravity of 1.048. The test reach was 2.5 m long, 35 cm high and 10 cm wide for the Louvain experiment and 1.2 m long, 70 cm high and 20 cm wide for the Taipei experiment. For the two experiments, the sluice gate is all placed in the middle of the flumes, and the water depths on the upstream and the downstream of it are 10 cm and 0 cm, respectively.

The dam-break flow wave of these experiments is released by rapidly lifting the sluice gate; however, it is difficult to determine the velocity profile during the wave propagation. Therefore, both experiments are simulated by applying different values of *σ*. The computational parameters are listed in Table 2. With *σ* = 0.1, 0.5, 1, the simulated water and bed surface elevations along the centreline of the channel are compared with the measured values, and shown in Figure 7 (Louvain) and Figure 8 (Taipei), respectively. In the Louvain experiment, the wave front locations, the hydraulic jump locations and the erosion magnitude are in fairly good agreement with the measurement in the experiment except for the erosion magnitude at *t* = 5*t _{o}*. The simulated results indicate that the value of

*σ*has a significant impact on the process of erosion and flow movement. With a larger value of

*σ*, more momentum can transfer between the water column and the erodible bottom boundary due to sediment exchange. Moreover, at different times, the value of

*σ*matching measurement is different. It can be considered that under the effect of erosion and friction, the velocity profile of wave is changing constantly during the wave propagation. Compared with measurement, the simulated erosion magnitude at

*t*= 5

*t*is small, which can be interpreted that the actual limited velocity of boundary layer is larger than the depth-averaged velocity of flow, and leads to a larger erosion magnitude. Although the prediction of bed surface elevations is not as good as the water surface, considering the small timescale and the small magnitude of the water elevation change, it can be considered satisfactory. In the Taipei experiment, the influence of

_{o}*σ*on the process of erosion and flow movement is more obvious. The mode predicted the wave front location and the erosion magnitude quite well with the measured data by using a value of

*σ*between 0 and 0.5. However, it fails to predict the locations of the hydraulic jump, and more studies are needed to verify this cause.

Coefficients . | n
. | m
. | υ
. | α
. _{s} | φ
. | θ
. _{c} |
---|---|---|---|---|---|---|

Values (Louvain) | 0.025 | 2 | 1.2 × 10^{−6} | 0.28 | 1.25 | 0.045 |

Values (Taipei) | 0.025 | 2 | 1.2 × 10^{−6} | 0.3 | 1 | 0.045 |

Coefficients . | n
. | m
. | υ
. | α
. _{s} | φ
. | θ
. _{c} |
---|---|---|---|---|---|---|

Values (Louvain) | 0.025 | 2 | 1.2 × 10^{−6} | 0.28 | 1.25 | 0.045 |

Values (Taipei) | 0.025 | 2 | 1.2 × 10^{−6} | 0.3 | 1 | 0.045 |

### Simulation of a dam-break flow in an erodible channel with a sudden enlargement

Experimental measurements for this test case were obtained in the laboratory of the Civil and Environmental Engineering Department of the Université Catholique de Louvain, Belgium (Palumbo *et al.* 2008; Goutiere *et al.* 2011). The structure of the experiment flume is 6 m long with a non-symmetrical sudden enlargement from 0.25 m to 0.5 m in width located 1 m downstream of the gate. The initial water depth is 0.25 m and 0 m in the upstream and downstream of the gate, respectively. There is an initial sediment layer of 0.1 m all over the flume. A simple sketch of this experiment is shown in Figure 9. The sediment consists of coarse sand with median diameter 1.82 mm and a density of 2,680 kg/m^{3}, deposited with a porosity *α _{s}* = 0.53. A free-slip boundary condition was applied on all side walls, with a free outﬂow boundary condition being used at the downstream outlet.

We also used different values of *σ* = 0, 0.5, 1 to perform this experiment, and a Manning *n* of 0.025 m^{−1/3} was set. Six gauges and two cross sections located downstream of the gate were set for recording the time history of water depth and bed surface to reflect the complicated process of water wave propagation and sediment transport. Figure 10 shows the model predicted and laboratory observed evolution of water surface at the measurement gauges, and the comparisons between the predicted and observed bed surface at two measured cross sections are shown in Figure 11. However, the difference of water surface with different values of *σ* can be neglected due to the calculation results, and therefore just the water surface with *σ* = 0.5 is shown. It shows that the model predicted evolution of water surfaces at the measurement gauges and cross-sectional bed surfaces are similar to the observed data. For the evolution of bed surface at *x* = 4.1 m, a comparatively obvious difference caused by *σ* is mainly focused on the position near the corner. It may be caused by the rapid and obvious change of bed surface, further changing the velocity profile of flow wave. In contrast, the change of bed surface *x* = 4.4 m is moderate, and the difference of erosion magnitude is also not obvious with different values of *σ*. It can be considered that the velocity profile of flow wave can be influenced by the bed surface, and further affects the transfer of erodible material.

## CONCLUSION

In this paper, we assumed that a thin layer connects the flow and the fixed bed, and the eroded material can instantly mix with the flow and transfer to the downstream as the boundary is fast evolving. Based on this assumption, a new term decided by entrainment and deposition was introduced to the momentum equation of the model. This term is used for considering the sediment effect on the momentum conservation of water–sediment mixture, which is related to the velocity profile of the flow. Several numerical tests were performed to verify the feasibility of the presented model to predict events such as sediment transport caused by tsunami and/or floods. Comparison of numerical results and experimental data indicates that the presented model can adequately describe the complex dynamic process, sediment transport and bed evolution. Through the analysis of simulated results according to the actual scales and some model functions, it was found that the effect of limited velocity on the process of erosion and flow movement should be considered. The value of limited velocity mainly depends on the velocity profile of flow wave. With the increase of the value of limited velocity, the erosion amount increases correspondingly. Moreover, the computed bed elevation changes are sensitive to the sediment particles. However, the velocity profile of flow wave can be changed during the flow movement, which may be related to the changes of sediment concentration, bed surface and friction resistance. The effect of the limited velocity on water level and bed equations should be further considered in a modification of the model. More extensive information about the velocity profile of flow wave and experimental data is also needed for calibrating the presented model terms.

## ACKNOWLEDGEMENTS

The authors thank two anonymous reviewers for helpful suggestions. Financial support from the NSFC (Grant No. 41272346, 41101008) and the National Science Fund for Distinguished Young Scholars (41225011) is acknowledged.

## REFERENCES

*.*