Modeling River Morphodynamic Process Using a Depth-Averaged Computational Model and an Application to a Mountain River

Abstract

Bank erosion is a dominant river morphodynamic process resulting in encroaching valuable farming land and channel migration. Prediction of bank erosion and channel migration requires understanding of the morphodynamics of the entire river system. Numerical modeling is an ideal method for this task. However, models with full capabilities and applications on complex real-world problems are rare. In this study the finite element-based computational model, CCHE2D, and its flow, sediment transport, and bank erosion modules are introduced. The model is capable of simulating unsteady flows with nonuniform sediment transport and cohesive/non-cohesive material bank erosion. The effects of helical secondary current on sediment transport induced by flow curvatures are reflected in both bed load and suspended sediment formulations. This model is validated using multiple sets of experimental data and applied to bank erosion problems of the Chuoshui River, a real-world mountain river in Taiwan. Characterized by typhoon floods, steep channel slopes, and high sediment load and mobility, this river often exhibits a braided pattern consisting of multiple curved channels. Channel bed change and bank erosion caused by 10 years of typhoon floods in a selected reach have been simulated, and the computed bank erosion results agreed with the field observation.

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.

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 mechanic-based and multilayered bank stability model of Simon et al. [33, 35]. The coupled 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 x}+v\frac{\partial u}{\partial y}=-g\frac{\partial \eta }{\partial x}+\frac{1}{h}\left(\frac{\partial h{\tau }_{\mathit{xx}}}{\partial x}+\frac{\partial h{\tau }_{\mathit{xy}}}{\partial y}\right)+\frac{{\tau }_{\mathit{\eta x}}-{\tau }_{\mathit{bx}}}{\mathit{\rho h}}+f_{\mathit{Cor}}v$$E1

$$\frac{\partial v}{\partial t}+u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y}=-g\frac{\partial \eta }{\partial y}+\frac{1}{h}\left(\frac{\partial h{\tau }_{\mathit{yx}}}{\partial x}+\frac{\partial h{\tau }_{\mathit{yy}}}{\partial y}\right)+\frac{{\tau }_{\mathit{\eta y}}-{\tau }_{\mathit{by}}}{\mathit{\rho h}}-f_{\mathit{Cor}}u$$E2

where u and v are depth-integrated velocity components in x and y directions, respectively; t is the time; g is the gravitational acceleration; $\eta$ is the water surface elevation; $\rho$ is the density of water; h is the local water depth; $f_{\mathit{Cor}}$ is the Coriolis parameter; ${\tau }_{\mathit{xx}}$, ${\tau }_{\mathit{xy}}$, ${\tau }_{\mathit{yx}}$, and ${\tau }_{\mathit{yy}}$ are depth-integrated Reynolds stresses; and ${\tau }_{\mathit{\eta x}}$, ${\tau }_{\mathit{\eta y}}$, ${\tau }_{\mathit{bx}}$, and ${\tau }_{\mathit{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 \mathit{uh}}{\partial x}+\frac{\partial \mathit{vh}}{\partial y}=0$$E3

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 x}\right)}^2+2{\left(\frac{\partial v}{\partial y}\right)}^2+{\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)}^2+{\left(\frac{\overline{\partial U}}{\partial z}\right)}^2}$$E4

$$\overline{l}=\frac{1}{h}\int \mathit{\kappa z}\sqrt{\left(1-\frac{z}{h}\right)}\mathit{dz}\approx 0.267\mathit{\kappa h}$$E5

$$\frac{\overline{\partial U}}{\partial z}=\frac{1}{h}\int \frac{\partial U}{\partial z}\mathit{dz}=C_m\frac{u_{\ast }}{\mathit{h\kappa }}$$E6

where u_∗ is the shear velocity, κ = 0.41 is the Karman constant, and $C_m\approx$ 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 \mathit{hc}}{\partial t}+\frac{\partial \mathit{uhc}}{\partial x}+\frac{\partial \mathit{vhc}}{\partial y}-\frac{\partial }{\partial x}\left[{\epsilon }_{s}h\frac{\partial c}{\partial x}\right]-\frac{\partial }{\partial y}\left[{\epsilon }_{s}h\frac{\partial c}{\partial y}\right]={\mathit{\alpha \omega }}_s\left(c_{\ast }-c\right)-S_r$$E7

