**Simulated Flow Velocity Structure in Meandering Channels: Stratification and Inertia E**ff**ects Caused by Suspended Sediment**

#### **Fei Yang 1, Xuejun Shao 1, Xudong Fu 1,\* and Ehsan Kazemi <sup>2</sup>**


Received: 28 April 2019; Accepted: 11 June 2019; Published: 15 June 2019

**Abstract:** In this study, the coupled effects of sediment inertia and stratification on the pattern of secondary currents in bend-flows are evaluated using a 3D numerical model. The sediment inertia effect, as well as the stratification effect induced by the non-uniform distribution of suspended sediment, is accounted for by adopting the hydrodynamic equations without the Boussinesq approximation. The 3D model is validated by existing laboratory experimental results. Simulation results of a simplified meandering channel indicate that sediment stratification effect enhances the intensity of secondary flow via reducing eddy viscosity, while sediment inertia effect suppresses it. The integrated effects result in an increase and a reduction in the secondary flow, respectively, at lower and higher concentrations (near-bed volumetric concentrations of 0.015 and 0.1 are, respectively, considered in this study). This suggests that the dominance of the suspended sediment effect depends on the sediment concentration profile. With the increase of concentration under a specific sediment size, the secondary flow rises to reach a maximum, and then decreases. Moreover, as the sediment concentration increases, an exponentially decaying rate has been found for the secondary flow. It is concluded that in the numerical simulation of flow in meandering channels, when concentration is high, the variable-density hydrodynamic equations without the Boussinesq approximation should be considered.

**Keywords:** stratification effect; inertia effect; secondary flow; meandering; sediment laden flows

#### **1. Introduction**

In fluvial process studies, suspension has always been known to have a fundamental role in defining the preferential bend migration, since the stratification brought by suspended load modifies the turbulence, flow structure and sediment entrainment. Early works relate sediment stratification effect on flow velocity to concentration profiles. Vanoni [1] and Einsterin and Chien [2] classically documented the sediment stratification effect on the velocity profile in turbulent open channel flows. Vanoni showed that the von Karman constant κ reduces at high sediment concentration. This was based on the hypothesis that the suspended sediment dampens turbulence and increases velocity gradient near the wall. Later, many laboratory investigations [3–5] confirmed this finding. Castro-Orgaz et al. [6] analytically indicated that the reduction of κ in the high concentration condition is a physically coherent approximation. Suspension stratification effect modulates the turbulence intensity which in turn reduces the vertical sediment diffusivity. A few attempts referred the suspension effect to heat or salinity induced stratification effect [7]. By analogy with thermally stratified boundary layers, the Flux Richardson number or the gradient Richardson number for sediment-laden flow is defined as a measure

of stratification effect or turbulence damping [8,9] and the Monin-Obukhov length scale is introduced to establish damping functions for the eddy viscosity [10]. Herrmann and Madsen [11] obtained the optimal parameters for damping eddy diffusivity from observations. Nezu and Azuma [12] presented simultaneous measurements of both sediment and fluid in dilute flows. They pointed out that flow turbulence intensity is slightly dampened far from the bed while it rises in the near-wall region owing to the velocity difference between fluid and sediment. The ratio of concentration to slope is a primary indicator of stratification effect (which is greater in low slope rivers) according to Wright and Parker [13]. For high suspension concentration, Lamb et al. [9] experimentally observed that the Richardson number critical value of ~1/4 for completely damped turbulence [7] is not applicable where the major source of suspension is transported turbulence; and Muste et al. [14] provided the priority of two-phase over traditional mixed-flow perspective on sediment laden flows from image velocimetry experiments.

Although sediment laden flows have been widely studied by using various numerical models, most studies concern low suspension concentration using the advection-diffusion equation. The fundamental assumptions of these models limit their applicability. The advection-diffusion theory regulates the upward transfer of suspension due to turbulent diffusion with the downward settlement. Winterverp [15,16] numerically explained the sediment stratification effect by a buoyancy destruction term in the standard k-ε turbulence model at low concentration. Later, concerning hindered settling, van Maren et al. [17] extended this explanation to hyperconcentrated flows. Byun and Wang [18] coupled the sediment transport model and the 3D tidal hydrodynamics model to examine the reduction of turbulence and bottom shear stress due to the sediment stratification effect. Amoudry and Souza [19] performed 3D simulations and showed that sediment stratification effect leads to non-negligible modifications for bed evolution. Cantero et al. [20] exploited Direct Numerical Simulation (DNS)s of sediment laden flows under the dilute flow assumption and Boussinesq approximation. For higher settling velocity, the vertical Reynolds fluxes could be completely suppressed locally and relaminarization occurs [20,21]. For the condition of complete turbulence suppression, the critical product of the Richardson number and sediment settling velocity has a logarithmic dependence on the Reynolds number [22]. Dutta et al. [23] further investigated the sediment stratification effect on the vertical sediment diffusivity by DNSs. Dallali and Armenio [24] used Large Eddy Simulation (LES) to investigate suspension transport, supporting that the buoyancy effect on flow momentum is not neglectable for small particles. Alternatively, two-phase flow theory addresses water and sediment phases separately [25]. Recently, several models based the two-phase flow theory have been developed to include small scale mechanics [26–28]. The drift velocity shows that the inertia of fluid, particle turbulence and collisions effect contribute to sediment dispersion [26].

