A Review of Numerical Simulations of Secondary Flows in River Bends

: River bends are one of the common elements in most natural rivers, and secondary ﬂow is one of the most important ﬂow features in the bends. The secondary ﬂow is perpendicular to the main ﬂow and has a helical path moving towards the outer bank at the upper part of the river cross-section, and towards the inner bank at the lower part of the river cross-section. The secondary ﬂow causes a redistribution in the main ﬂow. Accordingly, this redistribution and sediment transport by the secondary ﬂow may lead to the formation of a typical pattern of river bend proﬁle. It is important to study and understand the ﬂow pattern in order to predict the proﬁle and the position of the bend in the river. However, there are a lack of comprehensive reviews on the advances in numerical modeling of bend secondary ﬂow in the literature. Therefore, this study comprehensively reviews the fundamentals of secondary ﬂow, the governing equations and boundary conditions for numerical simulations, and previous numerical studies on river bend ﬂows. Most importantly, it reviews various numerical simulation strategies and performance of various turbulence models in simulating the ﬂow in river bends and concludes that the main problem is ﬁnding the appropriate model for each case of turbulent ﬂow. The present review summarizes the recent advances in numerical modeling of secondary ﬂow and points out the key challenges, which can provide useful information for future studies.


Introduction
Rivers rarely run on straight paths, and most of them run on a curved route. The flow patterns on river bends are very complicated and have special characteristics. Basically, flow patterns in river bends are affected by certain factors such as the centrifugal force caused by the flow curvature and the vertically varying velocity, the cross-sectional stress, in addition to the pressure gradient in the radial direction due to the lateral slope of the water's surface. All these factors create a flow called secondary flow (Ghobadian and Mohammadi, 2011) [1], which may significantly affect the primary flow motion and, in turn, affect pollutant dispersion and sediment transport. Therefore, it is important to predict and estimate the flow behavior of river bends for proper design of nearby hydraulic structures, environmental and ecological assessments, and sound planning of inhabitation projects.
Flow patterns and characteristics in river bends have been widely studied due to the major significance of this subject. The secondary flow, according to Rozovskii (1961) [2] and De Vriend (1979Vriend ( ,1980 [3,4], is a helical path moving toward the outer bank at the upper part of the river and moving toward the inner bank at the lower part of the river. Thus, a slight deviation occurs on the shear force, which follows the same direction as the local flow near the bed. This deviation is from the direction of the mean flow (Engelund and Skovgaard, 1973) [5].
The flow in river bends can be studied using either physical or numerical methods. Several experimental and laboratory studies were performed to investigate the flow fea- Figure 1. The mechanism of secondary flow (Jansen, 1994) [27].
The main flow of the curved channel is affected by the redistribution of the velocity. The velocity increases near the outer bank of the channel and decreases near the inner bank of the channel (De Vriend, 1980) [4]. As the river enters the bend, the flow particles are affected by the acceleration due to centrifugal force. Accordingly, the flow particles tend to move towards the external wall of the river. Then, an increase in the water level occurs near the outside wall and decreases near the inside wall. Non-uniform pressure distribution is generated in the cross-section due to the slope of the free surface. This pressure increases near the external wall. Due to the difference between the small vertical accelerations and the high gravitational acceleration, there is a distribution of the hydrostatic pressure over the vertical. Subsequently, there is approximate equality in the pressure gradient in the radial direction over the vertical. Thus, the magnitude of the longitudinal velocities is equal at all points in the vertical direction and the centrifugal accelerations are completely counterbalanced by the centripetal acceleration. Due to the viscosity forces, there is a non-uniform distribution of the longitudinal velocity, which, in turn, causes a non-uniform distribution of the centrifugal acceleration. Accordingly, the velocity and the centrifugal acceleration increase near the flow surface and decrease near the bottom. The helical movement that is vertical to the main flow occurs because of the imbalance between the separate points in a vertical direction as shown in Figure 2. The movement of the helical path depends on the centrifugal acceleration, which means it is directed towards the outer bank at the upper part of the river as the centrifugal acceleration increases and towards the inner bank at the lower part of the river as the centrifugal acceleration decreases at this area. As the flow enters the bend, the free surface slope must be built up. Therefore, there is a negative pressure gradient near the internal wall and another positive pressure gradient close to the external wall (Jongbloed, 1996) [28].  (Jansen, 1994) [27].
The main flow of the curved channel is affected by the redistribution of the velocity. The velocity increases near the outer bank of the channel and decreases near the inner bank of the channel (De Vriend, 1980) [4]. As the river enters the bend, the flow particles are affected by the acceleration due to centrifugal force. Accordingly, the flow particles tend to move towards the external wall of the river. Then, an increase in the water level occurs near the outside wall and decreases near the inside wall. Non-uniform pressure distribution is generated in the cross-section due to the slope of the free surface. This pressure increases near the external wall. Due to the difference between the small vertical accelerations and the high gravitational acceleration, there is a distribution of the hydrostatic pressure over the vertical. Subsequently, there is approximate equality in the pressure gradient in the radial direction over the vertical. Thus, the magnitude of the longitudinal velocities is equal at all points in the vertical direction and the centrifugal accelerations are completely counterbalanced by the centripetal acceleration. Due to the viscosity forces, there is a non-uniform distribution of the longitudinal velocity, which, in turn, causes a non-uniform distribution of the centrifugal acceleration. Accordingly, the velocity and the centrifugal acceleration increase near the flow surface and decrease near the bottom. The helical movement that is vertical to the main flow occurs because of the imbalance between the separate points in a vertical direction as shown in Figure 2. The movement of the helical path depends on the centrifugal acceleration, which means it is directed towards the outer bank at the upper part of the river as the centrifugal acceleration increases and towards the inner bank at the lower part of the river as the centrifugal acceleration decreases at this area. As the flow enters the bend, the free surface slope must be built up. Therefore, there is a negative pressure gradient near the internal wall and another positive pressure gradient close to the external wall (Jongbloed, 1996) [28].  (Jongbloed, 1996) [28].
At the exit of the bend, there is a positive pressure gradient near the interior wall and a negative one close to the exterior wall. The free surface slope decreases towards the bend exit until it may gradually become zero at the end of the bend.

Secondary Flow Cells
There are two regions that could be found in the bend of the curved channel: the center region cell and the outer bank cell. The center region cell is the classical helical or spiral motion of the secondary flow that moves toward the outer bank at the upper part of the river and moves toward the inner bank at the lower part of the river. The center region cell results from the interaction of the centrifugal force, the pressure gradient resulting from the inclined water level, and the roughness of the bed (Rozovskii, 1961) [2]. The outer bank cell is a counterrotating eddy at the outer bank of the channel (Bathurst et al., 1977) [29]. The outer bank cell is mainly caused by turbulence anisotropy and is a key feature of bend flow, and it has been found that the proximity of the water surface and the outer bank can affect its generation and characteristics (Blanckaert, 2009;Blanckaert, 2011) [30,31]. The outer bank cell is usually smaller and weaker than the center region cell. However, it should not be underestimated because it could be strongest in the sharp channel bends, and it has a stabilizing effect on the erosion of the outer bank and the mixing abilities in this region (Graf and Blanckaert, 2002) [32]. This is shown in Figure 3.  (Jongbloed, 1996) [28].
At the exit of the bend, there is a positive pressure gradient near the interior wall and a negative one close to the exterior wall. The free surface slope decreases towards the bend exit until it may gradually become zero at the end of the bend.

Secondary Flow Cells
There are two regions that could be found in the bend of the curved channel: the center region cell and the outer bank cell. The center region cell is the classical helical or spiral motion of the secondary flow that moves toward the outer bank at the upper part of the river and moves toward the inner bank at the lower part of the river. The center region cell results from the interaction of the centrifugal force, the pressure gradient resulting from the inclined water level, and the roughness of the bed (Rozovskii, 1961) [2]. The outer bank cell is a counterrotating eddy at the outer bank of the channel (Bathurst et al., 1977) [29]. The outer bank cell is mainly caused by turbulence anisotropy and is a key feature of bend flow, and it has been found that the proximity of the water surface and the outer bank can affect its generation and characteristics (Blanckaert, 2009;Blanckaert, 2011) [30,31]. The outer bank cell is usually smaller and weaker than the center region cell. However, it should not be underestimated because it could be strongest in the sharp channel bends, and it has a stabilizing effect on the erosion of the outer bank and the mixing abilities in this region (Graf and Blanckaert, 2002) [32]. This is shown in Figure 3.

Velocity Redistribution
The flow distribution is nearly symmetric in straight channels, but the main velocity of the flow in channel bends is significantly affected by the secondary flow, and accordingly, a considerable redistribution occurs with respect to the main velocity. The flow par-

Velocity Redistribution
The flow distribution is nearly symmetric in straight channels, but the main velocity of the flow in channel bends is significantly affected by the secondary flow, and accordingly, a considerable redistribution occurs with respect to the main velocity. The flow particles that are close to the free surface are transported towards the outer bank of the river affected by the secondary flow. The core of the high-velocity region is near the inner bend at the bend entrance due to the potential vortex effect (Zeng et al., 2008) [24] and gradually shifts towards the outer bank due to secondary flow caused by the combined effects of pressure gradients and centrifugal forces. The velocity and momentum of these particles are larger compared to the ones in the outer bank. Thus, an increase occurs in the velocity of the outer bank region. The process is just the opposite near the bed of the river, where the flow particles with smaller velocity and momentum move towards the inner bank of the river, causing a low velocity at this area. The vertical convective transport of main flow momentum is significant at the area close to the sidewalls, while it is almost negligible near the central part of the river. The main reason for the velocity increase near the outer bank and reduction near the inner bank is vertical convection, while outward lateral interaction is due to radial convection (Jongbloed, 1996) [28]. The effect of secondary flow on the main flow is explained in Figure 4, where the streamlines represent the secondary flow and the isovels represent the main flow velocity.

Velocity Redistribution
The flow distribution is nearly symmetric in straight channels, but the main velocity of the flow in channel bends is significantly affected by the secondary flow, and accordingly, a considerable redistribution occurs with respect to the main velocity. The flow particles that are close to the free surface are transported towards the outer bank of the river affected by the secondary flow. The core of the high-velocity region is near the inner bend at the bend entrance due to the potential vortex effect (Zeng et al., 2008) [24] and gradually shifts towards the outer bank due to secondary flow caused by the combined effects of pressure gradients and centrifugal forces. The velocity and momentum of these particles are larger compared to the ones in the outer bank. Thus, an increase occurs in the velocity of the outer bank region. The process is just the opposite near the bed of the river, where the flow particles with smaller velocity and momentum move towards the inner bank of the river, causing a low velocity at this area. The vertical convective transport of main flow momentum is significant at the area close to the sidewalls, while it is almost negligible near the central part of the river. The main reason for the velocity increase near the outer bank and reduction near the inner bank is vertical convection, while outward lateral interaction is due to radial convection (Jongbloed, 1996) [28]. The effect of secondary flow on the main flow is explained in Figure 4, where the streamlines represent the secondary flow and the isovels represent the main flow velocity.   Figure 4 presented above shows that the variance of velocities along the streamlines of the secondary flow near the bed of the river is not significant. The difference increases towards the inner bank, and then streamlines are almost vertical to the isovels near the surface of the river. There are some factors that affect the secondary flow, such as the distance between the isovels (the local main velocity gradient), the secondary flow strength, and the sine of the angle of intersection between the main flow isovels and the streamlines of the secondary flow. Note that the effect of the secondary flow in the wall district will be comparatively stronger in the upper part of the section where the velocity gradient is of the same order of magnitude throughout the vertical. As a result, the main velocity profile becomes flatter near the surface, which is confirmed by experimental research (Jongbloed, 1996) [28].
On the other hand, the longitudinal velocity component at the entry of a bend and in the early part of the internal wall area increases, accompanied by a steep fall of the free surface. In addition, extreme velocity is found at the free surface at this place in the crosssection (Booij and Tukker, 1996) [33]. It is noticed that the maximum velocity at the end of a 180° curved bend shifts outward-the outer part of the cross-section-due to the redistribution of the main velocity (De Vriend, 1980) [4], while the opposite situation is found at the outside wall area at the entry of a bend. Note that the longitudinal velocities decrease due to a positive pressure gradient.   Figure 4 presented above shows that the variance of velocities along the streamlines of the secondary flow near the bed of the river is not significant. The difference increases towards the inner bank, and then streamlines are almost vertical to the isovels near the surface of the river. There are some factors that affect the secondary flow, such as the distance between the isovels (the local main velocity gradient), the secondary flow strength, and the sine of the angle of intersection between the main flow isovels and the streamlines of the secondary flow. Note that the effect of the secondary flow in the wall district will be comparatively stronger in the upper part of the section where the velocity gradient is of the same order of magnitude throughout the vertical. As a result, the main velocity profile becomes flatter near the surface, which is confirmed by experimental research (Jongbloed, 1996) [28].
On the other hand, the longitudinal velocity component at the entry of a bend and in the early part of the internal wall area increases, accompanied by a steep fall of the free surface. In addition, extreme velocity is found at the free surface at this place in the crosssection (Booij and Tukker, 1996) [33]. It is noticed that the maximum velocity at the end of a 180° curved bend shifts outward-the outer part of the cross-section-due to the redistribution of the main velocity (De Vriend, 1980) [4], while the opposite situation is found at the outside wall area at the entry of a bend. Note that the longitudinal velocities decrease due to a positive pressure gradient.
) in the inner wall region of a shallow channel (d/B = 0.1) (De Vriend, 1980) [4]. Figure 4 presented above shows that the variance of velocities along the streamlines of the secondary flow near the bed of the river is not significant. The difference increases towards the inner bank, and then streamlines are almost vertical to the isovels near the surface of the river. There are some factors that affect the secondary flow, such as the distance between the isovels (the local main velocity gradient), the secondary flow strength, and the sine of the angle of intersection between the main flow isovels and the streamlines of the secondary flow. Note that the effect of the secondary flow in the wall district will be comparatively stronger in the upper part of the section where the velocity gradient is of the same order of magnitude throughout the vertical. As a result, the main velocity profile becomes flatter near the surface, which is confirmed by experimental research (Jongbloed, 1996) [28].
On the other hand, the longitudinal velocity component at the entry of a bend and in the early part of the internal wall area increases, accompanied by a steep fall of the free surface. In addition, extreme velocity is found at the free surface at this place in the cross-section (Booij and Tukker, 1996) [33]. It is noticed that the maximum velocity at the end of a 180 • curved bend shifts outward-the outer part of the cross-section-due to the Water 2021, 13, 884 6 of 32 redistribution of the main velocity (De Vriend, 1980) [4], while the opposite situation is found at the outside wall area at the entry of a bend. Note that the longitudinal velocities decrease due to a positive pressure gradient.

Mathematical Description of the Main Flow and Secondary Circulation
The differential equations, which include the continuity equation and Reynolds-Averaged Navier-Stokes system, can be used to describe the three-dimensional incompressible fluid flow. For uniformity, these equations can be rewritten as below: Conservation of mass: The momentum equations: where U is the velocity vector; ρ is the density; g is the gravity acceleration; p is the pressure; T is the viscous stress tensor.

Boundary Conditions
The previous section introduced the mechanisms of secondary flow, the secondary flow cells, the velocity redistribution, and the governing equations. To simulate bend flows, certain boundary conditions should be applied, and the present section introduces the typical boundary condition setup.
There are four main kinds of boundaries for the numerical models that should be considered, namely: walls, free surface, inlet, and outlet.

Wall Boundaries
As the velocities at the wall are equal to zero and have large gradients in the vicinity, a very small grid distance should be used at this position. Thus, the high-Reynolds number model is no longer viable, and therefore, the first grid-point is defined at some distance from the wall (Jongbloed, 1996) [28]. Considering the effect of the viscous shear stresses or the turbulent shear stresses or both, it is recommended to separate the near-wall area into three regions (Booij, 1992) [34]: wall region, buffer layer region, and the turbulent boundary layer. At the wall region, the turbulent shear stresses can be abandoned (the so-called viscous sublayer). A buffer layer region comes next to the wall region where neither of the viscous shear stresses nor the turbulent shear stresses can be abandoned. The third region, the turbulent boundary layer, is distinguished where the viscous shear stresses can be abandoned. Accordingly, the boundary conditions and the first grid point should preferably be positioned in the last zones (Pezzinga, 1990) [35]. The boundary conditions for k (turbulent kinetic energy) and ε (turbulent energy dissipation rate) can be found using the following equations (Jongbloed, 1996) [27]: where u * is shear velocity, c µ is a model constant ( ∼ = 0.09), and y is the distance to the wall. Flow variation properties at the area close to the solid walls and bed are steep. To simulate the wall effect and estimate the turbulence parameters at this area, the wall function method is employed (Rodi, 1984) [36]. With the wall function, a fairly coarse grid is used in the area close to the wall. It is also considered to be an economical method of turbulent flow modeling (Jian and McCorquadale, 1998) [37]. Below are some of the turbulence equations for some models that use the wall function approach (Wilcox, 1998) [38]: where δ np = distance from the solid boundary to the center of the closest element; k p = turbulence kinetic energy at the closest cell to the wall; and u + p = defined as follows: where y + p = δ np u * /υ and y + 0 = 11.53 (Melaaen, 1990) [39] and E is related to the roughness Reynold number k + s = u * k s /υ as follows (Wu et al., 2000) [40]: where k s = equivalent roughness height of bed; B = additive constant = 5.2; ∆B = roughness function of the roughness Reynold number, which is calculated by (Cebeci and Bradshaw, 1977) [41]: The equivalent roughness height k s quantified the influence of roughness elements.

Free Surface
For bend flow, the rigid-lid approximation is often adopted. The rigid-lid approximation is used when the Froude number is small, and a frictionless flat plane is assumed (Leschziner and Rodi, 1979) [26]. This approximation for this fictitious rigid lid does not neglect the changes in surface elevation but these changes are taken into account in an indirect manner. Accordingly, non-zero pressure gradients in the radial and vertical directions are expected at the boundary, which would result in some errors because the superelevated regions are implicitly considered hydrodynamically passive. However, experiments conducted by Mcguirk and Rodi (1977) [42] showed that the error associated with this approximation does not exceed approximately 10% of the total channel depth (Mcguirk and Rodi, 1977) [42]. Therefore, the following assumptions at the free surface can be formulated (Jongbloed, 1996) [28]: The derivative in the vertical direction of the radial and longitudinal velocities are zero.
On the water surface, in absence of wind, the boundary conditions for the turbulence parameters are: For the free surface, the turbulent kinetic energy and net fluxes of horizontal momentum are set as zero. The pressure normal to the surface is set to atmospheric, velocity is set to zero, and the dissipation rate is calculated from the relation given by (Rodi, 1993) [43]: where h = local water depth.

Inlet
Determining the distributions of velocities and the parameters of the turbulence across the boundary of the inlet is likely not possible in practice. Therefore, most of these distributions are estimated based on physical reasoning. The experimental distribution of k is rarely available and the ε-distribution can only be obtained indirectly from measurements of other turbulent quantities (Leschziner and Rodi, 1979) [26]. To smooth out the boundary conditions from the positions of interest, the boundaries must lie at a distance long enough for the flow if the errors caused by the assumptions for the boundary conditions cannot be neglected (Jongbloed, 1996) [28]. Below are some of the turbulence equations for some models that are used at the inlet: where K in = inlet turbulent kinetic energy; ε in = inlet turbulent kinetic dissipation rate; U in = average inlet velocity; T i = intensity of turbulence; and L and = characteristic length of channel and length scale of the flow field, respectively.

Outlet
The discussion in the previous section can be used to determine the outlet boundary conditions. If the outlet boundary is far enough away, a zero gradient boundary condition may be employed for most variables.

Advantages and Disadvantages of Various Numerical Methods
Examples of the widely used turbulence modeling approaches include direct numerical simulation (DNS), Reynolds-averaged numerical simulation (RANS), RSM (Reynolds Stress Model), large eddy simulation (LES), and detached eddy simulation (DES). DNS is the most reliable approach, but its computational costs are super high and are, thus, not further discussed in this paper. The other turbulence modeling approaches are discussed in the present section.

RANS Modeling
For RANS models, the standard k−ε turbulence model (Versteeg and Malalasekera, 2007) [44] is considered to be relatively simpler. A very good performance can be achieved from the k-ε model in many cases as the model is well established. The model has also been extensively validated. However, it is considered more expensive than some other turbulence models, such as the mixing length model. The performance of the k-ε model may not be satisfactory in some cases, such as unconfined flows, flows with large extra strains, rotating flows, and flows driven by the anisotropy of normal Reynolds stresses.
The realizable k-ε model has been found to be more precise than some other models (e.g., the standard k-ε model) in predicting flows such as discrete flows and flows with complicated secondary flow features (Gildeh, 2013) [45].
The blending function drives the standard k−ω model in close wall regions and drives the k-ε family models in areas further from the surface. These differences make the SST model more accurate for a greater variety of flows than the standard k-ω model (Menter, 1994) [46].

RSM Modeling
RSM models may be considered as a subcategory of RANS models while they are computationally more expensive. According to Davidson (2011) [47], the advantages of RSMs could be summarized as follows: (1) the ability of RSM to respond effectively to any sudden change in the main strain; (2) modelling is not necessarily required for the production terms, which are used to explain many phenomena; and (3) the transport equations, which contain the main physical mechanisms that derive the evolution of turbulence such as production, redistribution, turbulent transport, viscous diffusion, and dissipation, are solved without assuming any behavior for Reynolds stresses. However, implementing the RSMs is considered to be a challenging task due to their complexity and CPU time consumption. Additionally, the numerical analyses of RSMs are very sensitive to influencing factors because of the small stabilizing of the second order derivatives of the momentum equations.

LES Modeling
Compared with the two-equation turbulence models (e.g., k-ε, k-ω), LES is typically more expensive in terms of computational costs. However, when comparing the computational cost with the accuracy of LES models, its efficiency, and other features with regard to computational time, can also be reasonable and worthy. This type of turbulence model can be more universally used since it overcomes the issues associated with using small eddies, such as tending to be more isotropic and less dependent on geometry. However, it is more computationally expensive since it needs finer grid densities for LES models than most RANS models. Although the large CPU costs for high-Re flows can be reduced by using a coarse near-wall mesh coupled with wall functions, the discretization schemes should be carefully designed to obtain more accurate results (Cable, 2009) [48].

DES Modeling
The DES model is a cross between the LES and RANS approaches. It is quite beneficial to combine the RANS modeling approach with the LES approach if the usage of LES models is computationally unaffordable. Typically, the DES method is much less expensive than the LES approach but more expensive than the RANS. In regions where large turbulence scales play a dominant role, the LES model based on the sub-grid model by the DES model is used. In areas where viscous impacts predominate, such as near-wall areas, the RANS model is used (Shur et al., 1999) [49].

Research on River Bends Flows and Key Findings
The previous sections have introduced the fundamentals of secondary flow, governing equations, and boundary conditions. Previous researchers have conducted numerical studies on river bend flows based on these fundamentals and configurations. This section reviews and discusses some of these studies and their key findings. The reviewed papers were mainly selected based on their research focuses, and they are classified into four categories: secondary flow in the curved channel; secondary flow and pollutant dispersion; secondary flow and sediment transport; and secondary flow in bend over topography.

Secondary Flow in a Curved Channel
In this section, some previous numerical research studies on secondary flow in channel bends were discussed, taking into account different bend geometries, and different models of turbulence. Then, each numerical model was examined and evaluated. Table 1 summarizes these numerical research studies and their key findings. Table 1. Summary of selected studies about secondary flow in curved channels used in Section 5.1.

Shape of Channel Numerical Model Key Findings Reference
S-shape with two identical 90 • reverse bends Standard k-ε low-Reynolds number k-ω k-ω model is better than the k-ε model in simulating the transverse velocities, which form the secondary flow phenomenon. (Khosronejad et al., 2007) [20] U-shaped mildly curved channel having a U-shape Depth-averaged 2D numerical model Generally, there was an agreement between the simulated and measured surface elevations at the inner, outer bank, and centerline. (Duan, 2004) [50] U-shaped sharply curved channel in 180 • bend

Depth-averaged 2D numerical model
There is good agreement between the simulated and measured results. (Duan, 2004) [50] A channel in 90 • bend low-Reynolds number k-ω, k-ε (realizable), k-ε (RNG) and LES models The performance of k-ε (RNG), k-ω and k-ε (realizable) models was good in simulating the flow in the bend, while the LES model has not shown its capability in simulating the flow in channel bends, especially after the curve.
The results of the comparison between the numerical model and the experimental data indicated that there is an acceptable level of agreement. The flume of two successive 180 • bends Standard k-ε model, and the non-linear k-ε model The two models showed the counter-clockwise secondary cell inherited from the first bend, while only the non-linear k-ε model showed the counter-rotating secondary cell of the previous bend, which is coincident with the higher primary velocities. A channel in 135 • bend with side slopes Standard k-ε, (RNG) k-ε, realizable k-ε, k-ω, and (SST) k-ω) The realizable k-ε model was the best in capturing the outer bank cells among the tested models. The standard k-ε model and the low-Reynolds number k-ω model were used by Khosronejad et al. (2007) [20] to simulate the flow in a channel bend. To validate the hydrodynamic model results, the experimental data of Ghanmi (1999) [57] were used. In this experiment, an S-shaped flume with two identical 90 • reverse bends connected by a straight channel segment was used. The geometric parameters of the flume include the length of the straight channel segment (3.5 m), the rectangular channel cross-section (0.5 m wide and 0.4 m depth), and the inner radius of each channel bend (R 0 = 1.25 m). The characteristics and flow conditions considered in the experiment were as follows: the channel boundary was made with aluminum covered with a sand layer with a median diameter, d 50 , of 2 mm, the discharge was 0.023 m 3 /s, the flow depth was h = 0.15 m, the mean inflow velocity was U = 0.3 m/s, and 2D velocities (longitudinal and transverse) were measured with a laser Doppler anemometer at 20 flume sections. The boundary conditions applied in this analysis included: inlet, outlet, free surface, and wall. At each section, three verticals were measured with five points for each vertical as shown in Figure 5. The results showed that the supreme velocity at the bend happens near the inner bank and moves toward the opposite bank as the outlet of the bend is approached. The streamwise velocity component differs from zero at or close to the bed to a maximum value when approaching the water surface. This situation is expected because it represents frictional effects. Based on the previous point, it can be intuitively concluded that the centrifugal force is greatest near the water surface and decreases with depth. The flow vectors near the bed are directed towards the inner bank and at the free surface towards the outer bank. The analysis showed that the mean absolute errors of the k-ω model (10.4%) are lower than the k-ε model (12.6%) in simulating the transverse velocities. Although the difference is relatively low, it could be concluded that the k-ω model is better than the k-ε model in simulating the transverse velocities which form the secondary flow phenomenon (Falcon 1984) [58]. The experiment of De Vriend (1979) [3] was used by Duan (2004) [50] to simulate the flow in a mildly curved channel. The channel comprised a 1.7-m-wide flume having a Ushaped ground plan, with vertical sidewalls and horizontal bottom. The radius of curvature of the flume axis in the bend was 4.25 m, and the length of both straight upstream and downstream was about 6.0 m. The ratio of radius of curvature to channel width was 3.5. The discharge at the inlet was 0.0671 m 3 /s, the averaged velocity was 0.202 m/s, and the averaged flow depth was 0.1953 m. The boundary condition of the upstream was the inflow unit discharge, and for the downstream was the measured flow depth. By including the dispersion terms, a depth-averaged 2D numerical model was used to simulate the The experiment of De Vriend (1979) [3] was used by Duan (2004) [50] to simulate the flow in a mildly curved channel. The channel comprised a 1.7-m-wide flume having a U-shaped ground plan, with vertical sidewalls and horizontal bottom. The radius of curvature of the flume axis in the bend was 4.25 m, and the length of both straight upstream and downstream was about 6.0 m. The ratio of radius of curvature to channel width was 3.5.
The discharge at the inlet was 0.0671 m 3 /s, the averaged velocity was 0.202 m/s, and the averaged flow depth was 0.1953 m. The boundary condition of the upstream was the inflow unit discharge, and for the downstream was the measured flow depth. By including the dispersion terms, a depth-averaged 2D numerical model was used to simulate the flow in this curved flume, taking into account that various sets of longitudinal and secondary flow profiles for governing equations in Cartesian coordinates were used to derive the dispersion terms. The simulated velocity vector field using a 2D numerical model and surface elevation in shaded color are shown in Figure 6. As can be observed in this figure, the higher velocity is directed towards the inner bank of the bend, while the lower velocity is directed towards the outer bank of the bend. The main reason behind the increase in velocity at the inner bank is the transverse convection of momentum induced by the secondary flow. However, the effect of the secondary flow in the mild bend is very weak because the radius of channel curvature is much larger than the width of the channel. Furthermore, in the bend, the surface water level increases at the outer bank and decreases at the inner bank. The centrifugal force is the main reason behind the increased flow at the outer bank. Generally, there was an agreement between the simulated and measured surface elevations at the inner, outer bank, and the centerline (Duan, 2004) [50]. The experiment of Rozovskii (1961) [2] was used by Duan (2004) [50] to apply the depth-averaged 2D numerical model by including the dispersion terms 2D numerical model in a sharply curved open channel. It consisted of a 180° curved reach with a 6-mlong straight approach and a 3-m-long straight exit. The ratio of mean radius of curvature to width was 1.0. The width of the channel was 1.7 m. The cross-section of the channel was rectangular. Water depth at the downstream end was 0.053 m, and the discharge was 0.0123 m 3 /s. The surface elevation and velocity distribution are shown in Figure 7 below. There is an acceleration in the longitudinal flow close to the outer bank of the upstream part of the bend and a deceleration in the flow velocity near the outer bank. The surface elevation increased at the outer bank and decreased at the inner bank. There was good agreement between the simulated and measured results. The impact of secondary flow was observed on the depth-averaged flow distribution and water surface elevations, which became more visible with the increasing of the channel curvature. Due to the transverse bed slope in natural rivers, the velocity of secondary flow near the bed is higher. There is a significant role of the secondary flow in transporting bedload and suspended load in the transverse direction (Duan et al., 2001) [59]. Obviously, a 3D model is more accurate than that of a depth-averaged model. However, the depth-averaged two-dimensional model with the dispersion terms according to this study was capable of simulating the hydrodynamic flow in meandering channels. As the secondary flow could not be The experiment of Rozovskii (1961) [2] was used by Duan (2004) [50] to apply the depth-averaged 2D numerical model by including the dispersion terms 2D numerical model in a sharply curved open channel. It consisted of a 180 • curved reach with a 6-mlong straight approach and a 3-m-long straight exit. The ratio of mean radius of curvature to width was 1.0. The width of the channel was 1.7 m. The cross-section of the channel was rectangular. Water depth at the downstream end was 0.053 m, and the discharge was 0.0123 m 3 /s. The surface elevation and velocity distribution are shown in Figure 7 below. There is an acceleration in the longitudinal flow close to the outer bank of the upstream part of the bend and a deceleration in the flow velocity near the outer bank. The surface elevation increased at the outer bank and decreased at the inner bank. There was good agreement between the simulated and measured results. The impact of secondary flow was observed on the depth-averaged flow distribution and water surface elevations, which became more visible with the increasing of the channel curvature. Due to the transverse bed slope in natural rivers, the velocity of secondary flow near the bed is higher. There is a significant role of the secondary flow in transporting bedload and suspended load in the transverse direction (Duan et al., 2001) [59]. Obviously, a 3D model is more accurate than that of a depth-averaged model. However, the depth-averaged two-dimensional model with the dispersion terms according to this study was capable of simulating the hydrodynamic flow in meandering channels. As the secondary flow could not be taken into account in the depth-averaged turbulence model, the Schmidt number that correlates mass dispersion with turbulent diffusion can be adjusted as a calibration parameter. Finally, the accuracy of the 2D approach is still limited according to the results, and more measurements for calibration and verification are needed to obtain more accurate results. The channel of a 90° bend was used by Ghaneeizad et al. (2010) [51] to study the flow in a curved channel. The width of this channel was 0.4 m, the discharge was 13.55 lit/s, and the depth was 9 cm. The flow is subcritical, and the velocity on subcritical flows is at a high range. The boundary conditions were as follows: pressure-inlet, pressure-outlet, symmetry, and wall for the inlet, outlet, upper surface, and channel wall, respectively. Some turbulence models are used to simulate the flow in this channel including k-ω, k-ε (realizable), k-ɛ (RNG), and LES. Then, the simulated results are compared with the measured ones, and the ability of each model is examined with respect to flow analysis. Accordingly, the k-ε (RNG) model showed good performance in simulating the flow in the bend, and k-ω and k-ε (realizable) did not noticeably differ from the k-ε (RNG) model. However, LES has not shown its capability in simulating the flow in channel bends, especially after the curve (Ghaneeizad et al., 2010) [51].
The same channel of a 90° bend used by Ghaneeizad et al. (2010) [51] was also used by Gholami et al. (2014) [52] to study the flow patterns in a strongly curved open channel bend via the k-ɛ (RNG) model to predict the turbulence, and the volume of fluid (VOF) method to simulate the water's free surface. The channel is shown in Figure 8 below. The inlet-water boundary was set as a velocity inlet. The free surface boundary was set as the pressure-inlet condition on which atmospheric pressure was assumed. The boundary at the downstream outlet was defined as a pressure-outlet condition, and wall functions were used to estimate the effect of walls. The water surface is uniform before and after the bend, while along the bend and because of the centrifugal force, a transversal slope occurs The channel of a 90 • bend was used by Ghaneeizad et al. (2010) [51] to study the flow in a curved channel. The width of this channel was 0.4 m, the discharge was 13.55 lit/s, and the depth was 9 cm. The flow is subcritical, and the velocity on subcritical flows is at a high range. The boundary conditions were as follows: pressure-inlet, pressure-outlet, symmetry, and wall for the inlet, outlet, upper surface, and channel wall, respectively. Some turbulence models are used to simulate the flow in this channel including k-ω, k-ε (realizable), k-ε (RNG), and LES. Then, the simulated results are compared with the measured ones, and the ability of each model is examined with respect to flow analysis. Accordingly, the k-ε (RNG) model showed good performance in simulating the flow in the bend, and k-ω and k-ε (realizable) did not noticeably differ from the k-ε (RNG) model. However, LES has not shown its capability in simulating the flow in channel bends, especially after the curve (Ghaneeizad et al., 2010) [51].
The same channel of a 90 • bend used by Ghaneeizad et al. (2010) [51] was also used by Gholami et al. (2014) [52] to study the flow patterns in a strongly curved open channel bend via the k-ε (RNG) model to predict the turbulence, and the volume of fluid (VOF) method to simulate the water's free surface. The channel is shown in Figure 8 below. The inlet-water boundary was set as a velocity inlet. The free surface boundary was set as the pressure-inlet condition on which atmospheric pressure was assumed. The boundary at the downstream outlet was defined as a pressure-outlet condition, and wall functions were used to estimate the effect of walls. The water surface is uniform before and after the bend, while along the bend and because of the centrifugal force, a transversal slope occurs on the water's surface. The higher velocity appeared near the inner bank along the bend while the lower velocity appeared near the outer bank. This pattern was reversed after the bend. The results of the comparison between the numerical model and the experimental data indicated that there is an acceptable level of agreement. Major secondary flows and minor ones (near the water surface and in the outer wall) with two opposite rotation directions are formed in the numerical model. The flow pattern is impacted by these two secondary flows, which stay after the bend of the channel for some distance. As the flow in the bend is affected by secondary flows, the velocity increases accordingly in the area close to the bed, more than that near the water surface. Additionally, the higher velocity depth gradations cause an increase in the shear stress. close to the bed, more than that near the water surface. Additionally, the higher velocity depth gradations cause an increase in the shear stress. The flume data by Siebert (1982) [60] were used by Stoesser et al. (2010) [53] to simulate the flow in meandering open channels by using steady Reynolds-Averaged Navier-Stokes equations (RANS) using an isotropic turbulence closure, and Large Eddy Simulation (LES). The Reynolds stress term was modeled by the standard k-ɛ turbulence model or by the standard k-ω turbulence model, respectively. The flume consisted of two successive 180° bends joined by a 0.5 m-long straight cross-over section. There was a 4 mlong straight inlet before the first bend and a 2 m-long outlet. The channel's section was a rectangle of 0.25 m in width. For the boundary conditions, the mass-flux was kept constant to obtain a pressure gradient between the upstream and downstream sections of the domain. For the bed and side walls, a smooth boundary with the no-slip condition was used, and for the water's surface, a fixed rigid lid with zero gradients for all variables was used. The results showed that the higher momentum moves from the inner bank of the channel at the entrance of the bend towards the outer bank near the bend exit. Both the LES and the RANS calculations agreed well with the measurements and there is no significant difference between them. The higher velocity appears near the inner bank of the channel at the beginning of the bend, and then it moves towards the outer bank due to the centripetal forces acting in the previous bend. The maximum velocity is predicted well by the LES under the free surface. Additionally, the moving of higher velocity towards the inner bank of the channel at the beginning of the bend is predicted well by both RANS simulations. However, the suppression of the maximum streamwise velocity below the water surface is absent. Both RANS simulation and the LES model agreed well with the experiment in predicting the primary cell and an outer bank secondary cell at the peak of the bend. However, the RANS code failed to predict the persistence of the outer bank cell until the exit of the bend. The results also showed a good agreement between RANS and LES in the comparison of the computed bed and side-wall shear stresses.
The same flume data by Siebert (1982) [60] were also used by  The flume data by Siebert (1982) [60] were used by Stoesser et al. (2010) [53] to simulate the flow in meandering open channels by using steady Reynolds-Averaged Navier-Stokes equations (RANS) using an isotropic turbulence closure, and Large Eddy Simulation (LES). The Reynolds stress term was modeled by the standard k-ε turbulence model or by the standard k-ω turbulence model, respectively. The flume consisted of two successive 180 • bends joined by a 0.5 m-long straight cross-over section. There was a 4 m-long straight inlet before the first bend and a 2 m-long outlet. The channel's section was a rectangle of 0.25 m in width. For the boundary conditions, the mass-flux was kept constant to obtain a pressure gradient between the upstream and downstream sections of the domain. For the bed and side walls, a smooth boundary with the no-slip condition was used, and for the water's surface, a fixed rigid lid with zero gradients for all variables was used. The results showed that the higher momentum moves from the inner bank of the channel at the entrance of the bend towards the outer bank near the bend exit. Both the LES and the RANS calculations agreed well with the measurements and there is no significant difference between them.
The higher velocity appears near the inner bank of the channel at the beginning of the bend, and then it moves towards the outer bank due to the centripetal forces acting in the previous bend. The maximum velocity is predicted well by the LES under the free surface. Additionally, the moving of higher velocity towards the inner bank of the channel at the beginning of the bend is predicted well by both RANS simulations. However, the suppression of the maximum streamwise velocity below the water surface is absent. Both RANS simulation and the LES model agreed well with the experiment in predicting the primary cell and an outer bank secondary cell at the peak of the bend. However, the RANS code failed to predict the persistence of the outer bank cell until the exit of the bend. The results also showed a good agreement between RANS and LES in the comparison of the computed bed and side-wall shear stresses.
The same flume data by Siebert (1982) [60] were also used by Fraga et al. (2012) [54] to simulate the flow in open channel bends. Two types of models were used in this simulation to solve the unsteady Reynolds-Averaged Navier-Stokes equations (RANS). The first one was an isotropic turbulence closure (standard k-ε model), and the second one was a nonlinear k-ε model which takes into account anisotropic effects. The boundary conditions were considered as: open boundaries, wall and free surface. The results were compared with the LES data from Stoesser et al. (2010) [53] who used the same flume to simulate the flow by using the steady Reynolds-Averaged Navier-Stokes equations (RANS) with an isotropic turbulence closure, and Large Eddy Simulation (LES). Both models were similar to the LES model, which showed the shift of high-energy momentum from the inner bank to the outer bank. The flow moved under the effect of centripetal forces towards the outer wall and momentum was concentrated there at the first bend. As the flow reaches the second bend, the higher momentum is still under the impact of centripetal forces and moving towards the inner bank, which represents the outer bank of the first bend. Then, along the second bend, the higher momentum moves towards the outer bank of the channel due to the centripetal forces. The streamwise velocity is affected by the secondary flow as well. The highest streamwise velocity is located at some depth under the free surface instead of on the surface itself. The higher primary velocity is located at some depth under the free surfaces of the outer bank because of the centripetal forces. All three models are in agreement regarding the typical features of meandering channels. Additionally, all three models showed that the counter-clockwise secondary cell was inherited from the first bend, while only the non-linear k-ε and LES models showed the counter-rotating secondary cell of the previous bend, which is coincident with the higher primary velocities. However, the non-linear k-ε model's main secondary cell shape clearly differentiates it from the LES model. As the flow reaches the second bend, the higher momentum is still under the impact of centripetal forces and moving towards the inner bank, which represents the outer bank of the first bend. Then, along the second bend, the higher momentum moves towards the outer bank of the channel due to the centripetal forces. The streamwise velocity is affected by the secondary flow as well. The highest streamwise velocity is located at some depth under the free surface instead of on the surface itself. The higher primary velocity is located at some depth under the free surfaces of the outer bank because of the centripetal forces. All three models are in agreement regarding the typical features of meandering channels. Additionally, all three models showed that the counter-clockwise secondary cell was inherited from the first bend, while only the non-linear k-ε and LES models showed the counter-rotating secondary cell of the previous bend, which is coincident with the higher primary velocities. However, the non-linear k-ε model's main secondary cell shape clearly differentiates it from the LES model.  Four boundary conditions were applied in this case that include: the inlet and the outlet of the channel, the walls (including side walls and bottom), and the surface of the channel (symmetry). The simulation was performed using the OpenFOAM model and pisoFoam solver, mostly used for incompressible and turbulent flows in closed channel cases (also referred to as the Rigid-lid Model, or RLM). Two turbulence models were used to perform the simulation of flow in closed river bends: the standard k-ε model and the realizable k-ε model. The distribution of the velocity is shown in the figure below. The velocity increases near the inner bank at the beginning of the bend, and then moves towards the outer bank at the end of the bend (Figure 10). Four boundary conditions were applied in this case that include: the inlet and the outlet of the channel, the walls (including side walls and bottom), and the surface of the channel (symmetry). The simulation was performed using the OpenFOAM model and pisoFoam solver, mostly used for incompressible and turbulent flows in closed channel cases (also referred to as the Rigid-lid Model, or RLM). Two turbulence models were used to perform the simulation of flow in closed river bends: the standard k-ε model and the realizable k-ε model. The distribution of the velocity is shown in the figure below. The velocity increases near the inner bank at the beginning of the bend, and then moves towards the outer bank at the end of the bend (Figure 10). In order to examine the flow velocity, cross-sections were taken in the bend as shown in Figure 11. The comparison between the experimental data and the two models is shown in the figure, where: U is the depth-averaged velocity magnitude; and , , and are the longitudinal, vertical, and lateral velocity averaged over the depth, respectively. Note that the depth-averaged velocity is represented by the vertical axis and is non-dimensionalized by , which is the downstream mean velocity. The horizontal axis represents the radial distance ( − )/ , non-dimensionalized by the channel width, where B is the channel width and r and are the radius of curvature for the outer and inner banks, respectively. Both models performed well in the flow simulation of the bend. In addition, the comparison in the figure showed that the two models performed well compared to other numerical models that were used by other researchers. However, near the bend exit, there is a slight variation in the two models, particularly at the inner bank, due to the change from a curved to a straight channel, which makes the flow features in this region complex and not easy to capture, showing the weaknesses of the two models in this area (Shaheed et al., 2018) [55]. In order to examine the flow velocity, cross-sections were taken in the bend as shown in Figure 11. The comparison between the experimental data and the two models is shown in the figure, where: U is the depth-averaged velocity magnitude; and u 1 , u 2 , and u 3 are the longitudinal, vertical, and lateral velocity averaged over the depth, respectively. Note that the depth-averaged velocity is represented by the vertical axis and is non-dimensionalized by U 1 , which is the downstream mean velocity. The horizontal axis represents the radial distance (r − r i )/B, non-dimensionalized by the channel width, where B is the channel width and r and r i are the radius of curvature for the outer and inner banks, respectively. Both models performed well in the flow simulation of the bend. In addition, the comparison in the figure showed that the two models performed well compared to other numerical models that were used by other researchers. However, near the bend exit, there is a slight variation in the two models, particularly at the inner bank, due to the change from a curved to a straight channel, which makes the flow features in this region complex and not easy to capture, showing the weaknesses of the two models in this area (Shaheed et al., 2018) [55].
Different kinds of RANS turbulence models including standard k-ε, Re-Normalization Group (RNG) k-ε, realizable k-ε, k-ω, and shear stress transport (SST) k − ω were used by Yan et al. (2020) [56] to simulate the flow in a high-curvature 135 • channel bend with an assessment of side slopes and its effects on the secondary flows. The channel was rectangular, and it included a straight inlet section of 12.19 m followed by a 135 • bend and then a straight exit section of 2.44 m. The width was 1 m, and the height was 0.9 m with vertical acrylic walls. The performances of different turbulence models were compared against experimental measurements, and the best performing model was used to perform additional computations to assess the impact of different side slopes on the flow distribution. In addition to the velocity profiles, the models were also compared with more details in terms of the streamlines of secondary flow at the section 60 • from the bend. Two circulation cells appeared in the observations-the main, bigger one that was rotating clockwise, and the smaller one that was rotating counterclockwise, which was also called the outer bank cell. The realizable k-ε model was found to be the best among the tested models, especially in capturing some small flow features, such as the outer bank cells.
There are other numerical simulation models in river bends that have been used in previous research studies. The effect of secondary flow on depth-averaged equations using a stress diffusion matrix was studied by Lien et al. (1999) [61]. They used a 180 • bend with a rigid bed to examine the flow pattern using a two-dimensional numerical model. Accordingly, the effect of secondary flow on the maximum velocity path along the channel was indicated. Different kinds of RANS turbulence models including standard k-ε, Re-Normalization Group (RNG) k-ε, realizable k-ε, k-ω, and shear stress transport (SST) k − ω were used by Yan et al. (2020) [56] to simulate the flow in a high-curvature 135° channel bend with an assessment of side slopes and its effects on the secondary flows. The channel was rectangular, and it included a straight inlet section of 12.19 m followed by a 135° bend and LES and RANS (k-ε model) were used by Booij (2003) [7] to simulate the flow pattern in a mild 180 • bend with a rectangular section and rotating annular flumes. He concluded that the main characteristics of secondary flow and main flow could be produced by RANS computations. However, there were no satisfactory reproductions of the complicated secondary flow fields that could be given by RANS for the flow in rotating annular flumes and river bends. The secondary and main flow were reproduced well by LES for the mildly curved river. The turbulence was reproduced reasonably by LES and there was a low level of deviations in the computed secondary flow pattern of the rotating annular flume from the measured flow pattern. Additionally, the outer bank cell was reproduced well by the LES computations, unlike the RANS computations. LES computations had a small sensitivity towards the exact choice of modeling parameters.
Three types of turbulence models, including standard k-ε, k-ω and RSM, were used to model the flow pattern in a 180 • bend by Safarzadeh Gandoshmin and Salehi Neishbouri (2005) [62]. They found that the flow pattern could be predicted reasonably by all three models within the direct channel. However, during the curve, the prediction of RSM was superior. In addition, the k-ω model could predict the rotation and the separation effect of flow more precisely than the standard k-ε model.
The currents and mass transport in curved channels were simulated via a 3D hydrodynamic model of free surface turbulent flows by Jian and McCorquadale (1998) [37]. They modified the standard k-ε model and applied it to two typical curved open channel flows: a 270 • channel bend with a sloped outer bank and a meandering channel with mass transport. Then, they compared the model results with available data and obtained generally good agreement. According to their simulation, the secondary flows of the single bend were much stronger than that of the meanders. The reason might be related to the opposite spiral movement generated by the alternate bends of the meandering channels. Additionally, there could be a large lateral mass transport in the curved channel because of the secondary flows. A 3D numerical model was applied by Lu et al. (2004) [63] to simulate the flow pattern in a channel with a 180 • bend using the standard k-ε turbulence model. The results of the comparison between the numerical model and experimental data showed acceptable agreement. However, the model was not able to predict the minor secondary rotational cell in the outer wall.
A numerical simulation of turbulent free-surface flow was presented by Bodnar and Prihoda (2006) [64] using the k-ω turbulence model to analyze the non-linearity of the water surface slope at sharp river bends. The predictions of the model for the turbulent quantities and free-surface position were qualitatively correct.
The flow pattern in two types of bends, including a sharp one with a 180 • bend and a mild one with a 270 • bend, was simulated using a two-dimensional depth-averaged model by Zhou et al. (2009) [65]. The simulation was performed with and without the consideration of secondary flow. The results showed a better agreement in the first state of the effect of the secondary flow.
The Detached Eddy Simulation (DES) model was used by Constantinescu et al. (2010) [66] to simulate the flow in a strongly curved channel with a 193 • bend. They concluded that the prediction of the secondary flow and velocity components distribution by the model was satisfactory.
A sharply curved open channel bend (width of channel =1.3 m, the radius of curvature =1.7 m) with a flat horizontal bed was used by Blanckaert in 2009 [29] to measure the secondary flow using an acoustic Doppler vertical profiler (ADVP). The main target of this study was to examine the impact of relative curvature (the ratio of the flow depth (H) to the radius of curvature (R)) on secondary flow. Blanckaert's results indicated that there was a strong increase in the turbulent kinetic energy in the curve.
Van Balen et al. used the RANS and LES models in 2010 [67] to carry out numerical simulations of the sharp curve studied by Blanckaert in 2009 [30]. They found that the main secondary flow was underestimated by the RANS model, while the friction loss was overestimated by the same model. According to van Balen et al. (2009) [22], the outer bank cell was not shown by the RANS model, while it was shown by the LES model. As a result, they concluded that secondary flow in curved channels was the result of both centrifugal effects and turbulence anisotropy.
A rectangular channel with two linked curves (channel width = 0.25 m, the radius of curvature = 1.0 m) was used by Stoesser et al. (2010) [52] to compare LES and RANS models with experimental results. They concluded that the generation of an outer bank cell could be predicted by the RANS model, which contrasts with observations by van Balen et al. (2009;2010) [22,66].
The flow in a 60 • compound meandering channel with seminatural cross-sections was studied by Jing et al. (2009) [68] experimentally and numerically. The 3D numerical turbulence model that was applied to simulate the flow was the Reynolds stress model (RSM). The water depth impacts on wall shear stresses and secondary flow were studied as well. The calculations of velocity fields, wall shear stresses, and Reynolds stresses were performed for a range of input conditions. The comparison results between the measurements and simulated velocity fields and Reynolds shear stresses were in a reasonably good agreement. This indicated that the complex flow could be successfully predicted by RSM.
Additionally, the Reynolds stress model (RSM) was used by Jing et al. (2011) [69] to investigate the structure of turbulent flow in meandering compound open channels with trapezoidal cross-sections. The water depth impacts on secondary flows were studied as well. The results of the simulated and measured data showed that there was a considerable difference between the secondary flows of the overbank flow and that of inbank flow. Additionally, the simulation indicated that the magnitude and direction of the secondary flows were affected significantly by the water depth. The comparison results between the numerical simulation and the experimental data for the calculated velocity fields and Reynolds shear stresses were in good agreement, indicating the capability of this model on simulating the complex turbulent flow in meandering compound open channels.