where c is the depth-integrated sediment concentration. The diffusivity coefficient for suspended sediment ${\epsilon }_s={\nu }_t/{\sigma }_c$ with the Schmidt number $0.5\le {\sigma }_c\le 1$. c_* and ω_s are the sediment transport capacity and settling velocity, and α is a coefficient. The source term S_r 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 \left(\delta c_b\right)}{\partial t}+\frac{\partial q_{\mathit{bx}}}{\partial x}+\frac{\partial q_{\mathit{by}}}{\partial y}+\frac{1}{L}\left(q_b-q_{\ast b}\right)+S_{\mathit{\text{bank}}}=0$$E8

where c_b and q_b 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. S_bank represents sediment input from bank erosion. Bed change is computed with the combined effect of suspended and bed-load transport ([42]):

$$\left(1-p^{\prime }\right)\frac{\partial z_b}{\partial t}={\mathit{\alpha \omega }}_s\left(c-c_{\ast }\right)+\left(q_b-q_{b\ast }\right)/L$$E9

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

$${\phi }_{\mathit{bk}}=0.0053{\left[{\left(\frac{n^{\prime }}{n}\right)}^{3/2}\frac{{\tau }_b}{{\tau }_{\mathit{ck}}}-1\right]}^{2.2}$$E10

where ${\phi }_{\mathit{bk}}=q_{b\ast k}/\left[p_{\mathit{bk}}\sqrt{\left({\gamma }_s/\gamma -1\right){\mathit{gd}}_k^3}\right]$ is a nondimensional bed-load transport capacity, $q_{b\ast k}$ is the equilibrium transport rate of the kth size class of bed load per unit width (kg/m/s), $p_{\mathit{bk}}$ is the bed material gradation, n is the Manning’s roughness coefficient for channel bed, $n^{\prime }=d_{50}^{1/6}/20$ is the Manning’s coefficient corresponding to the grain roughness, ${\tau }_b$ is the bed shear stress, ${\tau }_{\mathit{ck}}$ is the critical shear stress determined by ${\tau }_{\mathit{ck}}=0.03\left({\gamma }_s-\gamma \right)d_k{\left(p_{\mathit{hk}}/p_{\mathit{ek}}\right)}^{0.6}$, and $p_{\mathit{hk}}$ and $p_{\mathit{ek}}$ are the hiding and exposure probabilities for the kth size class of bed material, defined as $p_{\mathit{hk}}={\sum }_{j=1}^N{p}_{\mathit{bj}}{d}_j/\left(d_k+d_j\right)$ and $p_{\mathit{ek}}={\sum }_{j=1}^N{p}_{\mathit{bj}}{d}_k/\left(d_k+d_j\right)$.

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:

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.

$$tan\varphi =\frac{sin\alpha -\frac{1}{G}\frac{\partial \zeta }{\partial y}}{cos\alpha -\frac{1}{G}\frac{\partial \zeta }{\partial x}}$$E11

where

$$G=f\left(\theta \right)=1.7\sqrt{\theta }$$E12

$\alpha$ is the angle between the flow direction and the x-axis of the Cartesian coordinate system and $\theta$ is the Shields parameter:

