*Article* **Numerical Simulation of the Transportation of Cohesive Bank-Collapsed Materials in a Sharply Curved Channel**

**Guosheng Duan 1,2, Haifei Liu 2,3,\*, Dongdong Shao <sup>2</sup> , Wei Yang <sup>2</sup> , Zhiwei Li <sup>4</sup> , Chen Wang <sup>5</sup> , Shuo Chang <sup>5</sup> and Yu Ding 1,2,\***


**Abstract:** This study presents an integrated analysis of cohesive bank-collapsed material transportation in a high-curvature channel with a non-cohesive riverbed. A numerical model was established to simulate the erosion and transportation of collapsed materials in a 180◦ U-bend channel after verification. The novel aspect of this study is that the quantities of the collapsed materials that transformed into suspended and bed loads were comprehensively analyzed. The results show that finer collapsed sediments were only transformed into suspended loads after being eroded, while the coarser particles transformed into both suspended loads and bed loads. When the flow charge was 30 L/s, the quantity of collapsed materials (S1 and S2) that transported downstream was smaller, and coarser materials transformed into suspended loads with a ratio of 88.12–99.86% and bed loads with a ratio of 11.18–0.14%. When the flow charge was 55 L/s, due to the greater shear stress, the quantity of collapsed materials (S1 and S2) that transported downstream was greater, and the ratio ranged from 46.65% to 49.88% and from 50.12% to 53.35%, respectively. This research provides theoretical and practical benefits that reveal the mechanisms of channel bend evolution.

**Keywords:** bank collapse; sediment transportation; numerical simulation; curved channel; cohesive

## **1. Introduction**

Riverbank collapse is a key process in river morphodynamics that can affect channel mobility, floodplain evolution, and pollution transportation. Large amounts of sediments come into alluvial rivers, leading to a series of social and environmental problems, including farmland loss, embankment destruction, river turbidity, and river eutrophication [1,2].

Given the importance of riverbank collapse, it is not surprising that many studies have been carried out on this subject in these past decades. For cohesive riverbanks, one focus of recent work has been the mechanism of riverbank collapse and the relative influence of the factors that control mass failure [3]. In such studies, collapse processes were divided into three steps: (1) bank toe erosion, (2) tension cracks generated on top of a bank, and (3) mass failure on flat or cambered planes [4–6]. Simultaneously, the respective roles of bank shapes [7], near-bank hydrology [8], positive and negative pore pressures [9], high confining water pressures [10], and riparian vegetation coverages [11], as well as bank materials were quantified in the modelization of riverbank collapse [12]. In addition, several bank stability models based on limit equilibrium were established to evaluate cohesive

**Citation:** Duan, G.; Liu, H.; Shao, D.; Yang, W.; Li, Z.; Wang, C.; Chang, S.; Ding, Y. Numerical Simulation of the Transportation of Cohesive Bank-Collapsed Materials in a Sharply Curved Channel. *Water* **2022**, *14*, 1147. https://doi.org/10.3390/w14071147

Academic Editors: Qiting Zuo, Xiangyi Ding, Guotao Cui and Wei Zhang

Received: 17 February 2022 Accepted: 29 March 2022 Published: 2 April 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

bank stability and predict collapse volumes [13,14]. These notable contributions present much benefit for predicting channel bend evolution processes, especially for the rivers with drastic riverbank collapse [15].

In a curved channel with drastic riverbank collapse, the evolution processes become more complicated under the interaction between collapsed materials and near-bank hydrology [16]. After bank collapse occurs, collapsed sediments that accumulated at the bank toe will change the original channel topography, which can affect the velocity distribution and provide a sediment source [17]. Previous studies noted that collapsed materials can reduce the near-bank shear stress [18,19], further increase the near bank resistance, and make the high velocity area shift away from the riverbank [20,21]. Yu et al. [22] found that the presence of a collapsed block can cause greater downstream bank retreat, while a smaller near-bank velocity can protect the bank against erosion occurring upstream of the block end. Xie et al. [23] noted that the average wall shear force between the collapse body and the toe decreases when the collapsed body is located upstream of the apex of bend. The opposite situation occurs when the collapse body is located downstream of the apex of bend. Based on these qualitative analyses, some mathematical models were developed by changing bank erosion parameters to reflect protection from collapsed materials [15,24–27]. Nevertheless, these studies mainly emphasized the influence of collapsed materials on flow distributions.

In fact, channel hydrology reacts on collapsed materials at the same time. One important manifestation is the accumulation and transportation of collapsed materials. Though some studies qualitatively describe the transportation based on flume experiments and theory of sediment movement [22,28,29], there has been no uniform approach for quantifying the transportation of collapsed materials directly until now because of its difficulty to observe either in natural rivers or laboratory flume experiments. Instead, many studies introduced assumptions about the transportation of collapsed materials when simulating channel evolution processes by coupling water–sediments equations, bed evolution equations, and bank stability models. Some studies considered collapsed cohesive sediment as wash loads that were carried away instantaneously, with none being accumulative [9]. Some studies classified the collapsed materials based on particle size. Particles finer than 0.062 mm were considered as wash load, while the coarser particles were considered as bed sediments that were distributed uniformly across the bed area between the bank toe and the boundary of the near bank sediment routing segment, a distance equal to twice the bank height. Zong et al. [30] considered 50% of collapsed materials as wash load; the others accumulated at the bank toe with triangular silting shapes. Duan et al. [31] proposed that the volume of accumulated sediments is decided by sediment carrying capacity and assumed that sediment accumulated at the sediment deposition angle. Then river evolution process was simulated through water–sediment and bed evolution equations. Though these assumptions were indirectly demonstrated by comparing simulated and measured results, the further fraction of reworked collapsed materials, transported either as bed load or as suspended load, was rarely involved. As this further fraction can influence channel downstream morphology significantly in natural alluvial rivers, it is essential to quantify the different transport patterns of collapsed materials.

In this study, scenarios with cohesive collapsed materials and non-cohesive sediments in a 180◦ U-bend channel were simulated by a numerical model established based in Delft3D to evaluate the transportation of the collapsed materials in a sharply curved channel and quantify the suspended and bed loads that are transformed from collapsed materials.

#### **2. Materials and Methods**

In this study, a three-dimensional mathematical model established in Delft3D was adopted [32], which is comprised of a flow model and a sediment transport model.