Secondary Flow and Pollutant Dispersion
One of the effects of secondary flow in river bends is pollutant dispersion. As the inflow enters the curve of the river, there is a redistribution of flow momentum that occurs because of the secondary flow. As a result, there is considerable lateral mixing of pollutants. Thus, the pollutant spreading in curved rivers is much stronger than that in a straight river. Additionally, the distribution of contaminants discharged at the centerline of the river along the cross-section is nonsymmetric in curved rivers. According to Holley and Abraham (1973) [70], the pollutant dispersion at the inner bank is less than at the outer bank. Table 2 summarizes the relevant research used in this section and their key findings. Table 2. Summary of studies about pollutant dispersion in curved channels used in Section 5.2.

Shape of Channel Numerical Model Key Findings Reference
The flume consists of a single meander of two 90 • bends in alternating directions

Depth-averaged 2D numerical model
There is an agreement between the simulated and measured results. (Duan, 2004) [50] Three different meander situations were used including wide, shallow channels with a smooth bed, and two narrow channels with one having a smooth and the other a rough bed.
Standard k-ε turbulence model The agreement was fairly good in the case of the smooth bed for the vertical profiles of longitudinal velocity, but it was not so good for the case of the rough bed. (Demuren and Rodi, 1986) [71] The experiments of Chang (1971) [17] were used by Duan (2004) [50] to apply the 2D numerical model and evaluate its ability to simulate mass transport in meandering channels. The flume consisted of a single meander of two 90 • bends in alternating directions connected by short tangents. The flume width was 2.34 m, the ratio of actual length to the straight length was 1.17, while the sinuosity in the natural meandering channel is usually more than 1.5 (Leopold and Wolman, 1957) [72]. The cross-section was rectangular with 20 of 32 nominally smooth bed and sidewalls. The convection and diffusion equation were solved to obtain the concentration of transported mass: where C is the concentration, D b and E b = rates of mass deposition and entrainment, respectively; and σ c = turbulent Schmidt number for mass diffusion, which represents the ratio of eddy viscosity to eddy diffusivity. In order to study the dispersion of a pollutant, a constant concentration of pollutant was introduced at the inner bank, the midpoint of an upstream cross-section, and the outer bank for three separate runs. The velocity vector field and the water surface elevation are shown in Figure 12 below for the dye injection at the midpoint of the channel. It can be observed that the higher velocity transfers from the left bank of the first bend towards the right bank of the second bend because of the secondary flow. At the inner bank, the simulated elevation of the water surface is much lower than that at the outer bank. The comparison between the simulated depth-averaged velocity and concentration distribution and the measured experimental results at the second bend is shown in Figure 13. It is clear that there is good agreement between the simulated and measured results. This can give the impression that this model is able to simulate the depth-averaged flow field in meandering channels. However, as the simulation of secondary flow by the depth-averaged model is limited, calibration was performed for the Schmidt number to achieve the results in the simulations (Duan, 2004) [50]. The flow and pollutant dispersion in meandering channels were calculated mathematically and numerically by Demuren and Rodi (1986) [71] taking into account the three-dimensionality of the flow and pollutant concentration fields. The new model was an extension of the Leschziner and Rodi (1979) [26] three-dimensional numerical model. The standard k-ε turbulence model was used to improve the numerical accuracy. Three different meander situations were used to apply the new model. The three cases include a wide, shallow channel (width-to-depth ratio B/h = 20) with a smooth bed, and two fairly narrow channels (B/h = 5) with one having a smooth bed and the other a rough bed. The boundary conditions of the model were used as solid walls (bed and banks) and the surface is treated as a "rigid lid". There was a large secondary motion eddy that appeared in the wide-flume case. The results of the new model for longitudinal and secondary velocity field compared to the measurement were generally good. The pattern of the secondary flow in the smooth narrow flume was more complex with several developing and decaying eddies. The agreement was fairly good in the case of the smooth bed for the vertical profiles of longitudinal velocity, but it was not so good for the case of the rough bed. However, the prediction of the development of the depth-averaged longitudinal profile was fairly good for all cases. The prediction of secondary flow was realistic in all cases as well. Lastly, the bed shear stress, which is the main reason for bed erosion, appeared at the areas close to the depth-averaged velocity distribution. The rate of mixing in general as proceeding the downstream from one bend to another was reproduced well. Particularly, the model predicted the mixing well, which was much faster in the central discharge case than the inner or the outer bank discharge cases, and this could be due to the diffusivity that is generally larger in the center of the channel compared to near the banks. The agreement between the predicted and measured pollutant concentration was generally very good.  The flow and pollutant dispersion in meandering channels were calculated m matically and numerically by Demuren and Rodi (1986) [71] taking into account the t dimensionality of the flow and pollutant concentration fields. The new model was a tension of the Leschziner and Rodi (1979) [26] three-dimensional numerical model standard k-ɛ turbulence model was used to improve the numerical accuracy. Three d ent meander situations were used to apply the new model. The three cases include a w shallow channel (width-to-depth ratio B/h = 20) with a smooth bed, and two fairly na channels (B/h = 5) with one having a smooth bed and the other a rough bed. The boun conditions of the model were used as solid walls (bed and banks) and the surface is tr as a "rigid lid". There was a large secondary motion eddy that appeared in the wide-f  The flow and pollutant dispersion in meandering channels were calculate matically and numerically by Demuren and Rodi (1986) [71] taking into account t dimensionality of the flow and pollutant concentration fields. The new model w tension of the Leschziner and Rodi (1979) [26] three-dimensional numerical mo standard k-ɛ turbulence model was used to improve the numerical accuracy. Thr ent meander situations were used to apply the new model. The three cases includ shallow channel (width-to-depth ratio B/h = 20) with a smooth bed, and two fairl channels (B/h = 5) with one having a smooth bed and the other a rough bed. The b conditions of the model were used as solid walls (bed and banks) and the surface Figure 13. Comparison of simulated and measured depth-averaged velocity and concentration field (Reproduced from Duan (2004) [50] with permission from ASCE).