$$\theta =\frac{u_{\ast }^2}{g\left(\frac{{\rho }_s}{\rho }-1\right)d_{50}}$$E13

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 semi-analytical 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}$$E14

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=\mathit{ds}/\mathit{d\theta }$. 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 $\beta$ with an angle $\varphi$ as estimated by Eq. (11). In the curved channel, the secondary flow pushes the particle moving against the transversal slope by an angle $\delta$ [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 ($\ell$) and transverse (r) directions ([29]), respectively:

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

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

The difference of the corresponding velocity distribution and the depth-averaged values are

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

$${\tilde{u}}_r-u_r=6u_l\frac{h}{r}\left(2\frac{z}{h}-1\right)$$E18

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{c}-c=\left[4.77{\left(\frac{{\eta }_{\delta }}{1-{\eta }_{\delta }}\right)}^{\frac{\omega }{\kappa u_{\ast }}}{\left(\frac{\omega }{\kappa u_{\ast }}+0.4\right)}^{1.77}{\left(\frac{1-\eta }{\eta }\right)}^{\frac{\omega }{\kappa u_{\ast }}}-1\right]c$$E19

where ${\eta }_{\delta }$ = 0.05 is the relative depth of $\delta$. Eqs. (17)–(19) are employed for computing the dispersion terms; $S_r=D_{\mathit{xx}}+D_{\mathit{yy}}$,

$$D_{\mathit{xx}}=\frac{\partial }{\partial x}{\int }_{Z_b+\delta }^{\eta }\left(\tilde{u}-u\right)\left(\tilde{c}-c\right)\mathit{dz}=\frac{\partial }{\partial x}\left(I_l\cdot l_x+I_r\cdot r_x\right)$$E20

$$D_{\mathit{yy}}=\frac{\partial }{\partial y}{\int }_{Z_b+\delta }^{\eta }\left(\tilde{v}-v\right)\left(\tilde{c}-c\right)\mathit{dz}=\frac{\partial }{\partial y}\left(I_l\cdot l_y+I_r\cdot r_y\right)$$E21

$$I_l={\int }_{Z_b+\delta }^{\eta }\left({\tilde{u}}_l-u_l\right)\left(\tilde{c}-c\right)\mathit{dz}$$E22

$$I_r={\int }_{Z_b+\delta }^{\eta }\left({\tilde{v}}_l-v_l\right)\left(\tilde{c}-c\right)\mathit{dz}$$E23

numerically for the entire domain and applied in Eq. (7). l_x, l_y, r_x, and r_y are direction vectors of I_l and I_r, respectively. The integrals I_l and I_r 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

Figure 2.

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

$$f_s=\frac{F_r}{F_d}$$E24

where $F_r$ and $F_d$ are the resisting and driving forces, respectively. When $f_s<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

$$\mathit{BW}=\frac{H-H^{\prime }}{tan\beta }$$E25

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

$$\frac{H}{H^{\prime }}=\frac{1}{2}\left[\frac{{\omega }_2}{{\omega }_1}+\sqrt{{\left(\frac{{\omega }_2}{{\omega }_1}\right)}^2+4}\right]$$E26

$${\omega }_1=cos\beta sin\beta -{cos}^2\beta tan\phi$$E27

$${\omega }_2=2\left(1-K\right)\frac{c}{{\gamma }_s{H}^{\prime }}$$E28

will be used to test if a mass failure occurs: if the ratio of computed H and H′ 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 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, $\tau$, and the bank critical stress, ${\tau }_c$, normalized by the critical stress:

$$\epsilon =k\frac{\tau -{\tau }_c}{{\tau }_c}$$E29

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}$$E30

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}$$E31

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

Figure 3.

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

3. Validation and application results of the model

3.1. Validation of bed morphological change and bank erosion models

The sediment transport and bed morphological change simulation model was tested using physical model data. Four experiment test cases with a different channel geometry, curvature, flow condition, and sediment size distribution ([36]) were simulated. The sediments of these three cases are uniform. Figure 4 shows the computed final water depth of the case with a 180° bend and comparison of the simulation to the measurement. Red and blue color indicates deep and shallow water depth, respectively. Table 1 shows the parameters of the experimental flume and flow. Channel plane geometry, slope, sediment size (in bed and from inlet) distribution, and Manning’s n are needed to run the simulations; the flow rate and water depth are used for the upstream and downstream boundary conditions. Because the initial bed is horizontal and flat and the water depth almost constant, the resulting water depth distribution indicated more erosion along the outer bank and deposition along the inner bank. Although differing, the computed bed elevations along the channel agree reasonably well with the observation. The magnitude of the predicted erosion and deposition in the channels agreed very well with the measurement; the second water depth peak of the bed variation along the outer bank has some difference from the observed.