#### *2.1. Flow Model*

The three-dimensional bend flow model (k-ε model) in the σ coordinate system was applied to the Reynolds-averaged Navier–Stokes equations for incompressible free surface flow to obtain the flow velocities in the *ξ*, *η* and *σ* directions.

#### (1) Definition of *σ*-co-ordinate.

The *σ* co-ordinate system is defined as:

$$
\sigma = \frac{z - \zeta}{d + \overline{\zeta}} = \frac{z - \zeta}{H} \tag{1}
$$

where *z* is the vertical coordinate in physical space; *ζ* is the free surface elevation above the reference plane (at *z* = 0), m; d is the depth below the reference plane, m; and *H* is the total water depth, m.

$$H = d + \mathcal{J} \tag{2}$$

(2) Continuity equations.

In coordinate system *σ*, continuity equations can be transformed as follows:

$$\frac{\partial \zeta}{\partial t} + \frac{1}{\sqrt{G\_{\xi\xi}}\sqrt{G\_{\eta\eta}}} \frac{\partial \left[ (d+\zeta)L\sqrt{G\_{\eta\eta}} \right]}{\partial \xi} + \frac{1}{\sqrt{G\_{\xi\xi}}\sqrt{G\_{\eta\eta}}} \frac{\partial \left[ (d+\zeta)V\sqrt{G\_{\xi\xi}} \right]}{\partial \eta} = Q \tag{3}$$

where *t* is the time, s; (*Gξξ* ) 1/2 and (*Gηη*) 1/2 are the coefficients used to transform curvilinear to rectangular coordinates, m; *ξ* and *η* are horizontal curvilinear coordinates; *U* is depthaveraged velocity in *ξ* direction, m/s; *V* is depth-averaged velocity in *η* direction, m/s; and *Q* is the contributions per unit area due to the water discharge, m/s.

$$Q = H \int\_{-1}^{0} (q\_{in} - q\_{out}) d\sigma + P - E \tag{4}$$

where *qin* is the local source of water per unit of volume, 1/s; *qout* is the local sink of water per unit of volume, 1/s; *P* is precipitation, m/s; and *E* is evaporation, m/s.

(3) Momentum equations.

$$\frac{\partial\mu}{\partial t} + \frac{u}{\sqrt{\mathsf{C}\_{\mathsf{V}}^{2}}} \frac{\partial u}{\partial\xi} + \frac{v}{\sqrt{\mathsf{C}\_{\mathsf{V}}}} \frac{\partial u}{\partial\eta} + \frac{\omega}{d + \zeta} \frac{\partial u}{\partial\sigma} - \frac{v^{2}}{\sqrt{\mathsf{C}\_{\mathsf{V}} \varsigma \sqrt{\mathsf{C}\_{\mathsf{V}}}}} \frac{\partial \sqrt{\mathsf{C}\_{\mathsf{V}}}}{\partial\xi} + \frac{uv}{\sqrt{\mathsf{C}\_{\mathsf{V}} \varsigma \sqrt{\mathsf{C}\_{\mathsf{V}}}}} \frac{\partial \sqrt{\mathsf{C}\_{\mathsf{V}}}}{\partial\eta} - fv = -\frac{1}{\rho\_{0} \sqrt{\mathsf{C}\_{\mathsf{V}} \overline{\varsigma}}} \mathrm{P}\_{\overline{\xi}} + \mathrm{F}\_{\overline{\xi}} + \frac{1}{(d + \overline{\zeta})^{2}} \frac{\partial}{\partial\tau} (\mathrm{v} \eta \frac{\partial u}{\partial\tau}) + \mathrm{M}\_{\overline{\xi}} \quad \text{(5-2.32)}$$

$$\frac{\partial v}{\partial t} + \frac{u}{\sqrt{\text{C}\_{\text{f}}^{2}}} \frac{\partial v}{\partial \xi} + \frac{v}{\sqrt{\text{C}\_{\text{f}}}} \frac{\partial v}{\partial \eta} + \frac{u}{d + \xi} \frac{\partial v}{\partial r} + \frac{uv}{\sqrt{\text{C}\_{\text{f}} \xi \sqrt{\text{C}\_{\text{f}}}}} \frac{\partial \sqrt{\text{C}\_{\text{f}}}}{\partial \xi} - \frac{u^{2}}{\sqrt{\text{C}\_{\text{f}} \xi \sqrt{\text{C}\_{\text{f}}}}} \frac{\partial \sqrt{\text{C}\_{\text{f}}}}{\partial \eta} + fu = -\frac{1}{\rho\_{0} \sqrt{\text{C}\_{\text{f}} \eta}} P\_{\eta} + F\_{\eta} + \frac{1}{(d + \xi)^{2}} \frac{\partial}{\partial r} (v\nu \frac{\partial v}{\partial r}) + M\_{\eta} \tag{6}$$

where *u*, *v*, and *ω* are the velocities in the *ξ*, *η* and *σ* directions, respectively, m/s; *f* is the Coriolis parameter, 1/s; *ρ<sup>0</sup>* is the reference density of water, kg/m<sup>3</sup> ; *P<sup>ξ</sup>* and *P<sup>η</sup>* are the gradient hydrostatic pressures in the *ξ* and η directions, kg/(m<sup>2</sup> s 2 ); *F<sup>ξ</sup>* and *F<sup>η</sup>* are the turbulent momentum fluxes in the *ξ* and *η* directions, m/s<sup>2</sup> ; *M<sup>ξ</sup>* and *M<sup>η</sup>* are the sources or sinks of momentum in the *ξ* and *η* directions, m/s<sup>2</sup> ; and *ν<sup>V</sup>* is the vertical eddy viscosity, m2/s.

The vertical velocity *ω* in the adapting σ-coordinate system is computed from the continuity equation:

$$\frac{\partial \zeta}{\partial t} + \frac{1}{\sqrt{\mathsf{G}\_{\xi\xi}^{\mathsf{x}}} \sqrt{\mathsf{G}\_{\eta\eta}}} \frac{\partial \left[ (d+\zeta) u \sqrt{\mathsf{G}\_{\eta\eta}} \right]}{\partial \xi} + \frac{1}{\sqrt{\mathsf{G}\_{\xi\xi}^{\mathsf{x}}} \sqrt{\mathsf{G}\_{\eta\eta}}} \frac{\partial \left[ (d+\zeta) v \sqrt{\mathsf{G}\_{\xi\xi}} \right]}{\partial \eta} + \frac{\partial \omega}{\partial \sigma} = H(q\_{\mathrm{fin}} - q\_{\mathrm{out}}) \tag{7}$$

