1. Introduction

Alluvial rivers often have lateral movements: meandering or channel migration. The instability of the river channel flow tends to develop a curved channel pattern, in which the flow is forced to follow the channel's curvature; the centrifugal force thus created pushes the flow toward the outer bank, and the associated superelevation of the water surface drives the flow near the bed back toward the inner bank. The balance of these two forces creates a vertical recirculation, known as the helical flow or secondary current, in a channel bend. The upper part of the helical flow (near the water surface) is toward the outer bank, and the lower part (near bed) of

the flow is toward the inner bank. Sediment transport in a curved channel is strongly affected by such a helical flow system. Since more sediment particles are distributed near the bed in a vertical profile, the helical current distributes more sediment load to the inner bank and less sediment to the outer bank. As a result, erosion would occur along the outer bank and deposition along the inner bank. Inevitably, a skewed channel cross section is developed in channel bends with a lower bed near the outer bank and higher bed near the inner bank. In turn, more and more flow would be distributed along the outer bank, causing erosion and mechanic instability of the outer bank. The fundamental theory of the curved channel fluid dynamics has been established by Rozovskii [31]. This meander migration phenomenon had been observed in the field by Hickin and Nanson [13, 14], Parker [28], Begin [2, 3], as well as in the laboratory by Friedkin [10] and Chang et al. [4], among others. When multiple sub-channels coexist, in braided rivers, each of the curved sub-channels would develop under the influence of the same mechanism.

model was tested using field data of Goodwin Creek in Mississippi with promising results. Onda et al. [25] studied bank erosion process in a curved experimental channel using a 2D depth-averaged model using non-equilibrium sediment transport method. Waterman and Garcia [40] reported the development of a bank erosion submodel for banks with two-layered soil structure: a cohesive upper layer and non-cohesive lower layer. It was found bank slope was reduced by large flow events and steepened with lower flows. Iwasaki et al. [17] studied morphodynamic process of a densely vegetated meander river during a large flood event using a 2D model. The bar formation in the river was found contributed strongly to meander

Modeling River Morphodynamic Process Using a Depth-Averaged Computational Model…

To simulate the channel migration process using a depth-averaged 2D model, the

In this paper, a bank erosion model is developed based on a general hydrodynamic

CCHE2D is a depth-integrated 2D model for simulating free-surface turbulent flows, sediment transport, and morphological change. This is a finite element-based model of the collocation method using quadrilateral mesh ([19, 20]). The governing equations solving the flow are two-dimensional depth-integrated Reynolds equa-

> ∂hτxx ∂x þ ∂hτxy ∂y

∂hτyx ∂x þ ∂hτyy ∂y

where u and v are depth-integrated velocity components in x and y directions, respectively; t is the time; g is the gravitational acceleration; η is the water surface

<sup>þ</sup> τη<sup>x</sup> � <sup>τ</sup>bx

<sup>þ</sup> τη<sup>y</sup> � <sup>τ</sup>by

<sup>ρ</sup><sup>h</sup> <sup>þ</sup> <sup>f</sup> Cor<sup>v</sup> (1)

<sup>ρ</sup><sup>h</sup> � <sup>f</sup> Cor<sup>u</sup> (2)

model should be capable of capturing the following mechanisms in addition to general sediment transport: (1) the effect of the helical motion on the sediment transport in meandering channels, (2) bank erosion including mass failure, and (3) the moving boundary problem due to bank retreat. Nagata et al. [24], Duan et al. [6], and Jia et al. [21] developed 2D channel meandering models that adopted

and sediment transport model, CCHE2D ([19, 20]). Bank surface erosion, basal erosion, and mass failure are simulated based on the approaches of Osman and Thorne [26, 27] and Hanson and Simon [14]. The secondary helical current effects on suspended sediment and bed-load sediment transport have been considered. Since this is a two-dimensional model, computational mesh has to be adjusted when the bank boundaries move due to erosion. The processes of flow, sediment transport, bed change, and bank erosion are simulated on a mesh at each time step. After the bank lines have been moved by erosion, a new mesh conforming to the new bank lines is created, and the flow field and bed topography are interpolated from the current mesh to the new one. The computations of flow, sediment transport, bed change, and bank erosion are then continued on the new mesh for the next time step. Numerical tests using data of fixed bank experiments are conducted to validate the secondary current effect. Bank erosion capabilities are tested using hypothetical cases, and the

model has been applied to a field case of Chuoshui River in Taiwan.

∂η ∂x þ 1 h

∂η ∂y þ 1 h

development and thus back erosion [12].

DOI: http://dx.doi.org/10.5772/intechopen.86692

the moving grid technique.

2. Materials and methods

∂u ∂t þ u ∂u ∂x þ v ∂u <sup>∂</sup><sup>y</sup> ¼ �<sup>g</sup>

∂v ∂t þ u ∂v ∂x þ v ∂v <sup>∂</sup><sup>y</sup> ¼ �<sup>g</sup>

33

2.1 Hydrodynamic, sediment transport model

tions in the Cartesian coordinate system:

The transversal component of secondary flow velocity near the bed is always toward the center of curvature, and it deviates from the longitudinal direction of the total velocity near the bed. Empirical functions to describe the transversal component of secondary flow velocity have been formulated based on experimental data [7, 9, 22]. To simulate bank erosion, additional processes have to be considered. The bank erosion process is generally more complicated and has two categories ([38, 39]): basal erosion and geotechnical bank failure. The former (also called toe erosion) is caused by shear stress of the flow constantly eroding the base of the channel bank. When the basal erosion takes too much material away from the toe, a bank soil mechanic failure will take place. Basal erosion is a general process for both cohesive and non-cohesive banks.

