*2.3. Cohesive Sediment Transport Equations*

The cohesive sediment transport equation is similar to the noncohesive sediment transport equation [48], except the method to calculate the net exchange rate (*Db* − *Eb*) [14].

$$\frac{\partial(\text{HC})}{\partial t} + \frac{1}{J} \left[ \frac{\partial(h\_1 \text{HC} I)}{\partial \xi} + \frac{\partial(h\_1 \text{HC} V)}{\partial \eta} \right] - \frac{H\varepsilon}{J} \left[ \frac{\partial}{\partial \xi} \left( \frac{h\_2}{h\_1} \frac{\partial \text{C}}{\partial \xi} \right) + \frac{\partial}{\partial \eta} \left( \frac{h\_1}{h\_2} \frac{\partial \text{C}}{\partial \xi} \right) \right] + D\_b - E\_b = 0 \tag{13}$$

where *<sup>C</sup>* <sup>=</sup> the sediment concentration, kg·m<sup>−</sup>3; *Db* and *Eb* <sup>=</sup> the erosion and deposition rate respectively, kg·m−2·s<sup>−</sup>1, which are calculated [14] as follows

$$D\_{\mathfrak{b}} = \alpha \omega\_{\mathfrak{s}} \mathbb{C} \tag{14}$$

where <sup>α</sup> <sup>=</sup> deposition coefficient calculated by Equation (15); <sup>ω</sup>*<sup>s</sup>* <sup>=</sup> settling velocity, m·s<sup>−</sup>1.

$$\alpha = \begin{cases} 1 - \frac{\tau\_b}{\tau\_{cd}}, \tau\_b \le \tau\_{cd} \\ 0, \tau\_b > \tau\_{cd} \end{cases} \tag{15}$$

where τ*<sup>b</sup>* = the bed shear stress, Pa; τ*cd* = the critical shear stress for deposition, Pa.

$$E\_b = M \left(\frac{\tau\_b}{\tau\_{cc}}\right)^n \tau\_b > \tau\_{cc} \tag{16}$$

where *<sup>n</sup>* <sup>=</sup> an empirical coefficient; *<sup>M</sup>* <sup>=</sup> the erosion coefficient, kg·m−2·s<sup>−</sup>1; and <sup>τ</sup>*ce* <sup>=</sup> the critical shear stress for erosion, Pa.

Most of the model parameters used for cohesive sediment calculation (Table 2) have been calibrated and validated by the sedimentation process of the Three Gorges Reservoir on Yangtze River [48], where the study area of this paper is located. A larger value of settling velocity is chosen from measurements by Li et al. [11,12] because only the medium diameter of the sediment is considered in this study.

**Table 2.** Model parameters used for cohesive sediment calculation.


The morphological evolution due to cohesive sediment transport is calculated by the net sediment exchange rate (*Db* − *Eb*), in the same way as noncohesive suspended sediment calculation does. The flow and sediment modules are solved in an uncoupled way. Details of the numerical method can be found in Wang et al. [42]. Central difference explicit scheme is applied to Equation (5), and

Equations (12) and (13) are solved by QUICK (Quadratic Upstream Interpolation for Convective Kinematics) finite difference scheme.

#### **3. Study Case**

The Hunghuacheng reach (HHC, Figure 2), located 364 km upstream from the Three Gorges Project (TGP), is approximately 13 km long, consisting two sharply curved bends with a center bar named "Huanghuacheng" splitting the reach into two branches. It belongs to the back water zone of TGP. The large mean annual discharge (32,000 m3·s<sup>−</sup>1) makes it a mega river reach [49]. Measurements of bed topographies and bed material size are taken at nine cross-sections from S201 to S209 twice each year. Due to huge amount of cohesive sediment siltation, the left branch of this reach has been blocked in September 2010 [13]. Secondary flow models are applied to this reach because the secondary flow caused by the upstream bend of this reach plays an important role in channel morphodynamics [50,51]. Also, it has been shown that similar models perform well in confluence [38] and braided rivers [25], which justify the application of these models in this reach.

**Figure 2.** Planform geometry, bed elevation (*Z*0) on March 2012 and nine cross-sections measured in HHC reach (S209 and S201 are the inlet and outlet boundaries, respectively; incoming flow discharge and sediment concentration used as inlet boundaries are interpolated from Qingxichang and Wanxian gauging station, located upstream 476.46 km and 291.38 km from TGP, respectively; and the outlet boundary applies the water stage measured at Shibaozhai station, located upstream 341.35 km from the TGP).