(4) Turbulent model: k-ε model.

$$\frac{\partial \mathbf{k}}{\partial t} + \frac{u}{\sqrt{\mathbf{G}\_{\xi\overline{\xi}}}} \frac{\partial \mathbf{k}}{\partial \overline{\xi}} + \frac{v}{\sqrt{\mathbf{G}\_{\eta\eta}}} \frac{\partial \mathbf{k}}{\partial \eta} + \frac{\omega}{d + \zeta} \frac{\partial \mathbf{k}}{\partial \sigma} = \frac{1}{(d + \zeta)^{2}} \frac{\partial}{\partial \sigma} \left( D\_{\mathbf{k}} \frac{\partial \mathbf{k}}{\partial \sigma} \right) + P\_{\mathbf{k}} + P\_{\mathbf{k}\overline{\nu}} + B\_{\mathbf{k}} - \varepsilon \tag{8}$$

$$\frac{\partial \varepsilon}{\partial t} + \frac{u}{\sqrt{G\_{\xi\xi}}} \frac{\partial \varepsilon}{\partial \xi} + \frac{v}{\sqrt{G\_{\eta\eta}}} \frac{\partial \varepsilon}{\partial \eta} + \frac{\omega}{d + \xi} \frac{\partial \varepsilon}{\partial \sigma} = \frac{1}{(d + \xi)^2} \frac{\partial}{\partial \sigma} \left( D\_{\varepsilon} \frac{\partial \varepsilon}{\partial \sigma} \right) + P\_{\varepsilon} + P\_{\varepsilon\varpi} + B\_{\varepsilon} - c\_{2\varepsilon} \frac{\varepsilon^2}{k} \tag{9}$$

where *k* is turbulent kinetic energy, m2/s<sup>2</sup> ; *P<sup>k</sup>* is production term in transport equation for turbulent kinetic energy, m2/s<sup>3</sup> ; *B<sup>k</sup>* is buoyancy flux term in transport equation for turbulent kinetic energy, m2/s<sup>3</sup> ; *ε* is dissipation in transport equation for turbulent kinetic energy, m2/s<sup>3</sup> ; *Pε* is production term in transport equation for the dissipation of turbulent kinetic energy, m2/s<sup>4</sup> ; *Bε* is buoyancy flux term in transport equation for the dissipation of turbulent kinetic energy, m2/s<sup>4</sup> ; *c1<sup>ε</sup>* , *c2<sup>ε</sup>* and *c3<sup>ε</sup>* are constant coefficients; *D<sup>k</sup>* and *D<sup>ε</sup>* are eddy diffusivities of *k* and *ε*.

$$D\_k = \frac{v\_{mol}}{\sigma\_{mol}} + \frac{v\_{3D}}{\sigma\_k} \tag{10}$$

$$D\_{\varepsilon} = \frac{v\_{3D}}{\sigma\_{\varepsilon}}\tag{11}$$

$$P\_k = 2\eta\_{3D} \left[ \frac{1}{2(d+\frac{\tau}{\xi})^2} \left\{ \left(\frac{\partial u}{\partial \sigma}\right)^2 + \left(\frac{\partial v}{\partial \sigma}\right)^2 \right\} + \left(\frac{1}{\sqrt{G\_{\eta\eta}}} \frac{\partial u}{\partial \xi}\right)^2 \right] + 2\eta\_{3D} \left[ \frac{1}{2} \left(\frac{1}{\sqrt{G\_{\eta\eta}}} \frac{\partial u}{\partial \eta} + \frac{1}{\sqrt{G\_{\xi\xi}}} \frac{\partial v}{\partial \eta} \right)^2 + \left(\frac{1}{\sqrt{G\_{\eta\eta}}} \frac{\partial v}{\partial \eta}\right)^2 \right] \tag{12}$$

$$v\_{\Im D} = c\_{\mu}^{\prime} L \sqrt{k} = c\_{\mu} \frac{k^2}{\varepsilon} \tag{13}$$

$$\mathfrak{c}\_{\mu} = \mathfrak{c}\_{D} \mathfrak{c}'\_{\mu} \tag{14}$$

where *L* is the mixing length, m; *c<sup>µ</sup>* is the calibration constant, *c<sup>µ</sup>* = 0.09; *c<sup>µ</sup>* 0 is a constant in Kolmogorov–Prandtl's eddy viscosity formulation; and *c<sup>D</sup>* is a constant relating the mixing length, turbulent kinetic energy and dissipation in the k-ε model.

$$B\_k = \frac{v\_{3D}}{\rho \sigma\_\rho} \frac{g}{H} \frac{\partial \rho}{\partial \sigma} \tag{15}$$

$$P\_{\varepsilon} = \mathfrak{c}\_{1\varepsilon} \frac{\varepsilon}{k} P\_k \tag{16}$$

$$B\_{\varepsilon} = c\_{1\varepsilon} \frac{\varepsilon}{k} (1 - c\_{3\varepsilon}) B\_k \tag{17}$$

$$L = c\_D \frac{k\sqrt{k}}{\varepsilon} \tag{18}$$

$$c\_{1\varepsilon} = 1.44\tag{19}$$

$$c\_{2\varepsilon} = 1.92\tag{20}$$

$$
\mathscr{L}\_{\mathfrak{A}} = \mathbf{0} \tag{21}
$$

$$
\varepsilon|\_{\sigma=-1} = \frac{\mu\_{\ast b}^3}{\kappa z\_0} \tag{22}
$$

$$\varepsilon|\_{\sigma=0} = \frac{u\_{\ast b}^3}{\frac{1}{2}\kappa \triangle z\_s} \tag{23}$$

*2.2. Sediment Transport Equation*

(1) Suspended sediment transport equation

