Secondary Flow

Introduction

In river bends, the balance between the outward centrifugal acceleration and the inward pressure gradient produces a secondary flow which flows towards the outer bank near the surface and towards the inner bank near the bed (Boussinesq 1868, Bowker 1988). The most accurate way of modeling the secondary currents at river bends is using a 3D model. However, the effects of the secondary flow can be approximated to some extent in depth-averaged 2D models by approximating the vertical current velocity profiles in the streamwise and transverse directions to estimate dispersive stresses (Duan 2004; Wu and Wang 2004 2005; Song et al. 2012). The streamwise velocity profile is typically assumed to be either power-law (e.g. . Odgaard 1986; Wu et al. 2004) or logarithmic (e.g. Duan 2004; De Vriend 1977). Several modeling approaches have been proposed to parameterize the vertical profile of the transverse current velocity which represents the secondary flow (e.g. Odgaard 1986; De Vriend 1977; 1981; Bernard and Schneider 1992). The most common vertical profile for the transverse current velocity is a linear one with zero mean (Odgaard 1986). The magnitude of the linear transverse velocity is parameterized as a function of the local flow curvature.

The purpose here is not to accurately resolve the 3D velocity field but rather estimate the simple and efficient relations for the momentum dispersion terms due to secondary flow with reasonable accuracy. The approach by Deltares (2024) is adopted due to it's simplicity and flexibility. Similar to the approach of Bernard and Schneider (1992), the Deltares (2024) approach only introduces momentum dispersive stresses due to flow curvature. This makes these approaches easier and faster to implement and calibrate in practice. The Deltares (2024) approach was favored here because it's simpler and easier to understand conceptually and to calibrate compared to the Bernard and Schneider (1992).

Secondary Flow Intensity

The secondary/spiral flow intensity, $\begin{array}{l}I\end{array}$ , is defined such that it is related to the transversal velocity at the free-surface by the simple expression

1)	$\begin{array}{l}\displaystyle \boldsymbol{u}_{ns} = \boldsymbol{n}_n b_s I\end{array}$

where

$\begin{array}{l}\boldsymbol{u}_{ns}\end{array}$ : surface transverse secondary flow velocity [L/T]
$\begin{array}{l}\boldsymbol{n}_n\end{array}$ : unit vector perpendicular to the streamwise direction [-1,1] [-]
$\begin{array}{l}b_s\end{array}$ : an empirical non-dimensional coefficient [-]
$\begin{array}{l}I\end{array}$ : secondary flow intensity [L/T]

The spiral flow intensity is defined using the right-hand-rule such that positive values indicate counter-clockwise rotations and negative values indicate clockwise rotations.

Various formulations have been used in literature for the coefficient . Here the formulation of Kalwijk and Booij (1986) is adopted:

2)	$\begin{array}{l}\displaystyle b_s = \frac{3}{\kappa^2} \left(\frac{1}{2} - \frac{1}{m} \right)\end{array}$

The secondary flow intensity can be either assumed to be in equilibrium or non-equilibrium. The non-equilibrium approach requires solving a transport equation and is therefore much more computationally expensive.

The intensity of the secondary flow can therefore be described by the magnitude of velocity magnitude at the surface and bed. In general, there are two approaches for estimating the secondary flow intensity: (1) Non-equilibrium, and (2) equilibrium. In the equilibrium approach, the secondary flow intensity is estimated based solely on the local hydrodynamic conditions, and the production and dissipation of secondary flow are assumed to be equal. In the non-equilibrium approach, a transport equation is used to simulate the unsteady production, dissipation, advection, and diffusion of the secondary flow intensity. Here two models are implemented which parameterize the secondary flow intensity as a function of a transport variable.

Equilibrium Approach

In the equilibrium approach the secondary flow intensity is computed simply as

3)	$\begin{array}{l}\displaystyle I = I_e = \beta_I h \|\boldsymbol{U}\| K_c\end{array}$

where

$\begin{array}{l}I\end{array}$ : secondary flow intensity [L/T]
$\begin{array}{l}I_e\end{array}$ : equilibrium secondary flow intensity [L/T]
$\begin{array}{l}\beta_I\end{array}$ : user-specified non-dimensional calibration parameter approximately equal to 1 [-]
$\begin{array}{l}h\end{array}$ : water depth [L]
$\begin{array}{l}K_c\end{array}$ : flow curvature [1/L]
$\begin{array}{l}\boldsymbol{U}\end{array}$ : depth-averaged current velocity [L/T]

Non-equilibrium Approach

In the non-equilibrium approach, the secondary flow intensity is calculated by solving the following transport equation

4)	$\begin{array}{l}\displaystyle \frac{\partial(h I)}{\partial t} + \nabla \cdot \left( h \boldsymbol{U} I \right) = \nabla \cdot \left( \varepsilon_{Ih} h \nabla I\right) + \frac{h \| \boldsymbol{U}\| }{L_a} \left( I_e - I \right)\end{array}$