The year 2012 is selected to study the secondary flow effects on bed morphology variation in this reach because of the record amount of deposition that year. The inlet and outlet boundaries are S209 and S201, respectively (Figure 2). The observed flow and sediment discharges at Qingxichang (QXC) and Wanxian (WX) gauging stations have been depicted in Figure 3a,b, respectively. It clearly illustrates that the flow and sediment hydrographs are synchronous with each other at the two stations after the sediment discharges at WX station have been moved forward by one day. Considering the differences of hydrographs between the two stations and the contributions of tributary inflows are small, the interpolation method has been applied to calculate the incoming flow boundary condition at the HHC reach. The incoming suspended sediment concentration (SSC) boundary condition should be calculated through Equation (17). As the distance ratio of QXC-HHC to HHC-WX is equal to the ratio of the amount of deposition at QXC-HHC to that at HHC-WX in 2012, approximately 3:2 [52], and the flow discharges at the two stations are nearly the same, the interpolation method can be applied to approximate the SSC at the inlet boundary as well. The RMSE (Root Mean Square Error) value of the calculated SSC through the above two methods is 0.05 kg/m3, which is acceptable for sediment deposition is negligible when the SSC is less than 0.1 kg/m3. Besides, the SSC propagation is supposed to delay for one day from QXC station to the HHC reach. The water stage measured at Shibaozhai station is used as the outlet boundary condition (Figure 2). Only the flood season from May to November is simulated instead of a whole year because most sediment is transported during this period (Figure 4b), similar to the method applied by Fang and Rodi [53] to study the sedimentation of near dam region after TGP impoundment. This duration has been divided into six periods based on the water stage process (Figure 4a). It should be noted that the water stage rising during the last period of this process is resulted from the operation of TGP and the water stage and bed elevation data are both based on Wusong base level.

$$\mathcal{S}\_{\text{\\_}HHC} = \left[0.4(Q\_{\text{\\_}QX\mathbb{C}} - Q\_{\text{\\_}WX}) + Q\_{\text{\\_}WX}\right] / Q\_{\text{\\_}HHC} = \left(0.4Q\_{\text{\\_}QX\mathbb{C}} + 0.6Q\_{\text{\\_}WX}\right) / Q\_{\text{\\_}HHC} \approx 0.4S\_{\\_}\_{\text{\\_}WX} + 0.6S\_{\\_}\_{WX} \quad \text{(17)}$$

where *QS* = *Q* <sup>×</sup> *S*, kg/s; *Q* = flow discharge, m3/s; and *S* = sediment discharge, kg/s; 0.6 and 0.4 represent percentage of amount of sediment deposition at QXC-HHC and HHC-WX, respectively.

**Figure 3.** (**a**) Hydrograph at Qingxichang (QXC) and Wanxian (WX) gauging stations. (**b**) Sediment discharge (*QS*) measured at QXC and WX and calculated at HHC (The *QS* at WX station has been moved forward by one day).

**Figure 4.** (**a**) Hydrograph and water stage from May 1 to November 1 (*Q* and *ZS* represent discharge and water stage, respectively); the black filled circles divide the duration into several periods descripted clearly by the vertical black dash lines. (**b**) Suspended sediment concentration (SSC) as the inlet boundary in this duration.