$$\frac{\partial \mathbf{c}^{(l)}}{\partial t} + \frac{\partial \mathbf{u} \mathbf{c}^{(l)}}{\partial \mathbf{x}} + \frac{\partial \mathbf{c} \mathbf{c}^{(l)}}{\partial y} + \frac{\partial \left(\mathbf{w} - \mathbf{w}\_s^{(l)}\right) \mathbf{c}^{(l)}}{\partial z} - \frac{\partial}{\partial \mathbf{x}} \left(\varepsilon\_{\mathbf{s},x}^{(l)} \frac{\partial \mathbf{c}^{(l)}}{\partial \mathbf{x}}\right) - \frac{\partial}{\partial y} \left(\varepsilon\_{\mathbf{s},y}^{(l)} \frac{\partial \mathbf{c}^{(l)}}{\partial y}\right) - \frac{\partial}{\partial z} \left(\varepsilon\_{\mathbf{s},z}^{(l)} \frac{\partial \mathbf{c}^{(l)}}{\partial z}\right) = 0\tag{24}$$

where *c (l)* is the mass concentration of the sediment fraction *(l)*, kg/m<sup>3</sup> ; *u*, *v,* and *w* are the flow velocities in the *x*, *y,* and *z* directions, m/s; *εs,x (l)* , *εs,x (l)* and *εs,x (l)* are the eddy diffusivities of the sediment fraction *(l)*, m2/s; and *w<sup>s</sup> (l)* is the sediment settling velocity of the sediment fraction *(l)*, m/s.

(2) Bed load transport equation

The reference concentration is calculated as:

$$c\_a^{(l)} = 0.015 \rho\_s^{(l)} \frac{D\_{50}^{(l)} \left(T\_a^{(l)}\right)^{1.5}}{a \left(D\_\*^{(l)}\right)^{0.3}} \tag{25}$$

where *c<sup>a</sup> (l)* is the mass concentration at reference height a, kg/m<sup>3</sup> ; *D\* (l)* is the nondimensional particle diameter; and *T<sup>a</sup> (l)* is the nondimensional bed shear stress.

$$D\_\*^{(l)} = D\_{50}^{(l)} \left[ \frac{\left(s^{(l)} - 1\right)g}{v^2} \right]^{1/3} \tag{26}$$

$$T\_a^{(l)} = \frac{\left(\mu\_c^{(l)}\tau\_{b,cw} + \mu\_w^{(l)}\tau\_{b,w}\right) - \tau\_{cr}^{(l)}}{\tau\_{cr}^{(l)}}\tag{27}$$

$$|S\_{\theta}| = 0.006 \rho\_{s} w\_{s} D\_{50}^{(l)} M^{0.5} M\_{\text{e}}^{0.7} \tag{28}$$

where *ρ<sup>s</sup>* is the density of sediment, kg/m<sup>3</sup> ; *w<sup>s</sup>* is the sediment setting velocity, m/s; *M* is the sediment mobility number due to waves and currents; and *M<sup>e</sup>* is the excess sediment mobility number.

#### *2.3. Model Verification*

The numerical model was verified by simulating water flume experiments that studied the interaction between the riverbank and bed sediments, as described in the literature [33]. The flume experiments were constructed at the State Key Laboratory of Water Resources and Hydropower Engineering Sciences at Wuhan University to simulate riverbank collapse in Dengkou Reach of the Yellow River. The slope of flume was 1‰. The flow conditions were designed as: (1) discharge 30 L/s; (2) water depth 0.14 m; (3) water level 0.24 m. Water flume sizes, section settings, and sediment distributions were organized as shown in Figure 1. The erodible sediments were placed only on the shaded part of the water flume. The particle sizes of all the bed and bank sediments that collected from Dengkou Reach of the Yellow River are listed in Figure 2. Additionally, Table 1 lists the physical properties of cohesive sediments. The flow velocity was measured by an acoustic Doppler velocity (ADV) profiler; by mounting the probe to a rigid scaffold, any possible vibration of the probe due to flow action was minimized. When the experiment was completed, the

*Water* **2022**, *14*, x FOR PEER REVIEW 6 of 19

*Water* **2022**, *14*, x FOR PEER REVIEW 6 of 19

bed level was measured by a transparent mesh plate at typical sections after the water was drained.

**Figure 1.** Schematic diagram of the flume experiment. **Figure 1.** Schematic diagram of the flume experiment. 6.75 37.2 17.2 21.5

**Figure 2.** Particle distribution of cohesive and non-cohesive sediments. **Figure 2.** Particle distribution of cohesive and non-cohesive sediments.



A basin of 188 × 24 computational cells, where each was approximately 0.01 m2 at 0.2 m long and 0.05 m wide, was positioned with an initial bed slope of 1‰ in the rectangular flume along the flow direction, as shown in Figure 3. A basin of 188 <sup>×</sup> 24 computational cells, where each was approximately 0.01 m<sup>2</sup> at 0.2 m long and 0.05 m wide, was positioned with an initial bed slope of 1‰ in the rectangular flume along the flow direction, as shown in Figure 3.

The inlet and outlet boundary conditions shown in Figure 3 were set to be the discharge (30 L/s) and water level (0.24 m), respectively. All cohesive and non-cohesive grains had a density of 2650 kg/m<sup>3</sup> , as the dry densities were set to 500 kg/m<sup>3</sup> and 1600 kg/m<sup>3</sup> , respectively, based on measurements. The settling velocities of cohesive sediments that

ranged from 0.04 to 0.46 mm/s in this research were calculated using the formula derived by Dou [34]. For cohesive sediments, the critical bed shear stress for erosion ranged from 0.2 to 1.8 N/m<sup>2</sup> ; here we calculated using the formula in the literature [35,36]. The critical bed shear stress for deposition can be obtained using the formula in the literature [37]. The background horizontal viscosity and diffusivity values were set to 0.0005 and 0.0001 m2/s, respectively. *Water* **2022**, *14*, x FOR PEER REVIEW 7 of 19

*Water* **2022**, *14*, x FOR PEER REVIEW 7 of 19

**Figure 3.** Morphologic grid in calculate areas. **Figure 3.** Morphologic grid in calculate areas. A bed stratigraphy model containing 5 layers was used to track the evolution of the