The Boussinesq approximation [29] is the basis for mathematical formulation of various flow problems. According to Kundu and Cohen [30], after subtracting a static reference state, i.e., ∇*p*<sup>0</sup> (= ρ<sup>0</sup> → *<sup>g</sup>*, where <sup>ρ</sup><sup>0</sup> and *<sup>p</sup>*<sup>0</sup> <sup>=</sup> reference mixture density and pressure, and <sup>→</sup> *g* = acceleration of gravity), the momentum equation can be expressed as

$$(1 + \frac{\rho'}{\rho\_0})\frac{\mathbf{d}\,\overrightarrow{\boldsymbol{u}}}{\mathbf{d}t} = -\frac{1}{\rho\_0}\nabla p' + \frac{\rho'}{\rho\_0}\overrightarrow{\mathbf{g}} + \nu\nabla^2\overrightarrow{\boldsymbol{u}}\tag{1}$$

where *t* = time; *p* and ρ = local deviations of pressure and mixture density from static reference state; → *u* = flow velocity vector; and ν = kinematic viscosity.

Both the inertia and buoyancy terms (the term on the left-hand side and the second term on the right-hand side of the above equation, respectively) contain the density ratio. With the Boussinesq approximation, the sediment stratification effect is shown as the buoyancy term while the inertia effect (the density ratio in the inertia term) is neglected. One then has

$$\frac{\mathbf{d}\overrightarrow{\boldsymbol{\mu}}}{\mathbf{d}t} = -\frac{1}{\rho\_0}\nabla p' + \frac{\rho'}{\rho\_0}\overrightarrow{\mathbf{g}} + \nu \nabla^2 \overrightarrow{\boldsymbol{\mu}}\tag{2}$$

Sediment concentration near the water surface is always much lower than that near the bed for the condition of heavy particles. A large vertical density gradient which has obvious inertia effect on flow can be a result of high bulk sediment concentration. The vertical concentration gradient is often very large with significant inertia variations. Therefore, the applicability of the Boussinesq approximation needs to be assessed in such a case. To the best of our knowledge, there have been no reported documents on the numerical comparison of sediment laden open-channel flows with and without considering the Boussinesq approximation.

Because the role of fluid inertia is highly important in the acceleration/deceleration processes of fluid flow, a 3D modelling method is adopted to thoroughly analyze the inertia effects of suspended sediment under complicated flow boundary conditions. The 3D model is substantially different from traditional hydrodynamic models under Boussinesq approximation which only consider the stratification effect of suspensions as the buoyance term. In the following sections, first, model setup is provided; second, a few cases are selected for model verification; third, numerical experiments with and without Boussinesq approximation are presented for comparison; and finally, discussions and suggestions with regards to the sediment inertia effect are given.

#### **2. Model Setup**

#### *2.1. Hydrodynamic Model*

The governing hydrodynamic equations are 3D shallow water equations based on the hydrostatic pressure. The mass and momentum equations are written as follows:

$$\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} + \frac{\partial w}{\partial z} = 0 \tag{3}$$

$$
\rho \frac{\mathrm{d}u}{\mathrm{d}t} = -\frac{\partial p}{\partial \mathbf{x}} + \frac{\partial}{\partial \mathbf{z}} \Big(\rho \nu\_{\mathbf{x}} \frac{\partial u}{\partial \mathbf{z}}\Big) + \frac{\partial}{\partial \mathbf{x}} \Big(\rho \nu\_{\mathbf{h}} \frac{\partial u}{\partial \mathbf{x}}\Big) + \frac{\partial}{\partial y} \Big(\rho \nu\_{\mathbf{h}} \frac{\partial u}{\partial y}\Big) \tag{4}$$

$$
\rho \frac{\mathrm{d}v}{\mathrm{d}t} = -\frac{\partial p}{\partial y} + \frac{\partial}{\partial z} \Big(\rho \nu\_x \frac{\partial v}{\partial z}\Big) + \frac{\partial}{\partial x} \Big(\rho \nu\_h \frac{\partial v}{\partial x}\Big) + \frac{\partial}{\partial y} \Big(\rho \nu\_h \frac{\partial v}{\partial y}\Big) \tag{5}
$$

$$-\rho \mathbf{g} + \frac{\partial p}{\partial z} = 0 \tag{6}$$

where *x*, *y* and *z* = Cartesian coordinates; *u*, *v* and *w* = components of velocity in the *x*, *y* and *z* directions; *p* = pressure; ρ = density of mixture; and ν<sup>h</sup> and ν<sup>z</sup> = eddy viscosity in the horizontal (*x* and *y*) and vertical (*z*) directions. The density of the mixture is determined by

$$
\rho = \rho\_w + (\rho\_s - \rho\_w)\mathfrak{c},
\tag{7}
$$

where *c* = volumetric concentration of sediment; ρ*<sup>w</sup>* = water density; and ρ<sup>s</sup> = sediment density.

Similar to layer-integrated models [31,32], the vertical velocity is calculated from the integrated continuity equation as follows

$$dw = w\_{z\_b} - \frac{\partial}{\partial x} \int\_{z\_b}^{z} \mu \mathrm{d}z - \frac{\partial}{\partial y} \int\_{z\_b}^{z} v \mathrm{d}z \tag{8}$$

where *z*<sup>b</sup> = bed elevation; and *wzb* is the vertical velocity at the bed. At the water surface,

$$\frac{\partial H}{\partial t} = w\_{z\_b} - \frac{\partial}{\partial x} \int\_{z\_b}^{H} \mu \mathrm{d}z - \frac{\partial}{\partial y} \int\_{z\_b}^{H} v \mathrm{d}z \tag{9}$$

where *H* = water surface elevation.

*Water* **2019**, *11*, 1254

The pressure in the horizontal momentum equations is given by integration of vertical momentum equation as

$$p = p\_d + g \int\_z^H \rho \,\mathrm{d}z \tag{10}$$

where *p*a = pressure at the water surface which is set to zero in the present simulations. The pressure gradient is obtained as follows

$$\frac{\partial p}{\partial y} = g\rho\_w \frac{\partial H}{\partial y} + g \frac{\partial}{\partial y} \int\_z^H (\rho - \rho\_w) \mathrm{d}z \tag{11}$$

Eddy diffusivity and viscosity are obtained from a two-equation *k-*ω turbulence closure model [33]. Complete extinction of local turbulence is unsustainable according to Cantero et al. [22]. To prevent the turbulence to be totally damped by the sediment stratification effect, a prescribed upper threshold 0.2 is set on the flux Richardson number. As the Reynolds-Averaged Navier–Stokes (RANS) model is incapable of modeling flow under conditions of strong stratification [34], an empirical distribution [25] for viscosity is reserved here as follows.

$$\nu\_{\rm h} = \nu\_z = \begin{cases} \kappa l L\_\* z (h - z) / h & y < 0.5h \\ \frac{1}{4} \kappa l L h & y \ge 0.5h \end{cases} \tag{12}$$

where *h* = water depth; and *U*\* = friction velocity.

#### *2.2. Sediment Transport Model*

Sediment transport is represented by the advection-diffusion equation,

$$\frac{\partial \mathbf{d}}{\partial t} = -\frac{\partial \mathbf{c}w\_{\mathbf{f}}}{\partial y} + \frac{\partial}{\partial z} \left( \varepsilon\_{x} \frac{\partial \mathbf{c}}{\partial z} \right) + \frac{\partial}{\partial \mathbf{x}} \left( \varepsilon\_{\mathbf{h}} \frac{\partial \mathbf{c}}{\partial \mathbf{x}} \right) + \frac{\partial}{\partial y} \left( \varepsilon\_{\mathbf{h}} \frac{\partial \mathbf{c}}{\partial y} \right) \tag{13}$$

where *w*<sup>f</sup> = sediment settling velocity; and ε<sup>h</sup> and ε<sup>z</sup> = eddy diffusivity in the horizontal (*x* and *y*) and vertical (*z*) directions.

#### *2.3. Boundary Condition*

There are two thin layers of δ at the water surface and the bed boundaries, thus the thickness of the calculation domain is *h* − 2δ. Vertical boundary conditions are implemented at these two interfaces. Normal gradients of horizontal velocities as well as the net vertical flux of sediment are zero at the upper interface. At the bottom one, we enforce the bottom frictional stress for momentum equations as

$$
\rho \nu\_x \left( \frac{\partial u}{\partial z'}, \frac{\partial v}{\partial z} \right)\_{z\_b} = \rho \text{C}\_D \sqrt{u\_b^2 + v\_b^2} (u\_{b\prime} v\_b) \tag{14}
$$

where *u*<sup>b</sup> and *v*<sup>b</sup> = near-bed components of horizontal velocity; and *CD* = drag coefficient to be calibrated. Besides, the near-bed net flux of sediment is

$$\left(\varepsilon\_{\rm x}\frac{\partial c}{\partial z} - cw\_{\rm f}\right)\_{z\_b} = q\_{\rm s
u} - c\_b w\_{\rm f} \tag{15}$$

where *q*su = sediment entrainment rate; and *c*<sup>b</sup> = near-bed volumetric concentration.

Flow condition at the inlet is given by flow discharge. In other words, the flow velocity profile at the inlet can be estimated through the transport equation for turbulent components and momentum equations with the assumption that flow is locally uniform at this boundary. Similarly, the suspended sediment concentration profile at the inlet is determined by assuming local equilibrium in the inlet sediment entrainment. At the outlet, the water level is set to prescribed values as the flow is assumed to be subcritical.

#### *2.4. Domain Discretization*

For discretization of the computational domain, a Spectral Element Method (SEM) is used in the vertical direction (due to the importance of the vertical variations in the present work) and the finite difference method is employed in the horizontal directions. Besides, for ease of implementation of the boundary conditions, vertical domain is discretized into three layers (due to the existence of two thin layers adjacent to water surface and bed) with unstretched z-coordinates. First, vertical space is divided into a number of sublayers with evenly thickness of *h*0. The top/bottom sublayer is combined with its adjacent sublayer, producing a layer with a thickness between *h*<sup>0</sup> and 2*h*0. Then, all the rest inner sublayers are combined to form an inner layer with its thickness being multiple of *h*0. Therefore, there will be three vertical layers throughout the whole calculation domain. All these three layers are variable according to water surface and bed changes during calculation process. Thickness thresholds *z*min and *z*max are imposed over top and bottom layers. When the thickness of the top/bottom layer is out of range (>*z*max or <*z*min), after the calculation of water surface/bed elevation at a certain time step, the layer interface will move a distance of *h*<sup>0</sup> inward or outward to ensure the thickness stays in the range. Let *h*<sup>0</sup> = (*z*max − *z*min)/2, after adjustment, its thickness will be close to (*z*max − *z*min)/2. Then, there is no need to readjust the thickness for small surface/bed changes. Every thickness adjustment is accompanied by local redistribution of physical quantities.

This kind of vertical discretization allows us to conveniently employ the finite dimensional space/polynomial of relatively higher degree for inner layer and lower degree for top/bottom layers. For the horizontal diffusion/viscosity calculations at the inner layer, if the horizontally adjacent inner layers are consistent, the calculations will be straight-forward as in the *z*-coordinates method, otherwise, additional interpolations need to be carried out separately for each calculation. For the calculations at the top/bottom layer, the horizontally adjacent layers are always inconsistent, and interpolations are required.

Horizontal domain is discretized using curvilinear coordinates into staggered grids. The governing equations need to be converted accordingly for implementation. Here, we present only the framework in the Cartesian coordinates for simplicity.

#### *2.5. Solution of Generalized Equation by SEM*

The horizontal momentum, sediment and turbulence transport equations can be generalized as

$$\begin{cases} \begin{array}{c} \frac{d\phi}{dt} = \frac{\partial}{\partial z} \Big( w\_l \phi + \frac{\nu\_l}{\sigma} \frac{\partial \phi}{\partial z} \Big) + f, z \in (z\_{b\*}H) \\ \frac{\nu\_l}{\sigma} \frac{\partial \phi}{\partial z} \Big|\_{H} = \mathcal{g}\_{1\prime} \cdot \frac{\nu\_l}{\sigma} \frac{\partial \phi}{\partial z} \Big|\_{z\_b} = \mathcal{g}\_2 \end{array} \tag{16}$$

where φ = the physical quantity to be solved; *w*<sup>f</sup> is set only for sediment transport; σ = the Prandtl–Schmidt number for φ; and *f* = other terms. The top and bottom boundaries are chosen according to the specific equation to be solved. At each vertical layer *i* with thickness *Li* = *zi*<sup>+</sup><sup>1</sup> − *zi*, generalized equation is conducted under local vertical coordinate ζ = 2 (*z* − *zi*)/*Li* − 1. We set an approximation such that φ belongs to the finite dimensional space of degree *N*. Choosing the Legendre polynomials as the local basis functions, the approximate solution is defined by

$$\phi = \sum\_{k=0}^{N} \phi\_k h\_k(\zeta) \text{ with, } h\_k(\zeta) = \frac{1}{N(N+1)L\_N(\zeta\_k)} \frac{(1-\zeta^2)L\_N'(\zeta)}{\zeta - \zeta\_k} \tag{17}$$

where φ*<sup>k</sup>* and *hk*(ζ) = the coefficients and the Lagrange form of Legendre polynomials, respectively; interpolation point ζ*<sup>k</sup>* = Legendre–Gauss–Lobatto point with its weight ω*k*; and *LN*(ζ) and *L <sup>N</sup>*(ζ) <sup>=</sup> the

Legendre polynomials and their derivatives. Generalized equation is multiplied by the test function and integrated over ζ (between [–1, 1]) and time step. We choose the local basis function as the test function. Then, the weak formulation after including boundary condition, can be derived as

$$\begin{split} & \frac{L\_i}{2} \sum\_{j=0}^N \phi\_j^{n+1} B\_{jk} + \frac{2}{L\_i} \Delta t \sum\_{m=0}^N \nu\_m \sum\_{j=0}^N \phi\_j^{n+1} D\_{jk,m} + \Delta t \sum\_{j=0}^N w\_l \phi\_j^{n+1} C\_{jk} \\ & = \frac{L\_i}{2} \sum\_{j=0}^N \phi\_j^n B\_{jk} + \frac{L\_i}{2} \Delta t \sum\_{j=0}^N f\_j^n B\_{jk}, \quad k = 0, \cdots, N. \end{split} \tag{18}$$

The integral coefficients are given as

$$B\_{jk} = \begin{cases} \omega\_j & j=k\\ o & j \neq k \end{cases}, \quad \mathcal{C}\_{j,k} = \frac{d h\_k(\zeta\_j)}{d \zeta} \omega\_{j\prime}, \; D\_{kj,m} = \frac{d h\_k(\zeta\_m)}{d \zeta} \frac{d h\_j(\zeta\_m)}{d \zeta} \omega\_m \tag{19}$$

where ω*<sup>j</sup>* = the weighting function of Gaussian quadrature for ζ*j*. The depth integrated continuity and horizontal momentum equations are solved in the coupled way.

Depth integrated continuity equation is discretized with a semi-implicit approach in curvilinear coordinates as

$$\begin{split} H^{n+1} &= z\_b + h^{n+1} \\ &= z\_b + h^n - f \Delta t \Big( \frac{\partial}{\partial \xi} \frac{h^n \overline{u}^n}{\overline{\phantom{v}}} + \frac{\partial}{\partial \eta} \frac{h^n \overline{u}^n}{\overline{\phantom{v}}} \Big) (1 - \theta) - f \Delta t \Big( \frac{\partial}{\partial \xi} \frac{h^n \overline{u}^{n+1}}{\overline{\phantom{v}}} + \frac{\partial}{\partial \eta} \frac{h^n \overline{u}^{n+1}}{\overline{\phantom{v}}} \Big) \theta \end{split} \tag{20}$$

where ξ and η = horizontal curvilinear coordinates; *u* and *v* = depth-averaged values of *u* and *v*; θ = the implicitness factor for temporal discretization; and *J* = Jacobi determinant. Substituting the discretized continuity equation into the horizontal momentum equations,

*Li* 2 -*N m*=0 *umn*+<sup>1</sup>*Bmk* + <sup>2</sup> *Li* Δ*t* -*N p*=0 ν*p* -*N m*=0 *umn*+<sup>1</sup>*Dmk*,*<sup>p</sup>* <sup>−</sup> <sup>Δ</sup>*<sup>t</sup>* 2 *Li* ν σ ∂*u* ∂ζ Γ2 −Δ*t* <sup>2</sup> *g Li* <sup>2</sup> <sup>β</sup>1θ<sup>2</sup> -*N m*=0 *Bmk* ∂ ∂ξ *J* ∂ ∂ξ *hnun*<sup>+</sup><sup>1</sup> *<sup>J</sup>* <sup>+</sup> *<sup>J</sup>* <sup>∂</sup> ∂η *hnvn*<sup>+</sup><sup>1</sup> *J* = *Li* 2 -*N m*=0 *um*<sup>∗</sup> *nBmk* + *Li* <sup>2</sup> Δ*t* -*N m*=0 *fmnBmk* <sup>−</sup>Δ*tg Li* <sup>2</sup> β<sup>1</sup> -*N m*=0 *Bmk* ∂ ∂ξ −Δ*t*θ(1 − θ)*J* ∂ ∂ξ *hnun <sup>J</sup>* <sup>+</sup> <sup>∂</sup> ∂η *hnvn J* + *H<sup>n</sup>* , *k* = 0, ··· , *N*. (21)

where β1=ξ<sup>2</sup> *<sup>x</sup>* + ξ<sup>2</sup> *<sup>y</sup>*; ν = ν*<sup>h</sup>* for simplification; σ = Schmidt number; and *um*\* *<sup>n</sup>* = the backtracked value of *u*m at time step *n*. Momentum equation for *v* has a similar form. Underlined terms are zero at the outlet for imposing water level condition for subcritical flows. In the momentum and sediment transport equations, advection terms are incorporated in the total derivatives. By backtracking characteristic lines accurately under Eulerian–Lagrangian context [35], the backtracked value at the foot of the characteristic lines at time step *n* can be approximated by interpolation. The final set of momentum equations can be solved by Conjugate Gradient. Equation (21) is semi-implicit with respect to *v*. Hence, *v* is firstly taken as an explicit quantity in the horizontal momentum for *u*, and then one more iteration is applied after the solution of the horizontal momentum equations. Once the horizontal velocities are known, water surface level and vertical velocity can be computed using the depth integrated continuity equation.

#### **3. Model Verification and Results**

The developed 3D numerical model is employed to simulate flow and suspended sediment transport in straight and meandering channels, and validated by comparing the measured and predicted distributions of flow hydrodynamics and sediment transport characteristics. To examine the effect of suspended sediment on the flow structure, a meandering channel flow with considerable suspended sediment is considered. As there was no experimental data available to us for suspended sediment motion in meandering channels, we firstly validated the model for the cases with suspended sediment at equilibrium and non-equilibrium conditions in straight channels with comparing the model results with a set of existing experimental data (Section 3.1); then applied the model for a meandering channel (but without presence of suspended sediment, i.e., only clear water) and validated the accuracy of the model against existing experiments (Section 3.2); and finally, numerically simulated flow and sediment motion in a meandering channel with investigating the effects of suspended sediment on the velocity structure (Section 3.3).

#### *3.1. Uniform and Steady Open Channel Flows*

A data set of laboratory experiment carried out by Coleman [36] was used to verify the model for flows with suspended sediment at equilibrium condition. The experiment was carried out in a 15 m long and 0.356 m wide recirculating flume with the slope of ~0.002. The measurements were taken at 12 m downstream of the inlet. The flow discharge and depth were ~0.064 m3/s and ~0.169 m, respectively. The velocity and sediment concentration were measured with a Pitot-static tube. Forty runs were performed with sediments of three different grain diameters and different concentrations. Runs 9, 20, 25, 31, and 40 are selected for comparisons with numerical results. Diameter of the sand grains were about 0.105, 0.210, and 0.420 mm, respectively, in runs {9 and 20}, {25 and 31}, and {40}. The Schmidt number for sediment vertical diffusivity is closely related to sediment diameter. Hence, we use the Schmidt numbers of 0.7, 1.2, and 2, respectively, for those runs.

Figure 1 presents the estimated velocity and sediment concentration profiles in comparison with the experimental ones. A good agreement between them is observed except in the upper part of the profiles, i.e., near the water surface. One reason could be that the effects of lateral walls on the flow are neglected in the present simulations, while those effects could be relatively high in the present case due to the low aspect ratio (width to depth) of the flume.

The model performance is tested for non-equilibrium condition, too. The non-equilibrium suspended sediment laboratory experiment of Wang and Ribberink [37] is considered. This experiment was carried out in a 30 m × 0.5 m × 0.5 m straight flume with a slope of 0.00097. The flow was 0.216 m deep and the discharge was 0.0601 m3/s. The flume was composed of a 10 m inflow section with a rigid bed, a 16 m test section with a perforated bed and a 4 m outflow section. Sediment was released at upstream section through the water surface. Nearly uniform sediment was used with settling velocity of *w*<sup>f</sup> = 0.007 m/s. The velocity was measured with a micropropellor and the sediment concentration was estimated by using the siphon method. The sediment profile was calculated at the inflow section as in an equilibrium state under given near-bed condition of *<sup>q</sup>*su <sup>=</sup> 1.19 <sup>×</sup> <sup>10</sup>−<sup>6</sup> <sup>m</sup>/s. In the test section, *q*su was set to zero. According to Wang and Ribberink [37], fluctuations were observed in the measured profile of the sediment concentration.

The numerical result is presented in Figure 2 in comparison with the experimental data. The comparison implies that the model performance is quite good for the non-equilibrium condition, too. The model is thus validated for the flow with suspended sediment at equilibrium and non-equilibrium conditions in straight channels.

**Figure 1.** Equilibrium condition: estimated profiles of velocity (**a**) and sediment concentration (**b**) in comparison with the experimental data.

**Figure 2.** Non-equilibrium condition: calculated sediment concentration profiles in comparison with the experimental data.

#### *3.2. Steady Flows in Meandering Channels*

The sharp bend experiment of Blanckaert [38] is used for model evaluation in curved bend flows. The experiment was carried out in a zero-slope flume consisted ofa9m straight inflow, a 193◦ curved bend and a 5 m straight outflow section (Figure 3). The curvature radius of the bend was R = 1.7 m. The width of the channel was B = 1.3 m, and the lateral walls were vertical and hydraulically smooth. The velocity was measured by an acoustic Doppler velocity profiler. For the present simulation, we selected one steady case with the water depth of 0.159 m and discharge of 0.089 m3/s.

**Figure 3.** Plan view of Blanckaert's laboratory flume.

The calculation was performed on a horizontal staggered grid with 130 × 20 cells. Calculation domain contains the curved bend and two 2.6 m length sections of inlet and outlet. The time step was set to 0.1 s. Inlet depth was identical to the measured one after adjusting the drag coefficient. Friction resistance is not taken into account for the side walls.

Figure 4a shows the distribution of the measured depth-averaged streamwise velocity. The short flow route leads to an acceleration close to the inner bank which is strengthened by transversal pressure gradient. At the same time, the streamline curvature induced a cross-sectional secondary flow, which shifts the larger velocity core toward the outer bank. These two mechanisms together generate the gradual outward shift of the larger velocity core at the channel bend. In the present calculation, magnitude of vertical viscosity is set to half of the horizontal one. Figure 4b shows the result of the numerical model. Outward shifting process of larger velocity core is closely consistent with the experiment. Figure 5 presents the streamwise velocity profiles in comparison with the experiments at two cross-sections shown on Figure 3, i.e., cross-sections 90◦ and 180◦. Although the present model is hydrostatic with isotropic turbulence model, only some small local differences with the experiment are observed near the side walls (Figure 5). Figure 6 presents the estimated transverse velocity profiles at the same cross-sections in comparison with the experimental data. It shows that the model is also capable of predicting transverse velocity with good accuracy.

**Figure 4.** Distribution of the depth-averaged streamwise velocity: (**a**) Experiment; (**b**) modelling.

**Figure 5.** Calculated streamwise velocity profiles in comparison with the measurements: (**a**) Cross-section 90◦; (**b**) cross-section 180◦. η is the position in the lateral direction from the central line towards the outer bank.

**Figure 6.** Calculated transverse velocity profiles in comparison with the measurements: (**a**) Cross-section 90◦; (**b**) cross-section 180◦. η is the position in the lateral direction from the central line towards the outer bank.

#### *3.3. E*ff*ects of Suspended Sediment on Velocity Structure*

Flows in curved open channels modulate the velocity components both in the streamwise and transverse directions as shown previously. Inertia of flow plays an important role in the acceleration and deceleration processes. In the dilute sediment laden flows, the mixture density can be approximately treated as a homogeneous quantity, thus the effect of sediment inertia can be neglected. Since sediment has a greater density than water, if concentration increases, how the sediment inertia affects the laden

flow is still unclear. To investigate this issue, here, we carry out a numerical study on the inertia effect in a meandering channel. Channel geometry and flow discharge in Section 3.2 are adopted. The parameters employed in the simulation are reported in Table 1. The simulations are carried out both with and without considering the Boussinesq approximation, with investigating the effect of near-bed volumetric concentration *c*b. Runs BAL and N0L are those with and without the Boussinesq approximation for *c*<sup>b</sup> of 0.015; and runs BAH and N0H are also with and without applying the Boussinesq approximation but for a higher concentration, i.e., *c*<sup>b</sup> of 0.1. Run CL, which is the clear water case presented in Section 3.2, is also included for comparison.


**Table 1.** Conditions of numerical simulations.

According to the results of the simulations, all the runs have similar water surface elevation distributions with different transversal slope. Here, we investigate the effect of suspended sediment on the flow structure. Figure 7 presents the results associated with the cross-section 90◦ for the lower concentration case (*c*<sup>b</sup> = 0.015), i.e., the runs BAL and N0L, in comparison with the clear water (CL) case. As can be seen, the streamwise velocity near the inner bank is smaller in the BAL and N0L cases compared to the CL case, while it is larger near the outer bank for BAL and N0L runs. That means the secondary flow is stronger in the BAL and N0L runs with respect to the CL case. Besides, it is clear that the suspended sediment has negligible inertia effect on the flow structure in present condition, since the velocity and sediment profiles almost overlap in runs BAL and N0L.

Existence of sediment with low concentration reduces the eddy viscosity so that the vertical gradient of the streamwise velocity is quite large close to the bed. In transversal direction, shear stress is balanced with outward centrifugal force and inward transversal pressure gradient. We can suppose that the direct impact of eddy viscosity on the shear stress is clearer than its indirect impact on the centrifugal force and pressure gradient. Lower eddy viscosity should cause larger transversal velocity gradient under nearly constant shear stress condition. Therefore, the transversal velocity will be larger, which leads to a stronger secondary flow.

Figure 8 shows the results for the higher concentration case (BAH and N0H) in comparison with the CL case. The secondary flow in the runs BAH and N0H seems weaker than the CL case. Similar to the lower concentration case, the eddy viscosity is reduced because of the sediment stratification effect, but the secondary flow is not strengthened in this case. This indicates that the sediment stratification effect on the suppression of the turbulence is not the key factor. As mentioned above, in the cases BAH and N0H, density is taken as a variable quantity in Equation (11). The buoyant force from the sediment stratification effect in the transversal momentum equation controls the secondary flow behavior. For example, consider the profile at centerline (η = 0) in Figure 8. The sediment concentration is larger close to the inner bank and transversal gradient of sediment concentration is very large. Transversal buoyant force direction is outward which is in opposite direction to the transversal surface gradient and the total transversal pressure gradient is smaller. The centrifugal force and the resultant transversal velocity are smaller than the clear water case. Run BAH and N0H have almost the same transversal surface gradient which has decisive role on the total transversal pressure gradient in the upper depth. Run N0H has smaller density in the upper depth compared to BAH. Therefore, in run N0H, the transversal velocity is larger at the upper part of the depth; and accordingly, the transversal surface gradient becomes larger, thus the transversal velocity in the lower depth increases in the reverse direction. Therefore, the secondary flow in run N0H is stronger than that in run BAH. This means the

effect of sediment inertia on the flow momentum is significant while the Boussinesq approximation cannot disclose the flow information at higher concentration in meandering channels.

**Figure 7.** Velocity and sediment concentration profiles at the cross-section 90◦ for lower concentration case: (**a**) Streamwise velocity; (**b**) transversal velocity; (**c**) sediment volumetric concentration. (velocity unit: m/s).

Runs N0H and BAH have higher concentration where sediment inertia effect is equally important as stratification effect. To illustrate the integrated effect on the secondary flow with different concentrations and sediment sizes, the sensitivity analysis on the size and concentration of sediment is discussed here. Two normalized quantities <*f* <sup>n</sup> <sup>2</sup>> and <*f* <sup>s</sup>*f* <sup>n</sup>> were adapted as characteristic parameters which represent the secondary flow strength and advecting flow momentum [39]. The sediment size is varied from 0.01 to 0.1 mm with depth averaged volumetric concentration from 0.01 to 0.2 at the inlet. The sediment buoyancy effect on turbulence is excluded here as the simple profile presented in Equation (12) for eddy viscosity was employed in this analysis. The results are summarized in Figure 9. The secondary flow strength decreases as the concentration increases. The finer sediment has a nearly uniform profile and only slightly affects the secondary flow strength. For 0.05- and 0.10-mm sediments, the non-uniformity of concentration profile at low concentration is much higher, vertical density gradient is large and the secondary flow strength declines efficiently to <sup>1</sup> <sup>4</sup> with concentration larger than 0.075. At higher concentration, the group hinder settling velocity is much smaller and the profile approaches to uniformity. Thus, the decrease retards and an exponentially decaying rate is found for the secondary flow here.

**Figure 8.** Velocity and sediment concentration profiles at the cross-section 180◦ for higher concentration case: (**a**) Streamwise velocity; (**b**) transversal velocity; (**c**) sediment volumetric concentration. (velocity unit: m/s).

**Figure 9.** <*f* n <sup>2</sup>> and <*f* <sup>s</sup>*f* <sup>n</sup>> variation with suspended sediment at the center of cross-section 90◦. The dot line represents the clear water (CL) case.

#### **4. Conclusions**

This paper presents 3D numerical simulations of sediment laden flows in open channel meanders concerning suspension inertia effect on the flow structure. The 3D governing equations without the Boussinesq approximation were used here and both the inertia and buoyant effects were considered. The velocity and concentration calculations performed in bend channels with comparison between the

cases with and without Boussinesq approximation. Sediment size and specific density of sediment were studied to demonstrate the integrated effect. Although the additional turbulence transport produced by the presence of sediment could not be simulated by the present model due to the limitations with the employed turbulence closure model, the modifications in the vertical profiles of concentration and velocity due to the sediment were still estimated successfully, similar to the previous studies.

The Boussinesq approximation that neglects all other density effects except buoyancy for dilute flows results in obvious discrepancy in flow structure and sediment transport for high sediment concentration in meandering channels. Both the inertia and stratification effects of sediment play important roles in river meandering processes. The stratification effect promotes secondary flow at lower concentration by reducing eddy viscosity and suppresses secondary flow at higher concentration by the buoyance force acting on transversal pressure. The inertia effect suppresses secondary flow which is important at higher concentration. Depending on the sediment concentration and diameter, the integrated effect results in an increase and a decrease of the secondary flow at lower and higher concentrations, respectively. This suggests that the dominance of the sediment effect depends on the concentration. With the increase of concentration under a specific sediment size, secondary flow increases to reach a maximum and then declines. Moreover, as the sediment concentration rises, an exponentially decaying rate has been found for the secondary flow subject to inertia effect (Figure 9). It is concluded that numerical simulations without considering the Boussinesq approximation is necessary when concentration is high in meandering channels.

The present qualitative predictions based on the 3D numerical simulations provide a new perspective on sediment laden flow. A quantitative analysis of the inertia and stratification performance in meandering channels provides a comprehensive understanding of the water-sediment interactions. The lower Yellow River carries a heavy sediment load with bed material mainly composed of fine sand and silt. The main stream channel migrates very fast during the occurrence of hyperconcentrated flows while the meandering mechanics is still an open question. Extensive research efforts, including numerical modelling, have been continuously carried out on meandering mechanisms. Highly efficient 2D models are still the main tools for river engineering and perform quite satisfactorily in straight or mild meandering rivers [40]. This study would give some cues to 2D and 3D modelling of fluvial process for sediment laden flows.

**Author Contributions:** Funding acquisition, writing—review, project administration and supervision, X.S. X.F. and E.K.; methodology, software, validation, and writing—original draft preparation, F.Y. and X.F.

**Funding:** This research was funded by National Key R&D Program of China (Grant No. 2018YFC0407402) and National Natural Science Foundation of China (NSFC, Grant Nos. 51525901 and 91747207).

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **References**


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