A median size of 0.008 mm is used to represent the inflow cohesive sediment composition of this reach [13]. A flood event on 16 July 2012 is chosen as a verification case for this river reach simulation. Table 3 lists parameters and conditions of it. Because the radius to width ratio (*r*/*w*) is in the range of 0.8 to 2.0 (Table 3), this river reach belongs to sharply curved bends. The computation domain of the river reach is divided into 211 × 41 grids in longitudinal and transverse directions, with time steps of 1.0 s and 60.0 s for flow and sediment calculation, respectively.


**Table 3.** Channel dimensions and flow condition of HHC reach.

#### **4. Results**

The flow simulation results of L, B, and NL models are verified for the discharge of 30,200 m3·s<sup>−</sup>1, and the model with best performances has been selected. The preferable model L and the reference model N are used to predict cohesive sediment deposition during an annual hydrograph. The basic parameters, such as eddy viscosity coefficient and roughness of flow module, and parameters of sediment module are calibrated in N model first and then applied to the other models.

#### *4.1. Verifications*

#### 4.1.1. Flow

Figure 5 shows simulated water stage at the right bank and the depth-averaged velocities across the channel width of the HHC reach. It can be seen that the results of the L model are more reasonable than those of the other models. The velocity shift due to secondary flow can be well predicted by the L model at the end of the bends (S202 and S206), especially at the exit of the second bend (S202), in contrast to other models. In addition, as the high velocity core shifts to the right bank at the end of the first bend (S206), velocity of the left branches (S205) has been reduced. That explains why the velocities predictions by B and L models are lower than those by N and NL models at S205 (Figure 5c). Overall, the differences among B, N and NL models are small, while the L model is preferable according to the flow simulation results of the HHC reach.

**Figure 5.** (**a**) Water stage of the right bank (downstream view). (**b**–**d**) Depth-averaged velocity distribution measured and predicted by N, B, L, and NL models at three cross-sections for discharge 30,200 m3/s.

To quantitatively assess the performances of different models in flow simulation of the HHC reach, the RMSE of water stage and velocities of different models at typical cross-sections are listed in Table 4. The L model with the smallest RMSE results outperforms the other models at the discharge of 30,200 m3/s.


**Table 4.** The RMSE of water stage (rows 1–3) and velocities (rows 4–9) of different models.

#### 4.1.2. Sediment

Based on the above flow simulation results, the L model has been applied to the HHC reach to investigate the secondary flow effects on cohesive sediment deposition. The results of N model serve as references.

The deposition module is verified by field measurements (Figure 6a) in terms of planar distribution of deposition (Figure 6b,c), bed elevation (Figure 7), and amounts of deposition. Figure 6a–c show that the simulated planar distribution of deposits by the L and N models agree with field measurements qualitatively, with the maximum thickness of deposits found at the convex bank of the first bend, and the majority of deposits located at the right bank of the inlet and the left branch of the reach. The predicted thickness of sediment deposits by the L model is approximately 1 m thicker than that by N model on the concave bank of the first bends (region 1, Figure 6d), which is much closer to the measurement 5–7 m (Figure 6a). Bed elevations simulated by the two models matches well with measurements at S204–S206 (Figure 7). Predictions of total amounts of deposition from S206 to S203 are 8.33 <sup>×</sup>10<sup>6</sup> m3 and 8.0 <sup>×</sup> <sup>10</sup><sup>6</sup> m3 by the N and L models respectively, while the field measurement during the same period is 8.18 <sup>×</sup>10<sup>6</sup> m3 [13]. The relative error is around 2%, which qualify the sediment module used in this paper. In general, the L model performs better than the N model in predicting the planar distribution of cohesive sediment deposition.

**Figure 6.** (**a**) Planar distribution of sediment thickness measured, the maximum is 7 m from March to August, 2012. (**b**) Sediment thickness simulated by the L model (**c**) and N model. (**d**) The difference between the L and N models.

**Figure 7.** Comparison of bed elevation at cross-sections between measurements and predictions.

#### *4.2. Secondary Flow E*ff*ects on Cohesive Sediment Deposition*

The differences in planar distribution and amounts of deposition predicted by the L and N models have been illustrated in Figures 6d and 8, respectively, which clearly suggest the secondary flow effects on cohesive sediment deposition. Due to its impacts, high velocity core shifts from the convex to the concave bank of the bend, leading to the redistribution of bed shear stress and the consequent morphological changes [9]. Shifts of high velocities predicted by the L model result in the more deposition in region 1, 5, and 6 and less deposition in regions 3 and 4. The increase of sediment deposits in region 1 reduces sediment transported to region 2, resulting in less deposition here. The difference of predicted amount of deposition between the two models is about 0.31 <sup>×</sup> <sup>10</sup><sup>6</sup> m3 from 11 September to 1 November, as is clearly shown in Figure 8. This difference is small compared to the total amount of deposition during the whole year, approximately 8.0 <sup>×</sup> 10<sup>6</sup> m3. However, this difference can accumulate if the water stage keeps rising due to the impoundment of TGP. In general, secondary flow effects on cohesive sediment deposition become more obvious in the last period of the annual hydrograph when the sediment load is low and water stage is high (Figure 4).

**Figure 8.** Differences in deposition volume during different periods (average SSC means the average suspended sediment concentration during each period).

The total deposition volume is calculated from S203–S206 during different periods of this year, because this part of the reach is seldom affected by the inlet and outlet boundaries. Deposition of this part is greatly impacted by the velocity redistribution at S206 (e.g., Figure 5c,d), which is controlled by the secondary flow produced in the upstream bend and the bed topography (transverse bed slope) there. In addition, the sediment load plays an important role in the deposition of this part. Therefore, the average of suspended sediment load during different periods has been shown in Figure 8 as well. When the sediment load is low, the velocity redistribution plays a dominate role resulting in more sediment transport downstream and less deposition due to the shift of high velocities to the right branch. Otherwise, the situation is just reversed, and more deposits can occur in the left branch resulting from the huge amount of sediment transported, despite of the fact that the velocities are higher in the right branch. These can qualitatively explain the difference in predicted amounts of deposition during different periods except the fifth period (1–11 September). In that period, the transverse bed slope at S206 is high enough to strengthen velocity redistribution further, thus surpasses the effects of higher sediment load and result in less predicted deposition by the L model than the N model. During the last period (11 September–1 November), the significant difference of predicted deposition volume is resulted from both the low sediment load transport and the large transverse bed slope.

Figure 9 shows the predicted depth gradients (a) and velocity distributions (b) by the L and N models at S206 on 5 June and 18 September (as typical days of the first and last periods), respectively, illustrating the effects of bed topography. It clearly reveals that the velocity redistribution on 18 September is resulted from the bed topography effects as the sediment load on the two days is ~0.1–0.3 kg·m<sup>−</sup>3. In all, the low sediment load and the velocity redistribution induced by secondary flow produced by upstream bend and the bed topography result in the difference deposition predictions by the two models.

**Figure 9.** (**a**) Depth gradient (represents bed topography effects). (**b**) Velocity distribution predicted by the N and L models at S206 on typical days of the first and last period, respectively.

#### **5. Discussion**

One of the most important physical processes in meandering rivers is the outward shifting of main flow velocity caused by secondary flow, which is driven by channel curvature or point bars bed topography [3]. The latter one is called topography steering [9], which plays a significant role in meander dynamics [3]. Whether and how the correction terms representing the secondary flow effects quantify this process and the performances of these models in meandering channels of different scales will be discussed in this part. Besides, secondary flow effects on the total amount of deposition of the aforementioned part of this reach (S203–S206) are controlled by the properties of cohesive sediment, which will be investigated as well.

#### *5.1. Secondary Flow E*ff*ects on Flow Field*

#### 5.1.1. Topography Effects

Equation (6) clearly reveals that the correction terms of the three models are directly proportional to the gradients of water depth (*H*). Due to the effects of bed topography, the longitudinal and transverse gradient of water depth in HHC reach is in the range of 0.01 to 0.001 and 0.01 to 0.1, respectively. Therefore the magnitudes of correction terms follow the same tendency as that of the gradients of water depths, in other words, the correction terms are able to reflect the topography effects. This finding has been justified by Lane [54] who pointed out that correction terms represent the gradients of the transport of momentum. Figure 10 depicts the distributions of (a) *S*<sup>ξ</sup> and (b) *S*<sup>η</sup> of the L, B and NL models along the channel. The orders of magnitude of them are within 0.01 to 0.001 in the longitudinal

direction, which is the same as the longitudinal gradient of water depth. In the transverse direction, the order of magnitude of the L model is 0.01–0.1, which is consistent with the transverse gradient of water depth, while those of the B and NL models are approaching to zero and in the range of 0.01 to 0.001, respectively. The smaller orders of magnitude of the two models are resulted from the methods of them. As to the B model [35], it only considers the longitudinal correction. As to the NL model [1], the sharpness of the HHC reach limits the growth of the secondary flow. Since the L model considers the corrections in both directions and has larger correction values than the other two models, it outperforms the other models in the flow simulation as shown in Figure 5. In addition, the simulation results shown in Figure 5 clearly indicate that 2D depth-averaged model that include secondary flow effects (e.g., the L and B models) should be given first priority when it comes to sharp meandering channels with bed topography, such as the HHC reach. This has been confirmed by de Vriend [55] who found that his mathematical model with considering secondary flow effects worked better for curved bend flow simulation over developed bed.

**Figure 10.** Correction terms (**a**) *S*ξ and (**b**) *S*η distributions of the L, B, and NL models along the channel.

#### 5.1.2. Applicability of Different Secondary Flow Models

The differences among these models are listed in Table 1, which mainly lie in whether considering the effects of phase lag (B and NL models), sidewall boundary conditions (B model), and bend sharpness (NL model). As the HHC reach is sharply curved bends, the saturation effect considered by the NL model has weakened the secondary flow effects, which result in the minor differences of simulation results between the NL and N models (Figure 5). The depth to width ratio (*H*/*w*) distinguishes between meandering channels of different scales. It is approximately 0.001–0.06 in the HHC reach at the discharge of 30,200 m3/s, while that in the laboratory bend channels and small meandering rivers are in the range of 0.05 to 0.25 [45] and 0.06 to 0.1 [56], respectively. Therefore, the effects of wall boundary conditions and phase lag have been reduced for such small value of *H*/*w*. Although B model has taken the bed topography effects into account in a similar way as the L model does, its correction terms only focus on the longitudinal direction. Consequently, the flow simulation results of the L model are better than that of the B model in the HHC reach. Overall, L model is preferable to flow simulation in meandering channels of mega scale, such as HHC reach. However, for laboratory scale curved bends with flat bathymetry, the B model obtains better results [45]. And for sharply curved bends of laboratory and small meandering rivers scales, the advantages of the NL model have been exhibited according to the flow simulation results by Blanckaert [1,2] and Ottevanger [57]. The *H*/*w* may play an important role, while the main reasons remain to be further investigated.