Figure 4.

Computed water depth and comparison of numerical results (curve) and experimental data (Case 4).

Parameters	Case 4 LFM flume
Discharge (m^3/s)	0.17
Flume width (m)	1.7
Water depth (m)	0.2
Flow velocity (m/s)	0.5
Water surface slope (^0/00)	1.8
Chezy coefficient (m^1/2/s)	26.4
Manning coefficient (n)	0.0280
D50 (mm)	0.78
Sediment transport rate (m^2/s)	13 x 10^−6
Bend radius (m)	4.25
Bend length (m)	13.35

Table 1.

The conditions and parameters of the physical model.

A qualitative study of sediment transport in conjunction with the bank erosion simulation is also presented (Figure 5). A river channel with the sine-generated shape, constant bed roughness, channel width, and longitudinal slope was generated for the simulation. The initial bed erosion simulation was performed with fixed banks. The bank erosion simulation started after the bed erosion has been performed for a while. When the bank erosion simulation is completed for one time step, the mesh bank line and internal points are shifted and the model re-discretized, and so on, as indicated in Figure 3.

Figure 5.

Simulated bank erosion and channel morphologic change using bank-full discharge (the color contour indicates flow velocity magnitude).

The colors in Figure 5 indicate the flow velocity magnitude. One sees the phase difference between the shape of the channel bends and the velocity distribution. The highest flow velocity shifts downstream. In the process of bank erosion, the outer bank line retreats gradually, and the main channel of this bend shifts accordingly; the cross-section form of the channel also changes particularly at the beginning stage, and the water depth near the outer bank becomes larger, while that near the inner bank becomes smaller. This change makes it possible to form a point bar near the inner bank (Figure 5); then the point bar later becomes dry. Although the distance of the two banks increases, the width of the wetted channel remained approximately the same. Another feature of the simulated results is that when the main channel moves toward the outer bank due to bank erosion, a small channel near the inner bank is formed behind the point bar (Figure 5d, e). This probably is because the small channel shortcuts from one bend to the next, the local water surface slope, and sediment transport capacity are relatively large. This phenomenon appears also in some natural rivers ([21]).

3.2. Application of the bank erosion model to a field case in Chuoshui River

Chuoshui River is in the middle of the Taiwan inland which is located in the South China Sea across the Taiwan Strait. Originated from the central mountains, the river forms a large alluvial fan and then empties into the South China Sea. The channel slope in the mountain area is very steep. The valley of the river is wide, and a typical braided river pattern with multiple curved sub-channels can be observed from aerial photos and satellite imagery. The study reach is situated at the connection part of the mountain and the alluvial fan of the river. The channel slope is about 0.0069 for the mountain part, and it reduces suddenly to about 0.0041 over the alluvial fan. The hydrology is dominated by seasonal typhoon events and a large amount of sediments from the mountain watershed. The characteristic of braided river varies downstream somewhat, and the number of sub-channels decreases over the alluvial fan.

One should recognize that the predictability of bank erosion is limited by the facts that (1) not all the processes are understood and formulated accurately, (2) collecting field data necessary for the analysis is extremely difficult and costly, and (3) the accuracy of the flow simulation is affected by the modeling methodologies and computer capacity. When a real-world bank erosion problem is modeled, one focuses on the dominant processes and carries out appropriate calibration and validation. These processes could be related to several parameters: sediment properties of bed materials and bank, including bank slope, height, and bank material erodibility, as well as the conditions of the flow in the river channel (shear stress, water depth, channel curvature, etc.).

The flow discharge increases greatly, particularly during typhoon seasons. The multiple channels become a single one only when the discharge is very large during typhoon seasons. Due to the nature of the channel pattern, the main channel and secondary channels in the study reach change courses randomly and quickly. Sediment transport is dominated by the pattern of the flow discharge. The computational model, CCHE2D, was applied to simulate the bank erosion process in one reach of the river from Mingchu Bridge (CS 106.5) to Zhongsha Bridge (CS 52), a 26 km stretch (Figure 6).