The inlet and outlet boundary conditions shown in Figure 3 were set to be the discharge (30 L/s) and water level (0.24 m), respectively. All cohesive and non-cohesive grains had a density of 2650 kg/m3, as the dry densities were set to 500 kg/m3 and 1600 kg/m3, respectively, based on measurements. The settling velocities of cohesive sediments that ranged from 0.04 to 0.46 mm/s in this research were calculated using the A bed stratigraphy model containing 5 layers was used to track the evolution of the bed sediment. A time step of 0.03 s was adopted to obey all stability criteria, and the simulation time was 3 h when riverbank collapse did not occur. To avoid numerical instabilities caused by supercritical flow in shallow areas, a grid cell was considered dry if its depth was shallower than 0.01 cm. Boundary fitted grids were used in this model so that all erosion and deposition fluxes were applied to the bottom cell face. bed sediment. A time step of 0.03 s was adopted to obey all stability criteria, and the simulation time was 3 h when riverbank collapse did not occur. To avoid numerical instabilities caused by supercritical flow in shallow areas, a grid cell was considered dry if its depth was shallower than 0.01 cm. Boundary fitted grids were used in this model so that all erosion and deposition fluxes were applied to the bottom cell face. (1) Velocity verification

formula derived by Dou [34]. For cohesive sediments, the critical bed shear stress for (1) Velocity verification Because the main object of this study was to quantify the transformation of collapsed

erosion ranged from 0.2 to 1.8 N/m2; here we calculated using the formula in the literature [35,36]. The critical bed shear stress for deposition can be obtained using the formula in the literature [37]. The background horizontal viscosity and diffusivity values were set to 0.0005 and 0.0001 m2/s, respectively. A bed stratigraphy model containing 5 layers was used to track the evolution of the bed sediment. A time step of 0.03 s was adopted to obey all stability criteria, and the simulation time was 3 h when riverbank collapse did not occur. To avoid numerical instabilities caused by supercritical flow in shallow areas, a grid cell was considered dry if its depth was shallower than 0.01 cm. Boundary fitted grids were used in this model so that all erosion and deposition fluxes were applied to the bottom cell face. Because the main object of this study was to quantify the transformation of collapsed materials, the depth-averaged velocity played a more important role in longitudinal direction. Thus, the depth-averaged velocities from numerical simulation and flume experiments were selected to verify the adopted model. The average velocities in Sections 1, 3–7 were adopted to verify the simulated velocities with the measured velocities in the water flume experiments, as shown in Figure 4. Figure 5 presents the comparison between simulated and measured velocities at the same location in every selected section. The simulated velocities were moderately consistent with the experimental velocities, supporting the validity of the flow model. materials, the depth-averaged velocity played a more important role in longitudinal direction. Thus, the depth-averaged velocities from numerical simulation and flume experiments were selected to verify the adopted model. The average velocities in Sections 1, 3–7 were adopted to verify the simulated velocities with the measured velocities in the water flume experiments, as shown in Figure 4. Figure 5 presents the comparison between simulated and measured velocities at the same location in every selected section. The simulated velocities were moderately consistent with the experimental velocities, supporting the validity of the flow model.

**Figure 4.** Comparison of simulated velocities and measured velocities.

**Figure 4.** Comparison of simulated velocities and measured velocities.

**Figure 5.** Comparison of the simulated and measured velocities in different sections. (**a**) Section 1; (**b**) Section 3; (**c**) Section 4; (**d**) Section 5; (**e**) Section 6; (**f**) Section 7. **Figure 5.** Comparison of the simulated and measured velocities in different sections. (**a**) Section 1; (**b**) Section 3; (**c**) Section 4; (**d**) Section 5; (**e**) Section 6; (**f**) Section 7.

#### (2) Bed level verification

enough to investigate the evolution of flume bottom variation.

(2) Bed level verification Bed sediment transportation was verified by comparing bed levels between simulated and measured results in Sections 1, 3, 4, 5, 6, and 7 separately, as shown in Figure 6. Due to the existence of the transverse effect, the water at the surface in the curve channel will flow towards the concave bank, whereas the water at the bottom will flow in the opposite direction, which will carry sediments. Thus, the bed levels near the convex bank were higher than those near the concave bank in both simulated and experimental sce-Bed sediment transportation was verified by comparing bed levels between simulated and measured results in Sections 1, 3, 4, 5, 6, and 7 separately, as shown in Figure 6. Due to the existence of the transverse effect, the water at the surface in the curve channel will flow towards the concave bank, whereas the water at the bottom will flow in the opposite direction, which will carry sediments. Thus, the bed levels near the convex bank were higher than those near the concave bank in both simulated and experimental scenarios, as depicted in Figure 6. It can also be seen that the adopted model is accurate enough to investigate the evolution of flume bottom variation.

narios, as depicted in Figure 6. It can also be seen that the adopted model is accurate

**Figure 6.** Comparison of the simulated and measured bed levels in different sections. (**a**) Section 1; (**b**) Section 3; (**c**) Section 4; (**d**) Section 5; (**e**) Section 6; (**f**) Section 7. **Figure 6.** Comparison of the simulated and measured bed levels in different sections. (**a**) Section 1; (**b**) Section 3; (**c**) Section 4; (**d**) Section 5; (**e**) Section 6; (**f**) Section 7.

#### *2.4. Simulated Scenario*

*2.4. Simulated Scenario*  2.4.1. Location of Collapsed Materials

2.4.1. Location of Collapsed Materials Two scenarios listed in Table 2 were simulated by using the adopted numerical model to determine the location of collapsed materials under different flow conditions (30 L/s and 55 L/s). The water shear stress distribution of the curved channel, as the main driving force on riverbank collapse [28], was simulated without collapsed materials to determine the probable location of collapsed material accumulation. As bank collapse Two scenarios listed in Table 2 were simulated by using the adopted numerical model to determine the location of collapsed materials under different flow conditions (30 L/s and 55 L/s). The water shear stress distribution of the curved channel, as the main driving force on riverbank collapse [28], was simulated without collapsed materials to determine the probable location of collapsed material accumulation. As bank collapse occurs in the downstream of the bend exit [38], the collapsed materials were set to accumulate at the location in situ, circled by a dotted line in Figure 7, where the water shear stress was greater.

cumulate at the location in situ, circled by a dotted line in Figure 7, where the water shear **Table 2.** Simulated scenarios settings.


1 30 0.24 0.03 8 2 55 0.24 0.03 8

occurs in the downstream of the bend exit [38], the collapsed materials were set to ac-