#### *5.2. Secondary Flow E*ff*ects on Deposition Amounts*

According to the deposition simulation results, secondary flow effects on the total deposition volume are small during an annual hydrograph (Figure 8). However, these effects vary with the changes of the cohesive sediment properties, such as settling velocity and critical shear stresses of cohesive sediment, which depend on the flow conditions and the process of bed consolidation. Series of numerical experiments are designed to investigate secondary flow effects on the deposition volume of cohesive sediment with different properties; these effects are reflected by the relative difference in deposition amounts (RD) predicted by the N and L models. Numerical experiments are conducted under the same flow condition (Table 3) to keep the strength of secondary flow constant in the HHC reach. The calculation time for each experiment is 33 days. Different properties of cohesive sediment (Table 5) are represented by the variation of settling velocity (ω*s*) and the critical shear stress for deposition (τ*cd*). Other parameters used in sediment module are the same as that of HHC reach.

**Table 5.** Settling velocity (ω*s*) and critical shear stress for deposition (τ*cd*) in numerical experiments and results.


<sup>1</sup> The relative difference in deposition amounts (RD) predicted by N and L models.

Calculated RD values are listed in Table 5. It is obtained by calculating the difference of the predicted amounts of deposition by L and N models, and then divided by the N model predictions. The negative value of it means the amount of deposition simulated by the L model is smaller than that by the N model. The relationships of RD against ω*<sup>s</sup>* and τ*cd* are shown in Figure 11. RD is in reverse linear proportion to ω*s*, which means the secondary flow effects on the deposition volume increase with the decrease of settling velocity of cohesive sediment. For τ*cd* is ~0.44 Pa, RD is approaching zero. It implies that secondary flow nearly has no effect on the total deposition volume while its effects on planar distribution can still exit (Figure 6d). As the τ*cd* increases, the secondary flow impacts on deposition become greater. In general, RD varies with the settling velocity and critical shear stress for deposition of cohesive sediment and the magnitudes of RD are within 11% based on the parameter values used here.

**Figure 11.** The relationships of relative difference in deposition volume (RD) predicted by the L and N models against (**a**) settling velocity (ω*s*) and (**b**) critical shear stress (τ*cd*).

#### *5.3. Future Reseach Directions*

1. As the study case is a reach of Yangtze River, which is classed as a mega river, secondary flow effects on bed morphology of meandering channels of different scales (natural rivers with different width to depth ratio) should be investigated. Besides, as the bank of HHC reach is nonerosional, the evolutions of natural rivers with floodplain consisting of cohesive sediment should be simulated by the 2D model developed here. In addition, long-term simulations, such as decadal timescales, should be considered in the future to research the cumulative effects of secondary flow.

2. As to the cohesive sediment transport, the values of parameters play important parts in the distributions and amounts of sediment deposition (Figure 11). The roles they played should be compared with that of secondary flow in bed morphology variations. More importantly, the erosion processes should be studied as these processes cannot be reflected obviously in the HHC reach.