Secondary Flow and Sediment Transport
Another effect of secondary flow in river bends is sediment transport. Table 3 summarizes the research studied in this section and their key findings. Table 3. Summary of studies about sediment transport in curved channels used in Section 5.3.

Shape of Channel Numerical Model Key Findings Reference
A channel in 90 • bend Standard k-ε model and the low-Reynolds number k-ω model The low-Reynolds k-ω turbulence model is better than the standard k-ε model (Khosronejad et al., 2007) [20] A channel in 135 • bend Standard k-ε model and the low-Reynolds number k-ω model The low-Reynolds k-ω turbulence model is better than the standard k-ε model (Khosronejad et al., 2007) [20] A channel in 180 • bend k-ε turbulence model There is good agreement between the measurements and the predicted results (Wu et al., 2000) [40] Khosronejad et al. (2007) [20] used the channels of Matsuura (2004) [73] to apply the standard k-ε model and the low-Reynolds number k-ω model in the simulations of two different river bends as follows: The first channel of Matsuura (2004) [73] consisted of a rectangular open channel with a 90 • bend. The length of the channel was 20 m, the bed slope was 0.001, the width was 46 cm, the depth was 25.4 cm, and the radius of the bend was 672 mm. The bend entrance velocity was 0.285 m/s, and the uniform flow depth was 10.2 cm. Matsuura's channel is shown in Figure 14 below:

Secondary Flow and Sediment Transport
Another effect of secondary flow in river bends is sediment transport. Table marizes the research studied in this section and their key findings. Table 3. Summary of studies about sediment transport in curved channels used in Section 5.3.

A channel in 90° bend
Standard k-ε model and the low-Reynolds number k-ω model The low-Reynolds k-ω turbulence model is better than the standard k-ε model (Khosronejad et al., [20] A channel in 135° bend Standard k-ε model and the low-Reynolds number k-ω model The low-Reynolds k-ω turbulence model is better than the standard k-ε model (Khosronejad et al., [20] A channel in 180° bend k-ε turbulence model There is good agreement between the measurements and the pre-  [20] used the channels of Matsuura (2004) [73] to ap standard k-ε model and the low-Reynolds number k-ω model in the simulations different river bends as follows: The first channel of Matsuura (2004) [73] consisted of a rectangular open chan a 90° bend. The length of the channel was 20 m, the bed slope was 0.001, the width cm, the depth was 25.4 cm, and the radius of the bend was 672 mm. The bend e velocity was 0.285 m/s, and the uniform flow depth was 10.2 cm. Matsuura's ch shown in Figure 14 below: A significant helical flow was noticed in this area as well. This helical flow caused transverse sediment transport directed towards the inner bank of the bend. It is reported that the transverse sediment transport is the main cause of the degradation occurring near the outer banks and aggradation occurring near the inner banks of alluvial channel bends (Dietrich and Whiting 1989) [74]. The comparison between the simulated and measured bed topography showed that the low-Reynolds k-ω model fits the bed changes data better than the standard k-ε model. The mean error was also calculated. The results show that the mean absolute difference (M.A.D.) of the predicted from the observed bed levels was significant at 0.05 level of risk (α < 0.05) for the k-ω model (6.6 mm) and the k-ε model (8.2 mm). Additionally, the M.A.D. of the maximum observed scour depth of the k-ω and k-ε models was 11 and 13% and the standard deviations of absolute differences were 4.6 and 5.3 mm, respectively. Furthermore, standard deviations of absolute differences (S.D.A.D.) were 4.6 and 5.3 mm, and the coefficient of determination (r 2 ) for k-ω model (0.87) was relatively higher than the k-ε model (0.84). All these statistical measurements give an indication that the low-Reynolds k-ω turbulence model is better than the standard k-ε model. Based on this, it could be concluded that the k-ω model can predict the bed-shear stress and friction factor better than the k-ε model. The bank erosion was not considered because the walls were rigid in both physical and numerical models.
The second channel of Matsuura (2004) [73] consisted of a rectangular open channel with a 135 • bend. The inflow velocity was set at 0.283 m/s and the flow depth was 10.2 cm. The maximum scour depth observed (12.7 cm), which is higher than the 90 • case, was noticed along the outer bank at the bend exit ( Figure 15). The observations showed strong secondary flow and caused transverse sediment transport and aggradation toward the inner bank of the bend. The 135 • bend case was modeled considering the two turbulence models (standard k-ε and the low-Reynolds number k-ω) and the results of the simulation runs were compared to the corresponding observed results. significant at 0.05 level of risk (α < 0.05) for the k-ω model (6.6 mm) and the k-ɛ mo mm). Additionally, the M.A.D. of the maximum observed scour depth of the k-ω models was 11 and 13% and the standard deviations of absolute differences were 5.3 mm, respectively. Furthermore, standard deviations of absolute differences (S. were 4.6 and 5.3 mm, and the coefficient of determination (r 2 ) for k-ω model (0. relatively higher than the k-ɛ model (0.84). All these statistical measurements giv dication that the low-Reynolds k-ω turbulence model is better than the stand model. Based on this, it could be concluded that the k-ω model can predict the be stress and friction factor better than the k-ɛ model. The bank erosion was not con because the walls were rigid in both physical and numerical models.
The second channel of Matsuura (2004) [73] consisted of a rectangular open c with a 135° bend. The inflow velocity was set at 0.283 m/s and the flow depth was 1 The maximum scour depth observed (12.7 cm), which is higher than the 90° ca noticed along the outer bank at the bend exit ( Figure 15). The observations showed secondary flow and caused transverse sediment transport and aggradation tow inner bank of the bend. The 135° bend case was modeled considering the two tur models (standard k-ɛ and the low-Reynolds number k-ω) and the results of the sim runs were compared to the corresponding observed results. Figure 15. Geometric parameters and characteristics of 135° bend channel used in Matsuur [72] study (Reproduced from Khosronejad et al. (2007) [20] with permission from ASCE).
Similar to the 90° case, statistical analyses were conducted for comparison pu The mean absolute difference (M.A.D.) was found to be significant at a 0.05 leve nificance (α < 0.05) for the k-ω (11.2 mm) model and the k-ε model (13.3 mm). In ad Figure 15. Geometric parameters and characteristics of 135 • bend channel used in Matsuura 2004 [72] study (Reproduced from Khosronejad et al. (2007) [20] with permission from ASCE).
Similar to the 90 • case, statistical analyses were conducted for comparison purposes. The mean absolute difference (M.A.D.) was found to be significant at a 0.05 level of significance (α < 0.05) for the k-ω (11.2 mm) model and the k-ε model (13.3 mm). In addition, mean errors of the maximum scour depth for k-ω and k-ε models were 8.8 and 11%, respectively. The standard deviations of absolute differences (S.D.A.D.) and the coefficients of determination (r2) for k-ω and k-ε models were 7.4 and 8.7 mm, and 0.91 and 0.88, respectively. Based on the above, and compared with the results obtained from the 90 • bend case, it could be concluded that the maximum depth of scour for the 135 • bend was slightly underpredicted with both of the turbulence models, and the position of the maximum scour depth was not predicted as well as for the 90 • bend.
The flow and sediment transport are simulated by Wu et al. (2000) [40] in open channels through the use of a 3D numerical model. The 3D model was based on the general-purpose FAST3D flow solver developed at the University of Karlsruhe. The full Reynolds-averaged Navier-Stokes equations with the k-ε turbulence model were solved to calculate the flow. In order to validate the 3D numerical model, it was applied in a 180 • open channel bend that was studied experimentally by Odgaard and Bergs (1988) [75]. The channel was 80 m in length and 2.44 m in width and has a 180 • bend with 20 m long straight sections before and after the bend. The cross-section of the channel was trapezoidal with vertical side walls. At the beginning of the experiment, the channel was filled with a 30 cm thick layer of sand with an initially flat surface. The median diameter of the sand was d 50 = 0.3 mm and the geometric standard deviation was 1.45. The water depth, cross-stream bed profiles, and the streamwise development of the water depth were calculated and compared with the measurements obtained from the experiment. The results showed that there is good agreement between the compared values that can be observed with regard to the bed features. The calculated streamwise depth-averaged velocity was compared with measured velocities in the channel bend. The results showed good agreement with experimental data, taking into consideration the appearance of maximum velocities near the outer bank and the minimum velocities near the inner bank. This mainly could be a consequence of the water depth, which is larger near the outer bank than near the inner bank. The secondary flow motion was predicted well by the model. This could be seen by the appearance of positive angles near the surface that indicate motion to the outer bank and negative angles near the bottom motion toward the inner bank in the profiles of streamwise velocity vector. The model predicted the transverse bed slope and its downstream development well in general. There was good agreement between the measurements and the predictions regarding the bed form and the transverse slope. Dietrich et al. (1989) [74] investigated boundary shear stress, sediment transport, and bed morphology in a segment of a river during high and low flow. The segment was a sand-bedded river meander. They mapped the velocity, boundary shear stress, water surface, bedload, bedform, and suspended load transport fields during the spring snowmelt season for 1976-1979. They found that secondary flow changed the distribution of bed shear stress, which, in turn, affected the sediment transport processes. For example, Dietrich et al. (1989) [74] noticed that there was an outward shifting zone of maximum boundary shear stress through the bend due to the secondary flow, and the zone of high bedload transport clearly followed this zone. Jamieson et al. (2010) [76] investigated the spatial variability of three-dimensional Reynolds stresses in a developing channel bend. Experiments on the mean flow field and turbulence characteristics for flow in a mobile sand bed were conducted. In the experiments, the three components of instantaneous velocities at multiple cross-sections for different tests at different stages of clear water scour conditions were measured using Acoustic Doppler velocimeters. The experimental results showed that the secondary flow clearly altered the Reynolds stresses and, in turn, affects the sediment transport. For example, it was found that the near-bed maximum positive streamwise cross stream Reynolds stress coincides with the leading edge of the outer bank scour hole.
These studies demonstrated that the secondary flow could alter the boundary shear stresses and turbulence shear stresses, which, in turn, affects the sediment transport processes.

Secondary Flow in Bend over Topography
The research studied in this section focused on the secondary flow in bends over topography and they are summarized in Table 4. The laboratory flume data of Blanckaert (2003) [78] were used by Zeng et al. (2008) [24] to study the flow in sharp open channel bends. The flume consisted of a 9-m-long straight inflow channel reach followed by a 193 • bend with a constant centerline radius of curvature R = 1.7 m and a 5-m-long straight outflow reach. The flume was 22.7 m in length and 1.3 m in width with vertical lateral walls. The bed was covered with quasi-uniform sand having diameters in the range of 1.6-2.2 mm with an average of about d 50 = 2 mm. Two types of experiments were carried out: the first one with a flat bed and the second one with a mobile bed. The simulation of these two types of flumes was performed via the finite-differences RANS code to predict flow, sediment transport, and bathymetry in openchannel geometries with loose beds. Two turbulence models were used to provide the eddy viscosity: the k-ω Shear Stress Transport (SST) model and the Spalart-Allmaras (SA) model in low-Reynolds number versions. The main characteristics of the complex 3D flow field over a fixed bathymetry were accurately simulated by the RANS model, such as the interaction between the deformation of the streamwise velocities profiles and the strength of the cross-stream circulation cells, the production of turbulent kinetic energy, and others. The errors for the transverse velocity profiles were larger compared to the streamwise profiles. The RANS model also captured the main features of the flow field and the mobile bed. However, the deviations between model predictions and measurements increased. Accordingly, it could be considered that the accuracy and predictive capacity of the model was satisfactory. The major differences between the simulation and the measurement happen in the area close to the outer bank of the flume. The outer bank cell of the circulation measured in the flat bed experiment could not be predicted by the model in the simulation. In addition, a significant difference between the simulation and the measurement is observed for the turbulent kinetic energy in the equilibrium experiment. The turbulence models required in the simulation of both processes according to Blanckaert and De Vriend (2004) [79] are the ones that resolve turbulence anisotropy and the kinetic energy transfer between mean flow and turbulence.
The laboratory flume data of Blanckaert (2010) [80] were used by van Balen et al. (2010) [77] to apply the LES model to the simulation of a curved open-channel flow over an erodible bed. The flume had a movable sand bed and vertical sidewalls. The movable bed exhibited both the small-scale dunes on the bed and the large-scale point bar-pool structure of the bed along the bend. In the numerical model, the inflow section and outflow section are shortened to a length of 3.8 m each to save computational costs (see Figure 16).
The boundary conditions are treated as follows: free surface as a horizontal, impermeable rigid lid, solid walls with the wall-function approach, and for the outflow boundary, a convective boundary condition is imposed. When the geometry of the channel changed from straight to curved, the flow pattern along the bend will vary according to the location. In such a case, the distribution of the streamwise velocities at the entry of the bend is skewed toward the inner bank due to the presence of the sudden favorable longitudinal pressure gradient. The simulations were run for three separate LES computations, including a different roughness height of the bottom as shown in Table 5 below. The results showed that there are small differences between the three simulations regarding the primary flow structure. In addition, the flow structure of the three cases is fully determined by the large-scale bed topography. The flow decelerates in the upstream part of the bend, while the main flow accelerates approximately halfway in the bend and the return flow (backflow) exists near the inner bank. The flow in the downstream part of the bend recovers from the strong curvature-induced effects and the recirculation zone in the upstream part of the bend. Accordingly, there is a good qualitative agreement between the LES model and the experimental findings. However, the impact of the complex bed topography and sharp channel bend made the flow pattern very complex. The existence of the small-scale dune forms could be explained by quantitative differences. LES could not directly account for these small-scale dune forms because of the coarse measurements of the bed topography. There could be a spatial lag in the streamwise development of the flow due to the existence of dunes on the bed. The differences were marginal between the LES model and the RANS model as carried out by Zeng et al. (2008) [24] for the bigger part of the flow. This may be due to the turbulence-related momentum transport that does not play a role in the complete pattern of the momentum transport, and this is because of the combination of strong curvature and the presence of complex topography. Additionally, in the curvature, the turbulence stresses tend to be strongly isotropic. The boundary conditions are treated as follows: free surface as a horizontal, im meable rigid lid, solid walls with the wall-function approach, and for the outflow bo ary, a convective boundary condition is imposed. When the geometry of the cha changed from straight to curved, the flow pattern along the bend will vary accordin the location. In such a case, the distribution of the streamwise velocities at the entry o bend is skewed toward the inner bank due to the presence of the sudden favorable lo tudinal pressure gradient. The simulations were run for three separate LES computat including a different roughness height of the bottom as shown in Table 5 below. Th sults showed that there are small differences between the three simulations regarding primary flow structure. In addition, the flow structure of the three cases is fully d mined by the large-scale bed topography. The flow decelerates in the upstream part o bend, while the main flow accelerates approximately halfway in the bend and the re flow (backflow) exists near the inner bank. The flow in the downstream part of the b recovers from the strong curvature-induced effects and the recirculation zone in the stream part of the bend. Accordingly, there is a good qualitative agreement between LES model and the experimental findings. However, the impact of the complex be pography and sharp channel bend made the flow pattern very complex. The existen the small-scale dune forms could be explained by quantitative differences. LES could

Discussion
Numerical models are considered one of the important tools for predicting flow in a curved channel, allowing for various environmental studies, such as sediment transport and pollutant dispersion.
Different kinds of models such as standard k-ε, k-ε (RNG), k-ω, k-ε (Realizable), non-linear k-ε, LES and others were used previously by researchers for simulating the flow in different kinds of curved channels in order to study the secondary flow and its effects on these channels. It is not easy to determine a model that is best for all types of curved channels because each model has to specifically be the best in some cases of the curved channel, while it could not be able to capture the flow features accurately in other cases. Additionally, the model could be able to capture some characteristics of curved channels accurately, while its performance could be less reasonable for other features such as secondary cells and outer cells. However, LES models have been generally reported to be more accurate than most of the other tested models.
In river bends, some variables should be taken into consideration such as the radius and skewness of curvature, the curve's length, and the asymmetry of the channel crosssection. These variables might have a significant impact on the secondary flow, which, in turn, could affect the flow characteristics of one side, and other fluvial processes, of the river (mixing, transport of sediment, bank erosion, and meandering) as compared to the other side (Church et al., 2012) [81].
Uniform channel sections have been adopted in some studies to isolate the impact of the channel bend on the dynamics of the flow, therefore decreasing the number of problem dimensions (Booij, 2003 [7,22,29,77]. These studies shed light specifically on the formation of a second counterrotating secondary flow cell due to its significance in determining shear stress and the outside wall erosion on the curved channel. The choice of turbulence model depends on the type of secondary circulation being sought by numerical simulations. Generally, two types of secondary flow circulation can be observed in curved channels (Nezu and Nakagawa, 1993) [82]. The first type occurs due to the generation of centrifugal forces and pressure gradients. Reproduction of this type of secondary circulation in the numerical models is usually more feasible using both RANS and LES methods. The second type of secondary circulation occurs at the outer bank of the curve. This type is usually weaker than the first one and circulates in the opposite direction (Almeida and Ota, 2020) [83]. This type is formed by the non-homogeneity and anisotropy of Reynolds stresses (Nezu and Nakagawa, 1993;Blanckaert and De Vriend, 2004) [79,82]. According to Demuren and Rodi (1984) [84] and van Balen et al. (2009) [22], this type is more difficult to reproduce correctly by isotropic turbulent models. Moreover, Booij and Tukker (1996) [34] mentioned that there is no justification to use equal eddy viscosities for the different momentum transport terms in a curved channel.
Another important concept to be considered is the interaction of mean flow and turbulent fluctuations and energy transfer. The turbulence-induced near-bank circulation cells are usually hard to be represented by the linear turbulence closure models, despite the certain anisotropy of the cross-stream turbulence that they may well produce. On the one hand, the kinetic energy could be transferred from mean flow to turbulence and the mean flow vorticity is dissipated by the turbulent stresses. On the other hand, there is a transfer of the kinetic energy from turbulence to the mean flow as the experimental results indicated. However, comparing with the total production of turbulent kinetic energy, the amount of resituated kinetic energy is considered small. This could be mainly attributed to the boundary friction, which plays an important role in the outer bank cell dynamics. Therefore, in order to reproduce the outer bank cell accurately, turbulence closures that contain the ability to transfer kinetic energy between the turbulence and the average flow in either direction are required (Blanckaert and De Vriend, 2004) [79]. For example, some dynamic LES models can reproduce this two-way transfer. Indeed, LES models aim at capturing the scales that are large and contain energy in turbulent flow and modeling the interaction between the small and large scales. Most commonly used models assume that the essential function of the unresolved scales is to remove energy from the large scales. Then, the energy is dissipated through the action of viscous forces. While the energy on average is transmitted from large scales to the small ones, the inverse transfer is also recognized from small to large scales which could be quite important, and the employed LES model should be able to represent this property in order to accurately simulate near-bank circulation cells (Iliescu and Fischer, 2004) [85].

Future Research Needs
Although studying the secondary flow mechanisms in river bends and their relationship with other river processes has developed significantly, there are still some matters in this area that should be focused on and understood. Some aspects of secondary flow and the interrelation of flow characteristics and other fluvial processes of the river are discussed below: The relation between the secondary current and turbulence is one of the important aspects that should be better understood, especially how they actually depend on each other. Despite the significant advances in studying flows in meandering channels, demonstrated by recent publications, knowledge on the effects of mean flow three-dimensionality and secondary currents on turbulence in more complex flows is less complete. There are different facets of turbulence that should be further studied, and they can be studied by using several conceptual frameworks, such as the Reynolds-averaging framework, the coherent structures concept, and the eddy cascade concept.
The influence of secondary flows on hydraulic resistance is another important matter that should be taken into consideration. Generally, the transverse distributions of mean velocities are modified by the secondary flows. The fluid shear stresses, and bed shear stress are modified as well. The bulk friction factor increases according to these modifications as often assumed, compared to the same situation when the secondary flows are absent. Sometimes, the bulk friction factor is not affected by secondary currents, while the boundary shear stress and near-bank velocities change significantly (Kean et al., 2009) [86].
As it is known from previous studies, channel deformation, bank stability, and sediment transport are affected significantly by secondary flows, and the mechanisms and regular patterns require further investigation. In addition, suspended sediments are involved as a key factor of secondary current generation (Vanoni, 1946) [87]. The sediment transport affected by secondary flows in river bends and the bend curvature has important impacts on the sediment dynamics. The hydrodynamic features of rivers including secondary currents are highly influenced by river bends, and accordingly, cause a higher erosion power and enhanced sediment transport rates leading to increased channel migration rates (Church et al., 2012) [81].
Vertical, transverse, and longitudinal mixing are also highly affected by secondary flows (Rutherford, 1994) [88], and the effects should be investigated further. The mixing rate may be enhanced by secondary flows or dampened in all directions or selectively according to the specific flow configuration. For the case of curved channels, the longitudinal dispersion increases in the bend, which, in turn, reflects an increase in the turbulence intensity. On the other hand, the dispersion decreases in the curvature by enhancing transverse mixing. Accordingly, there is a highly efficient longitudinal dispersion at the beginning of the bend, and then it reduces greatly after the bend apex (Boxall et al., 2003;Marion and Zaramella, 2006;Rowinski. et al., 2008) [89][90][91].
Some other matters should be considered and focused on in future research studies with respect to flow in river bends. For example, the flow characteristics at more positions should be studied in order to closely examine the redistribution of the main flow and the evolution of the secondary flow. The velocity at the bottom and near the side walls should be studied carefully due to its significant impact on sediment transport. In numerical models, boundary conditions, roughness, and turbulence should receive more attention in river bends in order to gain more accurate results.
The outer bank cell plays an important role in the erosion of the outer bank and the mixing abilities in this region and needs to be well predicted. Thus, it is important to develop new techniques or models to better predict the outer bank cell.
Although many numerical models have been used to study the flow in river bends, the complexity of flow structure in this area, in addition to difficulty in prediction and understanding, indicate the need for further studies of numerical modeling, especially those that can accurately simulate the main secondary cell and the outer bank cell. A wider variety of numerical models such as LES, Detached Eddy Simulation (DES), Delayed Detached Eddy Simulation (DDES), and anisotropic RANS models is needed to be further assessed in simulations of the flow in river bends.

Concluding Remarks
Flow in river bends is obviously three-dimensional in nature, and different turbulence models have been used in the literature to simulate the flow at these areas. This study provides a comprehensive review on secondary flows in river bends, which is useful for a better understanding of the flow features in nature. The following conclusions can be drawn from the outputs of the simulations.
Theoretically, the results provided by a 3D model are more likely to outperform a depth-averaged (2D) model. However, 2D models might be more preferable due to their computational cost-effectiveness. A reasonable approximation was obtained from 2D models used in the study of mildly and sharply curved bends. In order to obtain more accurate results, a large number of measurements should be made available for calibration and verification.
The velocity distribution is significantly affected by the radius of the curvature. Higher velocity at the strongly curved bends occurs near the inner bank at the beginning of the bend, then it moves gradually towards the outer bank. The maximum velocity of the secondary velocities occurs below the free surface. The secondary velocities decrease as the radius of curvature increases. In addition, the higher velocity is displaced from the entrance towards the central part of the bend. The longitudinal velocities in strongly curved bends give a rough approximation of reality; moving the higher velocity from the inner bank of the bend towards the outer bank is predicted well. However, in reality, this occurs over a shorter distance. The longitudinal velocities in smoothly curved bends were predicted well in the central region and near the outer bank.
In the strongly curved bend, turbulence plays a significant role in the flow. The boundary conditions also play an important role in the distributions of the velocity in the cross-section. Choosing the correct and appropriate boundary conditions may lead to a more accurate simulation and better results.
Based on all studies reviewed in this paper, it can be seen that most of the RANS models and LES models performed well in simulating the main secondary cell in bends with some differences between them in accuracy. The differences could be attributed to some factors such as the geometry of the channel, curvature, roughness, etc., on the one hand, and the factors that were studied in the curved channel by these models such as type of the secondary cells, sediment transport, pollutant dispersion, etc., on the other hand. However, LES and related models could be recommended to use in order to better capture secondary flow structure, particularly the outer bank cell.