Because river morphodynamics involves multiple processes such as turbulent flow, channel bed change, bank erosion, and sediment transport, numerical models can be used to handle most of the processes and associated parameters effectively. With the rapid development of computer technology and facilities, bank erosion has been studied with numerical simulations. Struiksma et al. [36], Shimizu and Ikekura [32], Jin and Steffler [18], Jia and Wang [19], and Wu and Wang [43] have developed the depth-averaged 2D models that considered the effect of helical flow. Finnie et al. [8] added the secondary flow effect to a depth-averaged model by solving a transport equation for stream-wise vorticity. Lien et al. [23] included the dispersion stresses due to integration into a depth-integrated model. Fang et al. [8] considered the influence of the helical current-induced vertical velocity on suspended sediment distribution and improved the calculation results. Simon et al. [33] proposed a sophisticated bank stability and toe erosion model, which considered wedge-shaped bank failure with several distinct bank material layers and irregular bank geometry. Their model is able to incorporate root reinforcement and surcharge effects of six vegetation species, including willows, grasses, and large trees, and can simulate saturated and unsaturated soil strength considering the effect of pore water pressure. Abdul-Kadir and Ariffin [1] summarized bank erosion capabilities of 25 numerical models including 1D, 2D, and 3D methods. Most models can simulate either cohesive or non-cohesive banks, a few is capable of layered banks, and no attempt has been found for heterogeneous banks.

Recently more influential factors are included in numerical models for better predictions. Xiao et al. [44] and Gholami and Khaleghi [11] studied bank erosion affected by instream vegetation. Rinaldi et al. [30] included groundwater effect. Lai et al. [22] coupled a 2D hydrodynamic and sediment transport with a soil mechanicbased and multilayered bank stability model of Simon et al. [33, 35]. The coupled

#### Modeling River Morphodynamic Process Using a Depth-Averaged Computational Model… DOI: http://dx.doi.org/10.5772/intechopen.86692

model was tested using field data of Goodwin Creek in Mississippi with promising results. Onda et al. [25] studied bank erosion process in a curved experimental channel using a 2D depth-averaged model using non-equilibrium sediment transport method. Waterman and Garcia [40] reported the development of a bank erosion submodel for banks with two-layered soil structure: a cohesive upper layer and non-cohesive lower layer. It was found bank slope was reduced by large flow events and steepened with lower flows. Iwasaki et al. [17] studied morphodynamic process of a densely vegetated meander river during a large flood event using a 2D model. The bar formation in the river was found contributed strongly to meander development and thus back erosion [12].

To simulate the channel migration process using a depth-averaged 2D model, the model should be capable of capturing the following mechanisms in addition to general sediment transport: (1) the effect of the helical motion on the sediment transport in meandering channels, (2) bank erosion including mass failure, and (3) the moving boundary problem due to bank retreat. Nagata et al. [24], Duan et al. [6], and Jia et al. [21] developed 2D channel meandering models that adopted the moving grid technique.

In this paper, a bank erosion model is developed based on a general hydrodynamic and sediment transport model, CCHE2D ([19, 20]). Bank surface erosion, basal erosion, and mass failure are simulated based on the approaches of Osman and Thorne [26, 27] and Hanson and Simon [14]. The secondary helical current effects on suspended sediment and bed-load sediment transport have been considered. Since this is a two-dimensional model, computational mesh has to be adjusted when the bank boundaries move due to erosion. The processes of flow, sediment transport, bed change, and bank erosion are simulated on a mesh at each time step. After the bank lines have been moved by erosion, a new mesh conforming to the new bank lines is created, and the flow field and bed topography are interpolated from the current mesh to the new one. The computations of flow, sediment transport, bed change, and bank erosion are then continued on the new mesh for the next time step. Numerical tests using data of fixed bank experiments are conducted to validate the secondary current effect. Bank erosion capabilities are tested using hypothetical cases, and the model has been applied to a field case of Chuoshui River in Taiwan.

### 2. Materials and methods

#### 2.1 Hydrodynamic, sediment transport model

CCHE2D is a depth-integrated 2D model for simulating free-surface turbulent flows, sediment transport, and morphological change. This is a finite element-based model of the collocation method using quadrilateral mesh ([19, 20]). The governing equations solving the flow are two-dimensional depth-integrated Reynolds equations in the Cartesian coordinate system:

$$\frac{\partial u}{\partial t} + u \frac{\partial u}{\partial \mathbf{x}} + v \frac{\partial u}{\partial \mathbf{y}} = -\mathbf{g} \frac{\partial \eta}{\partial \mathbf{x}} + \frac{\mathbf{1}}{h} \left( \frac{\partial h \tau\_{\mathbf{x}\mathbf{x}}}{\partial \mathbf{x}} + \frac{\partial h \tau\_{\mathbf{xy}}}{\partial \mathbf{y}} \right) + \frac{\tau\_{\mathbf{yx}} - \tau\_{bx}}{\rho h} + f\_{\text{Car}} v \tag{1}$$

$$\frac{\partial v}{\partial t} + u \frac{\partial v}{\partial \mathbf{x}} + v \frac{\partial v}{\partial \mathbf{y}} = -\mathbf{g} \frac{\partial \eta}{\partial \mathbf{y}} + \frac{\mathbf{1}}{h} \left( \frac{\partial h \tau\_{\text{yx}}}{\partial \mathbf{x}} + \frac{\partial h \tau\_{\text{yy}}}{\partial \mathbf{y}} \right) + \frac{\tau\_{\text{yy}} - \tau\_{\text{by}}}{\rho h} - f\_{\text{Cor}} u \tag{2}$$

where u and v are depth-integrated velocity components in x and y directions, respectively; t is the time; g is the gravitational acceleration; η is the water surface

the flow is toward the inner bank. Sediment transport in a curved channel is strongly affected by such a helical flow system. Since more sediment particles are distributed near the bed in a vertical profile, the helical current distributes more sediment load to the inner bank and less sediment to the outer bank. As a result, erosion would occur along the outer bank and deposition along the inner bank. Inevitably, a skewed channel cross section is developed in channel bends with a lower bed near the outer bank and higher bed near the inner bank. In turn, more and more flow would be distributed along the outer bank, causing erosion and mechanic instability of the outer bank. The fundamental theory of the curved channel fluid dynamics has been established by Rozovskii [31]. This meander migration phenomenon had been observed in the field by Hickin and Nanson [13, 14], Parker [28], Begin [2, 3], as well as in the laboratory by Friedkin [10] and Chang et al. [4], among others. When multiple sub-channels coexist, in braided rivers, each of the curved sub-channels would develop under the influence of the

Current Practice in Fluvial Geomorphology - Dynamics and Diversity