5)	$\begin{array}{l}\displaystyle I_e = \beta_I h \|\boldsymbol{U}\| K_c\end{array}$

in which

$\begin{array}{l}I\end{array}$ : secondary flow intensity [L/T]
$\begin{array}{l}\beta_I\end{array}$ : secondary flow intensity scale factor [-]
$\begin{array}{l}h\end{array}$ : water depth [L]
$\begin{array}{l}L_a\end{array}$ : adaptation length [L]
$\begin{array}{l}\boldsymbol{U}\end{array}$ : current velocity magnitude [L/T]
$\begin{array}{l}\varepsilon_{Ih}\end{array}$ : secondary flow intensity horizontal diffusion coefficient [L²/T]
$\begin{array}{l}K_c = 1/R_c\end{array}$ : local flow curvature [1/L]
$\begin{array}{l}R_c\end{array}$ : local radius of curvature [L]

Secondary Flow Adaptation Length

The adaptation length is components as Kalkwijk and Booij (1986)

6)	$\begin{array}{l}\displaystyle L_a = f_L \frac{(m-1)h}{2\kappa^2}\end{array}$

where

$\begin{array}{l}f_L\end{array}$ : user-specified scaling factor (default of 1)
$\begin{array}{l}m\end{array}$ : non-dimensional roughness parameter [-]
$\begin{array}{l}h\end{array}$ : water depth [L]
$\begin{array}{l}\kappa\end{array}$ : von Karman constant of 0.41

Secondary Flow Diffusion Coefficient

The secondary flow intensity diffusion coefficient is scaled by the turbulent eddy viscosity as

7)	$\begin{array}{l}\displaystyle \varepsilon_{Ih} = f_D \nu_t\end{array}$

where

$\begin{array}{l}\varepsilon_{Ih}\end{array}$ : secondary flow intensity diffusion coefficient [L²/T]
$\begin{array}{l}\nu_t\end{array}$ : turbulent eddy viscosity [L²/T]
$\begin{array}{l}f_{D}\end{array}$ : scaling factor O(1) [-]

Flow Curvature

There are two approaches for computing the flow curvature: (1) based on the local velocity field, and (2) based on a user-specified Reference Line. By default HEC-RAS uses the local velocity field. The sign of the flow curvature indicates the direction in which the velocity is rotating along a streamline. HEC-RAS utilizes the right-hand convention in which positive flow curvatures indicate a counter-clockwise rotation, and negative flow curvatures indicate a clockwise rotation. With this convention, the flow curvature is positive for river bends which turn to the left and negative for river bends which turn to the right.

Flow Curvature Based on the Velocity Field

In this approach the horizontal curvature of flow streamlines is represented by the flow curvature. The local flow curvature is estimated using the local depth-averaged velocities and gradients as (Theisel 1995)

8)	$\begin{array}{l}\displaystyle K_c = \frac{1}{R_c} = \frac{U^2\dfrac{\partial V}{\partial x} + UV\left( \dfrac{\partial V}{\partial y} - \dfrac{\partial U}{\partial x} \right) -V^2\dfrac{\partial U}{\partial y}}{\|\boldsymbol{U}\|^3}\end{array}$

where

$\begin{array}{l}K_c\end{array}$ : flow curvature [1/L]
$\begin{array}{l}R_c\end{array}$ : radius of curvature [L/T]
$\begin{array}{l}\boldsymbol{U} = (U,V)^T\end{array}$ : depth-averaged velocity [L/T]

Flow Curvature Based on a User-Specified Reference Line

In this approach, the flow curvature is computed based on a user-specified centerline. This approach works well for rivers with a single and well-defined channel. It is not applicable to rivers with multiple channels or complex geomorphologies. The approach estimates the centerline radius of curvature based on the Circle Method. The radius of curvature at each cell is then computed as

9)	$\begin{array}{l}\displaystyle K_c(d) = \dfrac{K_c}{1 - d}\end{array}$

where $\begin{array}{l}K_c\end{array}$ is the curvature along the centerline, $\begin{array}{l}K_c(d)\end{array}$ is the curvature at a distance $\begin{array}{l}d\end{array}$ from the centerline. The sign of the distance is negative if the point is inside the curve and positive if the point is outside the curvature

Flow Curvature Limits

To improve stability the radius of curvature is limited as

10)	$\begin{array}{l}\displaystyle R_c > R_{c,\text{min}} = \min ( 10h,3 \sqrt{\Delta x} )\end{array}$

where

$\begin{array}{l}h\end{array}$ : water depth [L]
$\begin{array}{l}\Delta x\end{array}$ : local grid resolution estimated as the square root of the cell area [L]

Momentum Dispersion Stresses

The momentum dispersion stresses can be calculated analytically assuming utilizing the streamwise and transverse velocity profiles described above (e.g. Wu and Wang 20024; Wu 2007). These equations are provided in the Appendix. Here, the momentum dispersion stresses are calculated following the approach by Deltares (2024) as