Figure 6.

Study reach of Chuoshui River.

Even in a braided river, each sub-channel is a curved one. Because sediment transport in curved channels is affected by the secondary current, creating a lateral sediment motion and channel change, the computational model should include this mechanism to reflect the realistic transport processes. Figure 6 shows the air photo of this reach in 2007, where the nature of the braided river is clearly seen. The flow discharges for this figure are unknown. It is certain the braided river process is very active. The upstream part of the Chuoshui River between Mingchu Bridge (CS 106.5) and Zhangyun Bridge (CS 86.5) was used as a test site for the bank erosion model.

Field bed material samples taken in July 2004 were used as the initial bed composition for the study reach. The average sediment compositions in three sub-reaches are shown in Figure 7 and Table 2. One notes that the sediment sizes range from 0.283 to 282 mm. The trend of sediment particle size decreases downstream (CS 55–CS 70) is quite significant. Particularly the portion of coarse particles decreases more (Figure 7). To gain computational efficiency, the original measured sediment distribution data was simplified from nine size classes to six.

Figure 7.

Initial bed compositions.

CS range	D10	D20	D30	D40	D50	D60	D70	D80	D90
86.5–106.5	0.396	0.94	1.91	5.029	13.69	29.424	62.532	143.921	281.966
70–86.5	0.404	0.981	1.826	3.452	8.784	22.682	46.998	124.463	217.7
55–70	0.283	0.597	1.096	1.91	3.348	5.775	11.509	27.495	82.315

Table 2.

Specified bed materials in three channel reaches (d/mm).

The flow discharge in the Chuoshui River is highly variable, ranging from almost zero in dry seasons to more than 20,000 cms in some typhoon seasons. Because sediment transport is insignificant to channel change when the flow discharges are low, the bank erosion study was performed only when the flow rate is high (>4000 cms). Considering the bank-full discharge is 5700 cms in this channel, this criterion of simplification has included most significant flows. Figure 8 shows the simplified hydrograph including most of the typhoon events from June 8, 1998, to October 8, 2007. The corresponding downstream boundary condition, water surface elevation at the Zhongsha Bridge (CS 55), is also shown. Rating curves for the sediment transport rate were used for sediment boundary conditions with wash load being removed. The accumulated total time of these high flows is approximately 4 days and 20 hours.

Figure 8.

Discharge and water stage hydrograph at Zhangyun and Zhongsha Bridge (CS 55), respectively.

The sediment discharge hydrograph at Mingchu (CS 106.5) and Zhangyun (CS 86.1) was also filtered accordingly to remove low flow events. Because there is little information about the sediment size composition for the estimated sediment load, it is assumed that (1) 95% of sediment load is suspended and 5% is bed load and that (2) 80% of suspended load coming from upstream is wash load. Considering that the sediment composition would be a function of the flow discharge, large sediment particles can be moved only when the discharge is large, and fine particle can be moved by any flow; the composition of bed sediment can be adjusted by the erosion and deposition process.

The time step for bank erosion was set to be 1.0 hour, while that for the flow and sediment transport was 30 seconds. The critical stress used was consistent with the field data ([34]) for low cohesive bank materials.

The bank erosion was estimated using the difference between 1998 and 2007 DEM data. Figure 9 shows the measured and computed bed morphologic change of the second reach (II) from 1998 to 2007. The measured bed change is presented over the 1998 aerial photo with the initial mesh boundaries indicated by white lines (Figure 9a). The computed bed change is also presented in a similar fashion (Figure 9b), except that the part outside the computational mesh is contour lines rather than color shading. The initial mesh boundaries are indicated with white lines, while the final mesh boundaries due to bank erosion are presented with purple lines. The difference between these two colored lines represents bank erosion. As indicated in the figure, the simulated bank erosions are very close to the observed, particularly for those on the left bank. The lateral movement of the bank line ranges in several 100 meters. The most significant bank erosions occurred at the left bank. The white circles indicate the two significant bank erosion zones. The general shape and area of the simulated bank erosions are similar to the observed. The maximum erosion distance normal to the left bank line is more than 800 meters.