The transversal component of secondary flow velocity near the bed is always toward the center of curvature, and it deviates from the longitudinal direction of the total velocity near the bed. Empirical functions to describe the transversal component of secondary flow velocity have been formulated based on experimental data [7, 9, 22]. To simulate bank erosion, additional processes have to be considered. The bank erosion process is generally more complicated and has two categories ([38, 39]): basal erosion and geotechnical bank failure. The former (also called toe erosion) is caused by shear stress of the flow constantly eroding the base of the channel bank. When the basal erosion takes too much material away from the toe, a bank soil mechanic failure will take place. Basal erosion is a general process for both

Because river morphodynamics involves multiple processes such as turbulent flow, channel bed change, bank erosion, and sediment transport, numerical models can be used to handle most of the processes and associated parameters effectively. With the rapid development of computer technology and facilities, bank erosion has been studied with numerical simulations. Struiksma et al. [36], Shimizu and Ikekura [32], Jin and Steffler [18], Jia and Wang [19], and Wu and Wang [43] have developed the depth-averaged 2D models that considered the effect of helical flow. Finnie et al. [8] added the secondary flow effect to a depth-averaged model by solving a transport equation for stream-wise vorticity. Lien et al. [23] included the dispersion stresses due to integration into a depth-integrated model. Fang et al. [8]

considered the influence of the helical current-induced vertical velocity on

banks, and no attempt has been found for heterogeneous banks.

suspended sediment distribution and improved the calculation results. Simon et al. [33] proposed a sophisticated bank stability and toe erosion model, which considered wedge-shaped bank failure with several distinct bank material layers and irregular bank geometry. Their model is able to incorporate root reinforcement and surcharge effects of six vegetation species, including willows, grasses, and large trees, and can simulate saturated and unsaturated soil strength considering the effect of pore water pressure. Abdul-Kadir and Ariffin [1] summarized bank erosion capabilities of 25 numerical models including 1D, 2D, and 3D methods. Most models can simulate either cohesive or non-cohesive banks, a few is capable of layered

Recently more influential factors are included in numerical models for better predictions. Xiao et al. [44] and Gholami and Khaleghi [11] studied bank erosion affected by instream vegetation. Rinaldi et al. [30] included groundwater effect. Lai et al. [22] coupled a 2D hydrodynamic and sediment transport with a soil mechanicbased and multilayered bank stability model of Simon et al. [33, 35]. The coupled

same mechanism.

32

cohesive and non-cohesive banks.

elevation; ρ is the density of water; h is the local water depth; f Cor is the Coriolis parameter; τxx, τxy, τyx, and τyyare depth-integrated Reynolds stresses; and τηx, τηy, τbx, and τby are shear stresses on the water surface and the bed. Free-surface elevation of the flow is calculated by the depth-integrated continuity equation:

$$\frac{\partial h}{\partial t} + \frac{\partial uh}{\partial x} + \frac{\partial vh}{\partial y} = 0 \tag{3}$$

where φbk ¼ qb <sup>∗</sup> <sup>k</sup>= pbk

phk <sup>¼</sup> <sup>∑</sup><sup>N</sup>

where

Figure 1.

35

secondary current on suspended sediment.

coefficient for channel bed, <sup>n</sup><sup>0</sup> <sup>¼</sup> <sup>d</sup>1=<sup>6</sup>

DOI: http://dx.doi.org/10.5772/intechopen.86692

<sup>j</sup>¼<sup>1</sup>pbjdj<sup>=</sup> dk <sup>þ</sup> dj

2.2. Secondary current effect

sediment particle due to the bed slope:

determined by <sup>τ</sup>ck <sup>¼</sup> <sup>0</sup>:<sup>03</sup> <sup>γ</sup>ð Þ <sup>s</sup> � <sup>γ</sup> dk phk=pek � �0:<sup>6</sup>

� � and pek <sup>¼</sup> <sup>∑</sup><sup>N</sup>

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi <sup>γ</sup>ð Þ <sup>s</sup>=<sup>γ</sup> � <sup>1</sup> gd<sup>3</sup>

� � q

k

Modeling River Morphodynamic Process Using a Depth-Averaged Computational Model…

capacity, qb <sup>∗</sup> <sup>k</sup> is the equilibrium transport rate of the kth size class of bed load per unit width (kg/m/s), pbk is the bed material gradation, n is the Manning's roughness

<sup>j</sup>¼<sup>1</sup> pbjdk<sup>=</sup> dk <sup>þ</sup> dj

� �.

to the grain roughness, τ<sup>b</sup> is the bed shear stress, τck is the critical shear stress

In curved open channels, the flow is forced to follow a curved path with a variable radius of curvature (Figure 1a). On a bed with a transversal slope (Figure 1b), the bed-load motion is different from that with a stream-wise slope only. The path of a near-bed sediment particle is affected by main flow shear, stream-wise slope, as well as by the gravity component on the transversal direction. Van Bendegom's formula ([37]) was applied to calculate the moving angle of the

tan <sup>ϕ</sup> <sup>¼</sup> sin <sup>α</sup> � <sup>1</sup>

Suspended load and bed-load motion affected by the secondary flow and the gravity. (a) Definition of longitudinal and secondary current velocities. (b) Effect of transverse bed slope and secondary flow. (c) Effect of

G ∂ζ ∂y

cos <sup>α</sup> � <sup>1</sup> G ∂ζ ∂x

and exposure probabilities for the kth size class of bed material, defined as

is a nondimensional bed-load transport

, and phk and pek are the hiding

(11)

<sup>50</sup> =20 is the Manning's coefficient corresponding

Turbulence eddy viscosity is computed with the depth-integrated mixing length eddy viscosity model:

$$v\_t = \overline{l}^2 \sqrt{2\left(\frac{\partial u}{\partial \mathbf{x}}\right)^2 + 2\left(\frac{\partial v}{\partial \mathbf{y}}\right)^2 + \left(\frac{\partial u}{\partial \mathbf{y}} + \frac{\partial v}{\partial \mathbf{x}}\right)^2 + \left(\overline{\frac{\partial U}{\partial \mathbf{z}}}\right)^2} \tag{4}$$