11)	$\begin{array}{l}\displaystyle D_{xx} = - D_{yy} = 2 \rho \dfrac{UV}{\|\boldsymbol{U}\|} I \beta\end{array}$

12)	$\begin{array}{l}\displaystyle D_{xy} = D_{yx} = - \rho \dfrac{U^2 - V^2}{\|\boldsymbol{U}\|} I \beta\end{array}$

in which

$\begin{array}{l}D_{ij} = \boldsymbol{D}\end{array}$ : momentum dispersion stresses [L/T]
$\begin{array}{l}I\end{array}$ : secondary flow intensity [L/T]
$\begin{array}{l}\beta = \beta_c \left( 5\alpha - 15\alpha^2 + 37.5\alpha^3 \right)\end{array}$ (Deltares 2024)
$\begin{array}{l}\beta_c \in [0,1]\end{array}$ : user-defined calibration parameter [-]
$\begin{array}{l}\alpha = 1/m = \dfrac{ n \sqrt g}{\kappa R^{1/6}}\end{array}$ [-]
$\begin{array}{l}\kappa\end{array}$ : von Karman constant [-]
$\begin{array}{l}n\end{array}$ : Manning’s roughness coefficient [T/L^1/3]
$\begin{array}{l}g\end{array}$ : gravity [L/T²]
$\begin{array}{l}R\end{array}$ : hydraulic radius [L]

The Deltares (2024) approach utilizes a reduced set of terms. The method is preferred here because it is generally more robust than other approaches and can be more readily calibrated, especially as implemented here with scaling factors for the spiral intensity, adaptation length, and diffusion coefficient. The dispersion terms are such that when the flow curvature is zero, the dispersion terms also go to zero, making model calibration simpler and easier to understand.

Velocity Profiles

Although, the streamwise and transverse velocity profiles are not utilized for computing the dispersion stresses, the are utilizes for computing surface and near-bed velocities. The surface velocities are only utilized for output whereas the near-bed velocities are utilized for output and sediment transport calculations.

Streamwise Velocity Profile

The streamwise component of the velocity can be assumed to have a power-law profile as

13)	$\begin{array}{l}\displaystyle \frac{u_s(\xi)}{\|\boldsymbol{U}\|} = \frac{m+1}{m} \xi^{1/m}\end{array}$

where

$\begin{array}{l}u_s(\xi)\end{array}$ : streamwise velocity [L/T]
$\begin{array}{l}\xi\end{array}$ : non-dimensional vertical coordinate [0, 1] [-]
$\begin{array}{l}m\end{array}$ : non-dimensional roughness parameter (approximately equal to 7) [-]
$\begin{array}{l}|\boldsymbol{U}|\end{array}$ : depth-averaged velocity magnitude [L/T]

The roughness parameter $\begin{array}{l}m\end{array}$ can be calculated from the Manning’s roughness coefficient by

14)	$\begin{array}{l}\displaystyle m = \frac{\kappa R^{1/6}}{g \sqrt{g}}\end{array}$

in which

$\begin{array}{l}\kappa\end{array}$ : von Karman constant [-]
$\begin{array}{l}n\end{array}$ : Manning’s roughness coefficient [T/L^1/3]
$\begin{array}{l}g\end{array}$ : gravity [L/T²]
$\begin{array}{l}R\end{array}$ : hydraulic radius [L]

Transverse Velocity Profile

The transverse or spanwise velocity is approximated by the linear profile (Odgaard 1986)

15)	$\begin{array}{l}\displaystyle \frac{u_n(\xi)}{u_{ns}} = 2\xi - 1\end{array}$

where

$\begin{array}{l}\xi\end{array}$ : non-dimensional vertical coordinate between 0 and 1 [-]
$\begin{array}{l}u_{ns}=u_n(\xi=1)\end{array}$ : surface transversal (i.e. spiral or helical) velocity [L/T]
$\begin{array}{l}u_{nb}=u_n(\xi=0)\end{array}$ : near-bed transversal (i.e. spiral or helical) velocity [L/T]

The linear transverse velocity profile has been found to produce good results compared to laboratory measurements of current velocity profiles and surface water surface elevations (Wu and Wang 2004; 2005; Begnudelli et al. 2010; Lazzarin and Viero 2023).

Note on Bed Shear Stresses with Secondary Flow

The bed shear stresses can be modified to account for the effect of the secondary flow near-bed velocities. In practice however, it is not useful to apply to apply the secondary flow correction to the total bed shear stress because the momentum dispersion terms act as additional energy dissipation terms producing slightly higher water surface elevations (Qin et al. 2016). If the secondary flow correction is applied to the total bed shear stress, the overall friction losses tend to be over-estimated and the Manning’s roughness coefficient needs to be recalibrated. For this reason, the total bed shear stress is not modified to account for the secondary flow. However, the secondary flow correction is applied to the grain-related shear stress, which is used in computing the bed-load sediment transport potential.