**Figure 7.** Water shear stress distribution under different flow charge conditions. (**a**) Q = 30 L/s; (**b**) Q= 55 L/s. **Figure 7.** Water shear stress distribution under different flow charge conditions. (**a**) Q = 30 L/s; (**b**) Q = 55 L/s.

### 2.4.2. Composition of Collapsed Materials

2.4.2. Composition of Collapsed Materials In natural rivers, the shape of collapsed materials accumulated at the bank toe is very difficult to measure and predict due to its irregularity. Based on the existing research [30,31], the shape of collapsed materials is assumed to be right-angled trapezoidal in the cross-sectional direction, as shown in Figure 8a. The length of the collapsed materials is determined by the area in the dotted line (Figure 7), and the total thickness of the collapsed materials was 50 mm (Figure 8a). Since the critical size of coarse and fine sediment in the Yellow River was 50μm [39], two typical particles were selected to represent coarser and fine sediment, i.e., 37 μm (S1) and 74 μm (S2), based on the particle size distribution of the riverbank (Figure 2). It can also be obtained that the ratio of particles coarser than 50μm to finer than 50μm approximately equals to 3; thus, the ratio of In natural rivers, the shape of collapsed materials accumulated at the bank toe is very difficult to measure and predict due to its irregularity. Based on the existing research [30,31], the shape of collapsed materials is assumed to be right-angled trapezoidal in the crosssectional direction, as shown in Figure 8a. The length of the collapsed materials is determined by the area in the dotted line (Figure 7), and the total thickness of the collapsed materials was 50 mm (Figure 8a). Since the critical size of coarse and fine sediment in the Yellow River was 50 µm [39], two typical particles were selected to represent coarser and fine sediment, i.e., 37 µm (S1) and 74 µm (S2), based on the particle size distribution of the riverbank (Figure 2). It can also be obtained that the ratio of particles coarser than 50µm to finer than 50µm approximately equals to 3; thus, the ratio of VS1 to VS2 equals 3, where VS1 and VS2 represent the volumes of S1 and S2, respectively. Bed sediments were made of one typical particle size, 490 µm (S3), and the thickness was 100 mm, as shown in Figure 8.

#### VS1 to VS2 equals 3, where VS1 and VS2 represent the volumes of S1 and S2, respectively. Bed sediments were made of one typical particle size, 490 μm (S3), and the thickness was 2.4.3. Scenario Settings in Simulation for Transportation of Collapsed Materials

100 mm, as shown in Figure 8. The riverbank was fixed except for the collapsed materials and the bed sediments, as the main purpose of this study was to evaluate the quantities of collapsed materials transforming into suspended loads and bed loads. Based on the determined location and composition of collapsed materials, the grid and bed levels in the numerical model are shown in Figure 8b. The simulated scenario settings in Table 2 were used for boundary conditions. New section labels named based on grid numbers along the flow direction were also used in the following article to analyze the simulated results more conveniently. To ensure that the flow reached a steady state, river geomorphology was incorporated after 30 min. The simulated results based on the verified numerical model were analyzed in Section 3.

**Figure 8.** Collapsed materials and section settings in simulated scenarios. (**a**) Schematic diagram; (**b**) grid and bed level in numerical model. **Figure 8.** Collapsed materials and section settings in simulated scenarios. (**a**) Schematic diagram; (**b**) grid and bed level in numerical model.

#### **3. Results**

#### 2.4.3. Scenario Settings in Simulation for Transportation of Collapsed Materials *3.1. Curved Channel Evolution Process*

The riverbank was fixed except for the collapsed materials and the bed sediments, as the main purpose of this study was to evaluate the quantities of collapsed materials transforming into suspended loads and bed loads. Based on the determined location and composition of collapsed materials, the grid and bed levels in the numerical model are shown in Figure 8b. The simulated scenario settings in Table 2 were used for boundary conditions. New section labels named based on grid numbers along the flow direction were also used in the following article to analyze the simulated results more conveniently. To ensure that the flow reached a steady state, river geomorphology was incorporated after 30 min. The simulated results based on the verified numerical model were analyzed in Section 3. **3. Results**  *3.1. Curved Channel Evolution Process*  The variation in the riverbed elevation is an important performance factor in the curved channel evolution process. Figure 9 presents the bed levels at different times. As mentioned in Section 2.4, the river topography module began to run after 30 min of hydrodynamic pre-operation with a flow charge of 30 L/s or 55 L/s. There were three obvious characteristics in the curved channel evolution process. First, the bed level at the bend inlet decreased because sediments were washed downstream. Second, the bed level of the convex bank was higher than that of the concave bank over time due to the existence of cross circulation and surface pressure differences. Third, the bed level changed more dramatically when the flow charge was larger. Taking bed levels at 02:00:00 as an example, when the flow charge was 55 L/s, the bed levels of the bend inlet and concave bank were obviously lower, while the bed levels of the convex bank downstream of the bend outlet were obviously higher. This is attributed to the larger water shear stress and sediment carrying capacity of a larger flow charge with the same constant water level.

The variation in the riverbed elevation is an important performance factor in the curved channel evolution process. Figure 9 presents the bed levels at different times. As mentioned in Section 2.4, the river topography module began to run after 30 min of hydrodynamic pre-operation with a flow charge of 30 L/s or 55 L/s. There were three obvious characteristics in the curved channel evolution process. First, the bed level at the bend inlet decreased because sediments were washed downstream. Second, the bed level of the convex bank was higher than that of the concave bank over time due to the existence of cross circulation and surface pressure differences. Third, the bed level changed more dramatically when the flow charge was larger. Taking bed levels at 02:00:00 as an example, when the flow charge was 55 L/s, the bed levels of the bend inlet and concave Figure 10 describes the bed level changes in typical sections (M114, 118 and 121) in collapsed materials reached at different moments. Whether the flow charge was 30 L/s or 55 L/s, the bed levels of the right bank, where collapsed materials accumulated, decreased as collapsed materials were eroded by water flow. For the left bank, the bed level changes were more complicated. When the flow charge was 30 L/s, the bed levels of the left bank increased as time progressed. When the flow charge was 55 L/s, the bed levels of the left bank increased in the first 2 h and then decreased because a larger flow charge was more conducive to sediment transportation [40]. After 2 h, there were fewer sediments upstream, and the sediments that were previously deposited on the left bank were transported downstream under the action of water flow when the flow charge was 55 L/s.