Figure 9.

Comparison of computed bed change and bank erosion with observed.

The simulated bed change has differences from the observed although they are generally consistent. Downstream of CS 70, the computed bed change is consistent with the measured. The location of the deposition and erosion is correctly predicted with the simulated results having a little more deposition. Near the entrance of the reach, the bed change is dominated by degradation. The computed results show two separated channels, one is being eroded and the other has deposition.

The big point bar indicated in Figure 9 was eroded in the numerical simulation, contributed a lot of sediment to its downstream, and affected the simulation results. However, it was found later that the sediment of this 6 meter high point bar could have been taken away by sand miners. To certain extend, this attributes why more deposition is predicted by the numerical simulations.

To illustrate more clearly the bank erosion simulation, the computed bank lines are plotted together with the measured bed change in two cross sections (Figure 10). The green lines represent the bed cross-section profile of 1998; the red lines represent the profile of 2007. The bank heights at the bank erosion zones were more than 6 meters. It is seen that the major bank erosion in the channel was reasonably predicted by the model. Although the simulated location and amount of the bank erosions do not match exactly to the observation, the general trend of the bank erosion simulation is quite satisfactory. The observed bank erosion at CS 70 and CS 78 were about 500 and 800 m, respectively. As discussed earlier, the model-predicted incision in the main channel was less than the observed. The error is mainly attributed to lacking of desirable data of sediment transport. Secondly, the Chuoshui River in the study reach is of braided pattern with several major branches. Even the general trend of channel aggradation/degradation can be simulated; the sedimentation trend in each branch is difficult to control. More research is necessary.

Figure 10.

Bank erosion comparison. The measured and simulated bed elevation change in CS 70 and CS 78 shows that the bank erosion was simulated well. The incision of the channel thalweg in these two sections was not captured.

4. Major outcomes and conclusions

Morphodynamics of fluvial systems is complex involving channel bed change, bank erosion, and channel migration, and it results in soil loss, water quality deterioration, and property damages. Numerical models can be applied to simulate the system behavior by considering involved key physical mechanisms and processes, such as main and secondary flow, sediment transport processes and bank slope mass failure, etc.

The capabilities for simulating the secondary helical flow effects on suspended sediment and bed-load sediment transport have been developed and implemented to the CCHE2D model. The vertical profile for the main velocity and the secondary helical current were assumed to be the power law and linear distribution, respectively. Rouse’s distribution for suspended sediment concentration was adopted. For general applications, the curvature of the flow instead of the channel was used for the helical flow calculation. The bank toe and surficial erosion and mass failure mechanisms have also been developed with the mass wasted bank materials being transported as bed load. The current model was designed for banks with cohesive and homogeneous materials. The mesh stretching technique was developed and used to adjust dynamically the moving boundary, internal mesh nodal position, and associated interpolation. These are important to simulate rivers with significant bank line movement due to erosion.

Several sets of curved channel experimental data with different channel geometries, flow rates, sediment sizes, etc. were utilized to validate the developed sediment transport and morphodynamic simulation capabilities in good agreement. Bank erosion capabilities were tested first using a sine-generated channel and then the field case of Chuoshui River, Taiwan. The developed dynamic meshing method handled the moving boundary problem satisfactorily. The simulated and observed bank retreats in the studied Chuoshui River reach can be 500–800 m, which is agreed reasonably well. Because bank erosion occurred mainly in typhoon seasons, simulations used only flow discharges larger than 4000 cms. The computed bed change and bank erosion in one reach of this highly mobile braided river were compared with reasonable agreements to observations.

Acknowledgments

This work is supported by the Water Planning Institute, Department of Water Resources, Taiwan, the project of USDA Agricultural Research Service under Specific Research Agreement No. 58-6060-8-008, monitored by the National Sedimentation Laboratory, and the University of Mississippi.