$$\tilde{I} = \frac{1}{h} \left[ \kappa z \sqrt{\left( 1 - \frac{z}{h} \right)} dz \approx 0.267 \kappa h \right. \tag{5}$$

$$\frac{\overline{\partial U}}{\partial \mathbf{z}} = \frac{1}{h} \int \frac{\partial U}{\partial \mathbf{z}} d\mathbf{z} = \mathbf{C}\_m \frac{u\_\*}{h\kappa} \tag{6}$$

where u<sup>∗</sup> is the shear velocity, κ = 0.41 is the Karman constant, and Cm ≈ 2.34375 is based on the vertical log distribution of flow velocity ([19]).

Nonuniform suspended and bed-load sediment transport can be simulated. The depth-integrated convection-diffusion equation is solved for the suspended sediment transport:

$$\frac{\partial hc}{\partial t} + \frac{\partial uhc}{\partial \mathbf{x}} + \frac{\partial vhc}{\partial \mathbf{y}} - \frac{\partial}{\partial \mathbf{x}} \left[ \varepsilon\_i h \frac{\partial c}{\partial \mathbf{x}} \right] - \frac{\partial}{\partial \mathbf{y}} \left[ \varepsilon\_i h \frac{\partial c}{\partial \mathbf{y}} \right] = a a \rho\_i (\mathbf{c}\_\* - \mathbf{c}) - \mathbf{S}\_r \tag{7}$$

where c is the depth-integrated sediment concentration. The diffusivity coefficient for suspended sediment ε<sup>s</sup> ¼ νt=σ<sup>c</sup> with the Schmidt number 0:5≤σ<sup>c</sup> ≤1. c\* and ω<sup>s</sup> are the sediment transport capacity and settling velocity, and α is a coefficient. The source term Sr represents the dispersion due to the vertical distribution of flow velocity and suspended sediment concentration. Bed load is computed with the mass conservation equation:

$$\frac{\partial(\delta c\_b)}{\partial t} + \frac{\partial q\_{bx}}{\partial x} + \frac{\partial q\_{by}}{\partial y} + \frac{1}{L} \left(q\_b - q\_{\ast b}\right) + \mathcal{S}\_{bank} = \mathbf{0} \tag{8}$$

where cb and qb denote bed-load concentration and transport rate and "\*" denotes capacity. δ is the bed-load layer thickness. The subscripts "bx" and "by" indicate the component of bed load in x and y directions. L is the adaptation length of the bed load representing the non-equilibrium effect. Sbank represents sediment input from bank erosion. Bed change is computed with the combined effect of suspended and bed-load transport ([42]):

$$(\mathbf{1} - p')\frac{\partial \mathbf{z}\_b}{\partial t} = a a \rho\_\* (\mathbf{c} - \mathbf{c}\_\*) + (q\_b - q\_{b\*})/L \tag{9}$$

Eqs. (7)–(9) are used for nonuniform sediment transport. The bed-load capacity is computed with ([41]):

$$\rho\_{bk} = 0.0053 \left[ \left( \frac{n'}{n} \right)^{3/2} \frac{\tau\_b}{\tau\_{ck}} - 1 \right]^{2.2} \tag{10}$$

Modeling River Morphodynamic Process Using a Depth-Averaged Computational Model… DOI: http://dx.doi.org/10.5772/intechopen.86692

where φbk ¼ qb <sup>∗</sup> <sup>k</sup>= pbk ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi <sup>γ</sup>ð Þ <sup>s</sup>=<sup>γ</sup> � <sup>1</sup> gd<sup>3</sup> k � � q is a nondimensional bed-load transport capacity, qb <sup>∗</sup> <sup>k</sup> is the equilibrium transport rate of the kth size class of bed load per unit width (kg/m/s), pbk is the bed material gradation, n is the Manning's roughness coefficient for channel bed, <sup>n</sup><sup>0</sup> <sup>¼</sup> <sup>d</sup>1=<sup>6</sup> <sup>50</sup> =20 is the Manning's coefficient corresponding to the grain roughness, τ<sup>b</sup> is the bed shear stress, τck is the critical shear stress determined by <sup>τ</sup>ck <sup>¼</sup> <sup>0</sup>:<sup>03</sup> <sup>γ</sup>ð Þ <sup>s</sup> � <sup>γ</sup> dk phk=pek � �0:<sup>6</sup> , and phk and pek are the hiding and exposure probabilities for the kth size class of bed material, defined as phk <sup>¼</sup> <sup>∑</sup><sup>N</sup> <sup>j</sup>¼<sup>1</sup>pbjdj<sup>=</sup> dk <sup>þ</sup> dj � � and pek <sup>¼</sup> <sup>∑</sup><sup>N</sup> <sup>j</sup>¼<sup>1</sup> pbjdk<sup>=</sup> dk <sup>þ</sup> dj � �.

#### 2.2. Secondary current effect

In curved open channels, the flow is forced to follow a curved path with a variable radius of curvature (Figure 1a). On a bed with a transversal slope (Figure 1b), the bed-load motion is different from that with a stream-wise slope only. The path of a near-bed sediment particle is affected by main flow shear, stream-wise slope, as well as by the gravity component on the transversal direction. Van Bendegom's formula ([37]) was applied to calculate the moving angle of the sediment particle due to the bed slope:

$$\tan \phi = \frac{\sin a - \frac{1}{G} \frac{\partial \zeta}{\partial \dot{\eta}}}{\cos a - \frac{1}{G} \frac{\partial \zeta}{\partial \dot{\eta}}} \tag{11}$$

where

elevation; ρ is the density of water; h is the local water depth; f Cor is the Coriolis parameter; τxx, τxy, τyx, and τyyare depth-integrated Reynolds stresses; and τηx, τηy, τbx, and τby are shear stresses on the water surface and the bed. Free-surface elevation of the flow is calculated by the depth-integrated continuity equation:

Turbulence eddy viscosity is computed with the depth-integrated mixing length

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

dz ¼ Cm

where u<sup>∗</sup> is the shear velocity, κ = 0.41 is the Karman constant, and Cm ≈ 2.34375

Nonuniform suspended and bed-load sediment transport can be simulated. The depth-integrated convection-diffusion equation is solved for the suspended

� ∂

where c is the depth-integrated sediment concentration. The diffusivity coefficient for suspended sediment ε<sup>s</sup> ¼ νt=σ<sup>c</sup> with the Schmidt number 0:5≤σ<sup>c</sup> ≤1. c\* and ω<sup>s</sup> are the sediment transport capacity and settling velocity, and α is a coefficient. The source term Sr represents the dispersion due to the vertical distribution of flow velocity and suspended sediment concentration. Bed load is computed with the

<sup>∂</sup><sup>y</sup> <sup>ε</sup>sh <sup>∂</sup><sup>c</sup> ∂y � �

<sup>L</sup> qb � <sup>q</sup> <sup>∗</sup> <sup>b</sup>

<sup>∂</sup><sup>t</sup> <sup>¼</sup> αωsð Þþ <sup>c</sup> � <sup>c</sup> <sup>∗</sup> qb � qb <sup>∗</sup>

Eqs. (7)–(9) are used for nonuniform sediment transport. The bed-load capacity

n � �<sup>3</sup>=<sup>2</sup> τ<sup>b</sup>

τck � 1

" #<sup>2</sup>:<sup>2</sup>

<sup>φ</sup>bk <sup>¼</sup> <sup>0</sup>:<sup>0053</sup> <sup>n</sup><sup>0</sup>

where cb and qb denote bed-load concentration and transport rate and "\*" denotes capacity. δ is the bed-load layer thickness. The subscripts "bx" and "by" indicate the component of bed load in x and y directions. L is the adaptation length of the bed load representing the non-equilibrium effect. Sbank represents sediment input from bank erosion. Bed change is computed with the combined effect of

u∗ hκ

� �<sup>2</sup> <sup>s</sup>

∂u ∂y þ ∂v ∂x � �<sup>2</sup>

þ

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi <sup>1</sup> � <sup>z</sup> h r� �

> ð ∂U ∂z

<sup>∂</sup><sup>y</sup> <sup>¼</sup> <sup>0</sup> (3)

∂U ∂z

dz≈ 0:267κh (5)

¼ αωsð Þ� c <sup>∗</sup> � c Sr (7)

� � <sup>þ</sup> Sbank <sup>¼</sup> <sup>0</sup> (8)

� �=L (9)

(10)

(4)

(6)

þ

∂h ∂t þ ∂uh ∂x þ ∂vh

Current Practice in Fluvial Geomorphology - Dynamics and Diversity

<sup>þ</sup> <sup>2</sup> <sup>∂</sup><sup>v</sup> ∂y � �<sup>2</sup>

eddy viscosity model:

sediment transport:

∂hc ∂t þ ∂uhc ∂x þ ∂vhc <sup>∂</sup><sup>y</sup> � <sup>∂</sup>

mass conservation equation:

is computed with ([41]):

34

<sup>∂</sup>ð Þ <sup>δ</sup>cb ∂t þ <sup>∂</sup>qbx ∂x þ <sup>∂</sup>qby ∂y þ 1

suspended and bed-load transport ([42]):

<sup>1</sup> � <sup>p</sup><sup>0</sup> ð Þ <sup>∂</sup>zb

vt ¼ l 2 <sup>2</sup> <sup>∂</sup><sup>u</sup> ∂x � �<sup>2</sup>

> <sup>l</sup> <sup>¼</sup> <sup>1</sup> h ð κz

> > ∂U <sup>∂</sup><sup>z</sup> <sup>¼</sup> <sup>1</sup> h

is based on the vertical log distribution of flow velocity ([19]).

<sup>∂</sup><sup>x</sup> <sup>ε</sup>sh <sup>∂</sup><sup>c</sup> ∂x � �

#### Figure 1.

Suspended load and bed-load motion affected by the secondary flow and the gravity. (a) Definition of longitudinal and secondary current velocities. (b) Effect of transverse bed slope and secondary flow. (c) Effect of secondary current on suspended sediment.

Current Practice in Fluvial Geomorphology - Dynamics and Diversity

$$G = f(\theta) = \mathbf{1}. \mathbf{7}\sqrt{\theta} \tag{12}$$

u~l ul

u~<sup>l</sup> � ul ¼ ul

averaged values are

value has been given by ([16])

<sup>~</sup><sup>c</sup> � <sup>c</sup> <sup>¼</sup> <sup>4</sup>:<sup>77</sup> ηδ

DOI: http://dx.doi.org/10.5772/intechopen.86692

1 � ηδ � � <sup>ω</sup>

Zbþδ

Zbþδ

Il ¼ ðη

Ir ¼ ðη

Zbþδ

Zbþδ

computing the dispersion terms; Sr ¼ Dxx þ Dyy,

Dxx <sup>¼</sup> <sup>∂</sup> ∂x ðη

Dyy <sup>¼</sup> <sup>∂</sup> ∂y ðη

source term.

37

2.3. Bank erosion model

u~<sup>r</sup> ¼ 6ul

u~<sup>r</sup> � ur ¼ 6ul

<sup>κ</sup><sup>u</sup> <sup>∗</sup> ω κu<sup>∗</sup>

<sup>¼</sup> <sup>1</sup> <sup>þ</sup> <sup>m</sup> m

Modeling River Morphodynamic Process Using a Depth-Averaged Computational Model…

The difference of the corresponding velocity distribution and the depth-

h r 2 z <sup>h</sup> � <sup>1</sup>

1 þ m m

> h r 2 z <sup>h</sup> � <sup>1</sup>

The sediment concentration distribution is assumed to be the Rouse profile [29]; a simplified model for the difference of the concentration profile and an average

> þ 0:4 � �<sup>1</sup>:<sup>77</sup> <sup>1</sup> � <sup>η</sup>

" #

where ηδ = 0.05 is the relative depth of δ. Eqs. (17)–(19) are employed for

ð Þ <sup>u</sup><sup>~</sup> � <sup>u</sup> ð Þ <sup>~</sup><sup>c</sup> � <sup>c</sup> dz <sup>¼</sup> <sup>∂</sup>

ð Þ <sup>~</sup><sup>v</sup> � <sup>v</sup> ð Þ <sup>~</sup><sup>c</sup> � <sup>c</sup> dz <sup>¼</sup> <sup>∂</sup>

numerically for the entire domain and applied in Eq. (7). lx, ly, rx, and ry are direction vectors of Il and Ir, respectively. The integrals Il and Ir are sediment flux in l and r directions; they are transformed in the x and y direction for computing the

A mass failure would likely occur if a stream bank is high and steep. The failed bank material deposits near the bank toe and then is eroded away by the flow. Depending on geometries and soil properties, river bank failure may have several types: planar, rotational, cantilever, piping-type, and sapping-type ([5]). Planar and rotational failures usually occur to homogeneous, non-layered banks; cantilever failures likely happen to banks with a cohesive top layer and sand and gravel lower layers, while piping- and sapping-type failures most likely occur to the heterogeneous banks where seepage is observed. Osman and Thorne ([26, 27]) analyzed the planar and rotational failures and developed an analytical bank failure model (Figure 2). The bank stability is determined by a factor of safety, defined as

z h � �1=<sup>m</sup>

� 1

η � � <sup>ω</sup> κu ∗ � 1

∂x

∂y

Il � ly þ Ir � ry

<sup>u</sup>~<sup>l</sup> � ulÞð~<sup>c</sup> � <sup>c</sup>Þdz � (22)

<sup>~</sup>vl � vlÞð~<sup>c</sup> � <sup>c</sup>Þdz � (23)

z h � �1=<sup>m</sup>

� � (16)

� � (17)

� � (18)

c (19)

Il � lx þ Ir ð Þ � rx (20)

� � (21)

(15)

α is the angle between the flow direction and the x-axis of the Cartesian coordinate system and θ is the Shields parameter:

$$\theta = \frac{\mathfrak{u}\_\*^2}{\mathfrak{g}\left(\frac{\rho\_\*}{\rho} - 1\right) d\_{50}} \tag{13}$$

The expression of the G function and the coefficient was determined using laboratory experimental data ([37]). When water flows along a curved channel with varying curvatures, the secondary current occurs due to the centrifugal force (Figure 1c). The secondary flow is toward the outer bank of a meander bend in the upper portion of the flow depth and toward the inner bank in the lower portion of the flow. It therefore contributes to moving the net sediment flux in the transversal direction from the outer bank toward the inner bank of the channel systematically. This action erodes the bed near outer bank and deposits on the bed near the inner bank. The main flow is in turn affected by the updated bed topography and the channel pattern. It is not possible to simulate the bed load and bed change in curved channels without considering this process. However, because the depth-integrated model has no direct information about the secondary current, empirical or semianalytical estimation of the secondary flow is used in order to better predict the bed-load motion. The most significant parameter of this problem is the angle between main flow and the near-bed shear stress direction. In the current model, this angle is approximated by ([7])

$$
\tan \delta = 7 \frac{h}{r} \tag{14}
$$

where r is the radius of curvature of the main flow. The error of this formula is about 3% according to [7].

In natural rivers, r is not a given value because it may change with the local flow conditions. In this study, r is computed using the local flow vector directions, the nodal distance, and the mathematical definition: r ¼ ds=dθ. Figure 1b shows the motion of a sediment particle on the bed with a side slope. The gravity pushes the moving particle to move down the transversal slope β with an angle ϕ as estimated by Eq. (11). In the curved channel, the secondary flow pushes the particle moving against the transversal slope by an angle δ [Eq. (14)]. The sediment movement direction computed under the flow and secondary current conditions will be used to determine the bed-load direction in Eq. (8). Equilibrium shall be reached when these two effects cancel each other, and the sediment particles move along the main flow (longitudinal) direction (Figure 1a, b).

Similar to the bed-load sediment, the secondary flow effect for the suspended sediment was also modeled by adding a source term taking into account the net lateral motion of the suspended sediment (Figure 1c). Eq. (7) is a depth-integrated model. In the processes of vertical integration, one has to either assume the vertical variation of the variables is negligible or model the dispersion term to preserve the effect of velocity and sediment profiles on sediment transport. In the second case, the source (dispersion) term in this equation should be non-zero. Computing the dispersion term is, however, complicated, requiring the knowledge of the vertical velocity and suspended sediment profiles. In this study, the vertical variation of the main flow and secondary current is approximated with the power law and linear distribution in the longitudinal (ℓ) and transverse (r) directions ([29]), respectively:

Modeling River Morphodynamic Process Using a Depth-Averaged Computational Model… DOI: http://dx.doi.org/10.5772/intechopen.86692

$$\frac{\tilde{u}\_l}{u\_l} = \frac{1+m}{m} \left(\frac{z}{h}\right)^{1/m} \tag{15}$$

$$
\tilde{u}\_r = \mathfrak{G}u\_l \frac{h}{r} \left( 2\frac{z}{h} - \mathbf{1} \right) \tag{16}
$$

The difference of the corresponding velocity distribution and the depthaveraged values are

$$
\tilde{u}\_l - u\_l = u\_l \left[ \frac{\mathbf{1} + m}{m} \left( \frac{\mathbf{z}}{h} \right)^{1/m} - \mathbf{1} \right] \tag{17}
$$

$$
\tilde{u}\_r - u\_r = 6u\_l \frac{h}{r} \left( 2\frac{z}{h} - 1 \right) \tag{18}
$$

The sediment concentration distribution is assumed to be the Rouse profile [29]; a simplified model for the difference of the concentration profile and an average value has been given by ([16])

$$\tilde{\boldsymbol{\sigma}} - \boldsymbol{\sigma} = \left[ 4.77 \left( \frac{\eta\_{\delta}}{1 - \eta\_{\delta}} \right)^{\frac{\nu}{m\_\*}} \left( \frac{\alpha}{\kappa u\_\*} + 0.4 \right)^{1.77} \left( \frac{1 - \eta}{\eta} \right)^{\frac{\nu}{m\_\*}} - 1 \right] \boldsymbol{\sigma} \tag{19}$$

where ηδ = 0.05 is the relative depth of δ. Eqs. (17)–(19) are employed for computing the dispersion terms; Sr ¼ Dxx þ Dyy,

$$D\_{\rm xx} = \frac{\partial}{\partial \mathbf{x}} \int\_{Z\_b + \delta}^{\eta} (\tilde{u} - u)(\tilde{c} - c)dz = \frac{\partial}{\partial \mathbf{x}} (I\_l \cdot l\_{\mathbf{x}} + I\_r \cdot r\_{\mathbf{x}}) \tag{20}$$

$$D\_{\mathcal{yy}} = \frac{\partial}{\partial \mathbf{y}} \int\_{Z\_{\tilde{\mathbf{z}} + \tilde{\boldsymbol{\sigma}}}}^{\eta} (\tilde{\boldsymbol{v}} - \boldsymbol{v})(\tilde{\boldsymbol{c}} - \boldsymbol{c})d\mathbf{z} = \frac{\partial}{\partial \mathbf{y}} \left( I\_{\mathcal{l}} \cdot l\_{\mathcal{\boldsymbol{y}}} + I\_{\mathcal{r}} \cdot r\_{\mathcal{\boldsymbol{y}}} \right) \tag{21}$$

$$I\_l = \int\_{Z\_b + \delta}^{\eta} (\tilde{u}\_l - u\_l)(\tilde{c} - c)dz\tag{22}$$

$$I\_r = \int\_{Z\_{b^\*} + \delta}^{\eta} (\tilde{\nu}\_l - \nu\_l)(\tilde{c} - c)dz\tag{23}$$

numerically for the entire domain and applied in Eq. (7). lx, ly, rx, and ry are direction vectors of Il and Ir, respectively. The integrals Il and Ir are sediment flux in l and r directions; they are transformed in the x and y direction for computing the source term.

#### 2.3. Bank erosion model

A mass failure would likely occur if a stream bank is high and steep. The failed bank material deposits near the bank toe and then is eroded away by the flow. Depending on geometries and soil properties, river bank failure may have several types: planar, rotational, cantilever, piping-type, and sapping-type ([5]). Planar and rotational failures usually occur to homogeneous, non-layered banks; cantilever failures likely happen to banks with a cohesive top layer and sand and gravel lower layers, while piping- and sapping-type failures most likely occur to the heterogeneous banks where seepage is observed. Osman and Thorne ([26, 27]) analyzed the planar and rotational failures and developed an analytical bank failure model (Figure 2). The bank stability is determined by a factor of safety, defined as

<sup>G</sup> <sup>¼</sup> <sup>f</sup>ð Þ¼ <sup>θ</sup> <sup>1</sup>:<sup>7</sup> ffiffiffi

<sup>θ</sup> <sup>¼</sup> <sup>u</sup><sup>2</sup>

g <sup>ρ</sup><sup>s</sup> <sup>ρ</sup> � 1 � �

The expression of the G function and the coefficient was determined using laboratory experimental data ([37]). When water flows along a curved channel with varying curvatures, the secondary current occurs due to the centrifugal force (Figure 1c). The secondary flow is toward the outer bank of a meander bend in the upper portion of the flow depth and toward the inner bank in the lower portion of the flow. It therefore contributes to moving the net sediment flux in the transversal direction from the outer bank toward the inner bank of the channel systematically. This action erodes the bed near outer bank and deposits on the bed near the inner bank. The main flow is in turn affected by the updated bed topography and the channel pattern. It is not possible to simulate the bed load and bed change in curved channels without considering this process. However, because the depth-integrated model has no direct information about the secondary current, empirical or semianalytical estimation of the secondary flow is used in order to better predict the bed-load motion. The most significant parameter of this problem is the angle between main flow and the near-bed shear stress direction. In the current model,

tan δ ¼ 7

where r is the radius of curvature of the main flow. The error of this formula is

In natural rivers, r is not a given value because it may change with the local flow conditions. In this study, r is computed using the local flow vector directions, the nodal distance, and the mathematical definition: r ¼ ds=dθ. Figure 1b shows the motion of a sediment particle on the bed with a side slope. The gravity pushes the moving particle to move down the transversal slope β with an angle ϕ as estimated by Eq. (11). In the curved channel, the secondary flow pushes the particle moving against the transversal slope by an angle δ [Eq. (14)]. The sediment movement direction computed under the flow and secondary current conditions will be used to determine the bed-load direction in Eq. (8). Equilibrium shall be reached when these two effects cancel each other, and the sediment particles move along the main

Similar to the bed-load sediment, the secondary flow effect for the suspended sediment was also modeled by adding a source term taking into account the net lateral motion of the suspended sediment (Figure 1c). Eq. (7) is a depth-integrated model. In the processes of vertical integration, one has to either assume the vertical variation of the variables is negligible or model the dispersion term to preserve the effect of velocity and sediment profiles on sediment transport. In the second case, the source (dispersion) term in this equation should be non-zero. Computing the dispersion term is, however, complicated, requiring the knowledge of the vertical velocity and suspended sediment profiles. In this study, the vertical variation of the main flow and secondary current is approximated with the power law and linear distribution in

the longitudinal (ℓ) and transverse (r) directions ([29]), respectively:

h r

nate system and θ is the Shields parameter:

Current Practice in Fluvial Geomorphology - Dynamics and Diversity

this angle is approximated by ([7])

flow (longitudinal) direction (Figure 1a, b).

about 3% according to [7].

36

α is the angle between the flow direction and the x-axis of the Cartesian coordi-

∗

d<sup>50</sup>

θ

<sup>p</sup> (12)

(13)

(14)

Current Practice in Fluvial Geomorphology - Dynamics and Diversity

$$f\_s = \frac{F\_r}{F\_d} \tag{24}$$

stress and a coefficient which is also related to the critical shear stress. In Osman and Thorne's model [26, 27], the bank surface erosion rate, ε, was proportional to the difference of the bed shear stress, τ, and the bank critical stress, τc, normalized

Modeling River Morphodynamic Process Using a Depth-Averaged Computational Model…

<sup>ε</sup> <sup>¼</sup> <sup>k</sup> <sup>τ</sup> � <sup>τ</sup><sup>c</sup> τc

Field data of almost 200 sites ([15]) indicated that k can be expressed by another

Following a bank failure and retreat event, the mesh lines representing the bank boundaries need to be moved to an updated bank location resulting in a moving boundary problem. Computational mesh should be stretched to widen the computing domain for the widened river. After a bank mesh line is moved, internal mesh line adjustment will be necessary to redistribute the internal nodes in the updated

where k is the bank erosion rate which is a function of the critical stress:

<sup>k</sup> <sup>¼</sup> <sup>223</sup> � <sup>10</sup>�4τce

computational domain (widened channel). Once a mesh is stretched, the

General model execution procedure. Bank erosion loop is computed less frequent.

discretization of the computational domain should be updated. This procedure is called dynamic meshing. One has to recompute all the numerical parameters and differential operators again every time a mesh stretching is performed. Interpolation of the computational results from the previous mesh to the stretched new one is required before recomputing the flow. Because bank erosion process is much slower than the flow, sediment transport, and bed change, it can be computed with a much

(29)

�0:13τ<sup>c</sup> (30)

<sup>k</sup> <sup>¼</sup> <sup>0</sup>:1τ�0:<sup>5</sup> <sup>c</sup> (31)

by the critical stress:

Figure 3.

39

function of critical shear stress:

DOI: http://dx.doi.org/10.5772/intechopen.86692

where Fr and Fd are the resisting and driving forces, respectively. When fs < 1, a bank mass failure is expected to occur.

In Osman and Thorne's model, a bank has an initial slope; after the first collapse occurs, a new slope will be established. The bank will then keep this slope, and the subsequent mass failures will not change the slope (parallel retreat). Considering that river banks for any study have been experiencing bank failures for a long time, the bank slope observed in the field is likely the bank mass failure slope. It is therefore assumed that the bank failure slope, β, is a known value and only the parallel retreat processes need to be simulated. Under this condition, the lateral bank retreat distance with a constant bank slope is calculated by

$$BW = \frac{H - H'}{\tan \beta} \tag{25}$$

The critical ratio of the new and old bank height determined by

$$\frac{H}{H'} = \frac{1}{2} \left[ \frac{\alpha\_2}{\alpha\_1} + \sqrt{\left(\frac{\alpha\_2}{\alpha\_1}\right)^2 + 4} \right] \tag{26}$$

$$
\rho\_1 = \cos \beta \sin \beta - \cos^2 \beta \tan \varphi \tag{27}
$$

$$
\alpha\_2 = \mathbf{2}(\mathbf{1} - K) \frac{c}{\chi\_s H'} \tag{28}
$$

will be used to test if a mass failure occurs: if the ratio of computed H and H<sup>0</sup> is higher than that from Eq. (26), a bank failure is computed. K is the tension crack index: the ratio of observed tension crack depth to bank height. Usually, the failed material deposits first near the bank toe and then is disaggregated and eroded away by the flow. In the current approach, the failed bank material is considered as a supply source to the bed load. Since the time step for the bank erosion is much larger than that of the sediment transport, the source term representing this sediment supply from bank erosion is set uniform through the next bank erosion time step. This supply will result in higher near-bank sediment concentration or bed load. If the bank erosion is too fast, near-bank bed elevation would increase to slow down the bank erosion. Cohesive material erosion is proportional to excessive shear

Figure 2. Mode of bank mass failure (after Osman and Thorne, [26, 27]).

Modeling River Morphodynamic Process Using a Depth-Averaged Computational Model… DOI: http://dx.doi.org/10.5772/intechopen.86692

stress and a coefficient which is also related to the critical shear stress. In Osman and Thorne's model [26, 27], the bank surface erosion rate, ε, was proportional to the difference of the bed shear stress, τ, and the bank critical stress, τc, normalized by the critical stress:

$$
\varepsilon = k \frac{\tau - \tau\_c}{\tau\_c} \tag{29}
$$

where k is the bank erosion rate which is a function of the critical stress:

$$k = 223 \times 10^{-4} \tau\_c e^{-0.13 \tau\_c} \tag{30}$$

Field data of almost 200 sites ([15]) indicated that k can be expressed by another function of critical shear stress:

$$k = 0.1\tau\_c^{-0.5} \tag{31}$$

Following a bank failure and retreat event, the mesh lines representing the bank boundaries need to be moved to an updated bank location resulting in a moving boundary problem. Computational mesh should be stretched to widen the computing domain for the widened river. After a bank mesh line is moved, internal mesh line adjustment will be necessary to redistribute the internal nodes in the updated computational domain (widened channel). Once a mesh is stretched, the discretization of the computational domain should be updated. This procedure is called dynamic meshing. One has to recompute all the numerical parameters and differential operators again every time a mesh stretching is performed. Interpolation of the computational results from the previous mesh to the stretched new one is required before recomputing the flow. Because bank erosion process is much slower than the flow, sediment transport, and bed change, it can be computed with a much

Figure 3. General model execution procedure. Bank erosion loop is computed less frequent.

fs <sup>¼</sup> Fr Fd

bank mass failure is expected to occur.

Figure 2.

38

Mode of bank mass failure (after Osman and Thorne, [26, 27]).

where Fr and Fd are the resisting and driving forces, respectively. When fs < 1, a

In Osman and Thorne's model, a bank has an initial slope; after the first collapse occurs, a new slope will be established. The bank will then keep this slope, and the subsequent mass failures will not change the slope (parallel retreat). Considering that river banks for any study have been experiencing bank failures for a long time, the bank slope observed in the field is likely the bank mass failure slope. It is therefore assumed that the bank failure slope, β, is a known value and only the parallel retreat processes need to be simulated. Under this condition, the lateral

BW <sup>¼</sup> <sup>H</sup> � <sup>H</sup><sup>0</sup>

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ω2 ω1 � �<sup>2</sup>

þ 4

3

tan <sup>β</sup> (25)

5 (26)

β tan φ (27)

<sup>γ</sup>sH<sup>0</sup> (28)

bank retreat distance with a constant bank slope is calculated by

Current Practice in Fluvial Geomorphology - Dynamics and Diversity

H <sup>H</sup><sup>0</sup> <sup>¼</sup> <sup>1</sup> 2

The critical ratio of the new and old bank height determined by

4

ω2 ω1 þ

2 s

<sup>ω</sup><sup>1</sup> <sup>¼</sup> cos <sup>β</sup> sin <sup>β</sup> � cos <sup>2</sup>

<sup>ω</sup><sup>2</sup> <sup>¼</sup> 2 1ð Þ � <sup>K</sup> <sup>c</sup>

will be used to test if a mass failure occurs: if the ratio of computed H and H<sup>0</sup> is higher than that from Eq. (26), a bank failure is computed. K is the tension crack index: the ratio of observed tension crack depth to bank height. Usually, the failed material deposits first near the bank toe and then is disaggregated and eroded away by the flow. In the current approach, the failed bank material is considered as a supply source to the bed load. Since the time step for the bank erosion is much larger than that of the sediment transport, the source term representing this sediment supply from bank erosion is set uniform through the next bank erosion time step. This supply will result in higher near-bank sediment concentration or bed load. If the bank erosion is too fast, near-bank bed elevation would increase to slow down the bank erosion. Cohesive material erosion is proportional to excessive shear

(24)

larger time than that for the flow and sediment. This strategy can save a lot of computing time. Figure 3 briefly illustrates the corresponding computation procedure.