bank were obviously lower, while the bed levels of the convex bank downstream of the bend outlet were obviously higher. This is attributed to the larger water shear stress and sediment carrying capacity of a larger flow charge with the same constant water level.

**Figure 9.** Bed level of the curved channel at different times. (**a**) Q = 30 L/s, T = 00:00:00; (**b**) Q = 55 L/s, T = 00:00:00; (**c**) Q = 30 L/s, T = 02:00:00; (**d**) Q = 55 L/s, T = 02:00:00; (**e**) Q = 30 L/s, T = 04:00:00; (**f**) Q = 55 L/s, T = 04:00:00; (**g**) Q = 30 L/s, T = 06:00:00; (**h**) Q = 55 L/s, T = 06:00:00; (**i**) Q = 30 L/s, T = 08:00:00; (**j**) Q = 55 L/s, T = 08:00:00. **Figure 9.** Bed level of the curved channel at different times. (**a**) Q = 30 L/s, T = 00:00:00; (**b**) Q = 55 L/s, T = 00:00:00; (**c**) Q = 30 L/s, T = 02:00:00; (**d**) Q = 55 L/s, T = 02:00:00; (**e**) Q = 30 L/s, T = 04:00:00; (**f**) Q = 55 L/s, T = 04:00:00; (**g**) Q = 30 L/s, T = 06:00:00; (**h**) Q = 55 L/s, T = 06:00:00; (**i**) Q = 30 L/s, T = 08:00:00; (**j**) Q = 55 L/s, T = 08:00:00.

left bank increased as time progressed. When the flow charge was 55 L/s, the bed levels of the left bank increased in the first 2 h and then decreased because a larger flow charge was more conducive to sediment transportation [40]. After 2 h, there were fewer sedi-

Figure 10 describes the bed level changes in typical sections (M114, 118 and 121) in collapsed materials reached at different moments. Whether the flow charge was 30 L/s or 55 L/s, the bed levels of the right bank, where collapsed materials accumulated, decreased as collapsed materials were eroded by water flow. For the left bank, the bed level

ments upstream, and the sediments that were previously deposited on the left bank were transported downstream under the action of water flow when the flow charge was 55 L/s.

**Figure 10.** Bed levels of typical sections in collapsed materials reach (**a**) M114, Q = 30 L/s; (**b**) M114, Q = 55 L/s; (**c**) M118, Q = 30 L/s; (**d**) M118, Q = 55 L/s; (**e**) M121, Q = 30 L/s; (**f**) M121, Q = 55 L/s. **Figure 10.** Bed levels of typical sections in collapsed materials reach (**a**) M114, Q = 30 L/s; (**b**) M114, Q = 55 L/s; (**c**) M118, Q = 30 L/s; (**d**) M118, Q = 55 L/s; (**e**) M121, Q = 30 L/s; (**f**) M121, Q = 55 L/s.

#### *3.2. Quantities of the Collapsed Material Transformation and Transportation 3.2. Quantities of the Collapsed Material Transformation and Transportation* 3.2.1. Quantities of the Eroded Collapsed Materials

3.2.1. Quantities of the Eroded Collapsed Materials The quantities of the eroded sediment fractions can be obtained by arranging the remaining sediments on the riverbed at different moments. Figure 11 shows the remaining quantities of S2 on the riverbed at the beginning and end of the simulation when the flow charge was 30 L/s. The eroded quantity of S2 can be obtained by calculating the difference. Tables 3 and 4 list the erosion quantities of different sediment fractions at different moments. When the flow charge was 55 L/s, the erosion quantities of both sediment fractions were larger because of the larger water shear stress. Whether the flow charge was 30 L/s or 55 L/s, the erosion quantity of S2 was larger than that of S1. The main driving force was the water shear stress in the cohesive sediment erosion process, while the resistance provided the cohesive force, friction, gravity of particles, and electrochemical force effects. When the particle size was finer, the electrochemical force was The quantities of the eroded sediment fractions can be obtained by arranging the remaining sediments on the riverbed at different moments. Figure 11 shows the remaining quantities of S2 on the riverbed at the beginning and end of the simulation when the flow charge was 30 L/s. The eroded quantity of S2 can be obtained by calculating the difference. Tables 3 and 4 list the erosion quantities of different sediment fractions at different moments. When the flow charge was 55 L/s, the erosion quantities of both sediment fractions were larger because of the larger water shear stress. Whether the flow charge was 30 L/s or 55 L/s, the erosion quantity of S2 was larger than that of S1. The main driving force was the water shear stress in the cohesive sediment erosion process, while the resistance provided the cohesive force, friction, gravity of particles, and electrochemical force effects. When the particle size was finer, the electrochemical force was greater and played a more important role than other resistances [41]. *Water* **2022**, *14*, x FOR PEER REVIEW 14 of 19

**Figure 11.** Remaining S2 on riverbed. (**a**) Q = 30 L/s, T = 00:00:00; (**b**) Q = 30 L/s, T = 08:00:00. materials were transported as both suspended and bed loads after erosion. When the **Figure 11.** Remaining S2 on riverbed. (**a**) Q = 30 L/s, T = 00:00:00; (**b**) Q = 30 L/s, T = 08:00:00.

 **1 h 2 h 4 h 6 h 8 h**  S1 (kg) 0.013 0.036 0.064 0.086 0.107 S2 (kg) 0.559 1.268 2.306 3.085 5.670 Total (kg) 0.561 1.274 2.320 3.111 5.777

 **1 h 2 h 4 h 6 h 8 h**  S1 (kg) 0.040 0.164 0.392 0.550 0.655 S2 (kg) 4.068 6.769 7.428 7.574 7.599 Total (kg) 4.108 6.933 7.820 8.124 8.254

After the collapsed materials were eroded, they were transported downstream as suspended loads or bed loads. Several typical sections (M122, M124, M126, M128, and M130) downstream of the collapsed reach were adopted to calculate the quantities, as listed in Tables 5 and 6. When the flow charge was 30 L/s, the quantities of both suspended and bed loads flowing through sections decreased downstream. When the flow charge was 55 L/s, the quantities of the suspended loads flowing through all sections

 **M122 M124 M126 M128 M130** 

 **M122 M124 M126 M128 M130** 

**Table 5.** Quantities of the suspended and bed loads through typical cross sections (30 L/s).

**Table 6.** Quantities of the suspended and bed loads through typical cross sections (55 L/s).

Suspended loads (kg) 5.224 2.587 0.097 0.0895 0.0825 Bed loads (kg) 0.553 0.192 0.013 0.0005 0.0005 Total (kg) 5.777 2.779 0.110 0.090 0.083

Suspended loads (kg) 2.542 2.361 2.310 2.276 2.182 Bed loads (kg) 2.423 2.537 2.490 2.378 2.327 Total (kg) 4.965 4.898 4.800 4.654 4.509

3.2.3. Percentage of Sediment Fractions in Suspended and Bed Loads Transforming from

Tables 7 and 8 list the percentage of different sediment fractions in suspended loads and bed loads that were transformed from the collapsed materials. No bed load was obtained from S1 in any section and the flow charge was either 30 L/s or 55 L/s—i.e., all the collapsed materials at S1 were transported as suspended loads, while the S2 collapsed

**Table 4.** Erosion quantities of the collapsed materials (55 L/s).

3.2.2. The Quantities of the Suspended and Bed Loads

were almost the same as that of the bed load.

the Collapsed Materials


**Table 3.** Erosion quantities of the collapsed materials (30 L/s).

**Table 4.** Erosion quantities of the collapsed materials (55 L/s).


3.2.2. The Quantities of the Suspended and Bed Loads

After the collapsed materials were eroded, they were transported downstream as suspended loads or bed loads. Several typical sections (M122, M124, M126, M128, and M130) downstream of the collapsed reach were adopted to calculate the quantities, as listed in Tables 5 and 6. When the flow charge was 30 L/s, the quantities of both suspended and bed loads flowing through sections decreased downstream. When the flow charge was 55 L/s, the quantities of the suspended loads flowing through all sections were almost the same as that of the bed load.

**Table 5.** Quantities of the suspended and bed loads through typical cross sections (30 L/s).


**Table 6.** Quantities of the suspended and bed loads through typical cross sections (55 L/s).


3.2.3. Percentage of Sediment Fractions in Suspended and Bed Loads Transforming from the Collapsed Materials

Tables 7 and 8 list the percentage of different sediment fractions in suspended loads and bed loads that were transformed from the collapsed materials. No bed load was obtained from S1 in any section and the flow charge was either 30 L/s or 55 L/s—i.e., all the collapsed materials at S1 were transported as suspended loads, while the S2 collapsed materials were transported as both suspended and bed loads after erosion. When the flow charge was 30 L/s, the percentage of S1 suspended loads ranged from 2.05% to 11.64%, whereas that of S2 ranged from 88.36% to 97.95%. When the flow charge was 55 L/s, the percentage of S1 ranged from 5.15% to 6.46%, whereas that of S2 ranged from 93.54% to 94.85%.


**Table 7.** Quantities and percentage of sediment fractions flowing through cross sections as suspended and bed loads (30 L/s).

**Table 8.** Quantities and percentage of sediment fractions flowing through cross sections as suspended and bed loads (55 L/s).


Tables 9 and 10 list the percentage of S2 collapsed materials transforming into suspended and bed loads in typical sections. When the flow charge was 30 L/s, approximately 88–99.8% of S2 collapsed materials was transported as suspended loads, and only approximately 0.2–12% was transported as bed loads. When the flow charge was 55 L/s, approximately 47–50% was transported as suspended loads, and approximately 50–53% was transported as bed loads.

**Table 9.** Percentage of S2 collapsed materials transforming into suspended and bed loads (30 L/s).


**Table 10.** Percentage of S2 collapsed materials transforming into suspended and bed loads (55 L/s).


#### **4. Discussion**

Quantifying the transformation and transportation of collapsed materials in a curved channel is challenging because there are many influencing factors, such as flow velocity, water level, topography, and sediment characteristics. In addition, the suspended and bed loads will also transform mutually in the process. Nevertheless, under specific flow conditions, when sediments transportation is at equilibrium state, the quantity of the suspended and bed loads across typical sections would barely change [42]. Based on this, the results listed in Tables 7 and 8 were reasonable.

This study could be considered as a classic attempt to quantify the transportation of cohesive collapsed materials with non-cohesive riverbed, since there is rarely literature on this problem. The results not only demonstrate the existing theory of sediment transportation but also provide the ratio of suspended and bed loads that transformed from collapsed materials under certain flow conditions. At the same time, there are also several assumptions based on previous literature, such as the shape of collapsed material accumulation and the composition of collapsed materials. These reasonable assumptions would bring benefits in the numerical simulation of sediment transportation incontestably. However, they can also make the experimental results reified and difficult to apply to general phenomena.

Finally, as a key factor influencing sediment characteristics, particle size distribution plays an important role in the transformation and transportation of collapsed materials. Under certain flow conditions, sediment particle transport patterns differ. In this study, finer particles (S1) were only transported as suspended loads, while coarser particles (S2) were transported as both suspended and bed loads. For the same particle size distribution, the transport patterns were also different under different flow conditions [43]. Thus, it is essential to consider flow conditions and particle size distributions comprehensively when predicting natural river evolution processes.

#### **5. Conclusions**

Numerical studies were conducted to investigate the transformation and transportation of cohesive collapsed materials in a 180◦ sharply bent flume. The findings can be summarized as follows:


Furthermore, additional outcomes must be highlighted in future studies. First, as the particle size distribution of riverbanks and beds in nature ranges widely, more sediment fractions of collapsed materials and riverbeds should be added in the following simulations. Second, since the quantities of riverbank collapse events largely vary in natural rivers, different quantities of collapsed materials should be considered in future studies.

**Author Contributions:** Conceptualization, G.D. and H.L.; funding acquisition, H.L.; investigation, C.W. and S.C.; methodology, H.L., Z.L. and Y.D.; project administration, H.L.; validation, D.S. and W.Y.; writing—original draft, G.D. and H.L.; writing—review and editing, G.D., H.L., D.S., W.Y. and Y.D. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was supported by the Major Scientific and Technological Innovation Projects in Shandong Province (2021CXGC011201) and the Open Research Fund of Key Laboratory of Hydro-Sediment Science and River Training, the Ministry of Water Resources, China Institute of Water Resources and Hydropower Research (Grant number: IWHR-JH-2020-A-03).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

