Suspended Sediment Concentrations

There are two common types of depth-averaged suspended sediment concentrations used in literature

1) \bar C _{sk}= \frac{1}{h-a} \int_{z_b +a}^{z_s} c_{sk}(z) dz
2) C _{sk}= \frac{1}{U(h-a)} \int_{z_b +a}^{z_s} u_s(z) c_{sk}(z) dz

where

h : water depth [L]
U : depth-averaged current velocity [L/T]
z_b : bed elevation [L]
z_s : water surface elevation [L]
a : reference elevation [L]
c_{sk} : local fraction suspended sediment concentration [M/L3]

There are advantages and disadvantages to both options.

Suspended-load Transport

The suspended load transport or unit discharge is defined as

3) q_{sk} = \beta_{sk} h U \bar C_{sk}
4) q_{sk} = h U C_{sk}

in which \beta_{sk} is referred to as the suspended-load correction factor (described below).

Suspended Sediment Concentration Profiles

HEC-RAS is a depth-averaged model, and does not explicitly solve for the suspended-sediment concentration profile. However, an approximation of the concentration profile is utilized in various calculations which help account for the non-uniformity of the sediment concentration profile. The parameters which utilize an approximation of the concentration profile are the load correction factor \beta_{sk} and the secondary flow correction function \chi_{sk} (described below). 

Various suspended-load concentration profiles can be derived by assuming different boundary conditions and distributions for the vertical sediment diffusion. The simplest of the profiles is the exponential profile which assumes a constant vertical diffusion coefficient

5) c_{sk} = c_{ska} \mathrm{e}^{-\phi_k \xi}

where

\xi = (z - z_b - a) / (h - a) : normalized vertical coordinate [0,1]
a : reference height [L]
z_b : bed elevation [L]
z_s : water surface elevation [L]
h : water depth [L]
\phi_k = \omega_{sk} h / \varepsilon_{vk} : non-dimensional vertical mixing parameter [-]
\omega_{sk} : sediment fall velocity [L/T]
\varepsilon_{svk} : vertical suspended sediment mixing coefficient [L2/T]
c_{ska} : suspended sediment concentration reference height a [M/L3]

The depth-averaged concentration can be computed for the exponential profile as

6) \bar C_{sk} = \frac{1}{h-a}\int_{z_b+a}^{z_s} c_{sk}(z)dz = \int_{0}^{1} c_{sk}(\xi)d \xi = \frac{c_{ska}}{\phi_k} \left( 1 - \mathrm{e}^{-\phi_k} \right)

The parameter \phi_k represents the ratio between the gravitational force and the vertical mixing.

An alternative common concentration profile can be derived by assuming a parabolic profile for the vertical mixing coefficient as

7) c_{sk} = c_{ska} \left( \frac{1/\xi - 1}{1/\xi_a - 1} \right)^{r_k}

in which

\xi : non-dimensional vertical coordinate [,1)
\xi_a : non-dimensional reference elevation [-]
r_k = \sigma_{sk} \omega_{sk} / (\kappa u_*)  : non-dimensional Rouse parameter [-]
\omega_{sk} : sediment fall velocity [L/T]
\sigma _{sk} : Schmidt number [-]
u_* : bed shear velocity [L/T]

The figure below shows a comparison of the exponential and Rouse concentration profiles.

Figure 1. Example non-dimensional concentration profiles as a function of the parameter \phi_k and r_k.

Vertical Mixing Coefficient

The vertical mixing coefficient is calculated as

8) \varepsilon_{svk} = \frac{\nu_{tv}}{\sigma_{sk}}
9) \nu_{tv} = c_{Mv} u_* h

where

c_{Mv} : empirical coefficient approximately equal to \kappa / 6 [-]
u_* : bed shear velocity [L/T]
h : water depth [L]
\sigma_k : Schmidt number [-]

Schmidt Number

The Schmidt number is calculated using the van Rijn (1984b) formula which is based on measurements from Coleman (1981)

10) \frac{1}{\sigma _{sk}}=1+2\left(\frac{\omega _{sk}}{u_{*}}\right)^{2}\,\,\,\,\,\mathrm{for} \,\,0.1 \lt \frac{\omega_{sk}}{u_{*}} \lt 1

where

ω_{sk} : sediment settling velocity  [L/T]
u_*: total bed shear velocity [L/T]

The applying the above equation the values of ω_{sk}/u_* values are limited to the range indicated above.

Suspended-load Transport Equation

The suspended-load transport equation may be written as

11) \frac{\partial }{\partial t} \left(\frac{h C_{sk}}{\beta _{sk}}\right) +\nabla \cdot \left( h \boldsymbol{U} C_{sk} \right) =\nabla \cdot \left(\varepsilon_{shk} h \nabla C_{sk} \right) + E_{sk} - D_{sk}

where

C_{sk}: suspended-load sediment concentration of the kth grain class [M/L3]
β_{sk}: suspended-load correction factor for the kth grain class [-]
U: depth-averaged current velocity in jth - direction [L/T]
h: water depth [L]
ε_{shk}: suspended-load horizontal diffusion (mixing) coefficient corresponding to the kth grain class [L2/T]
E_{sk} : suspended-load erosion rate [M/L2/T]
D_{sk}: suspended-load deposition rate [M/L2/T]

Suspended-load Correction Factor

The suspended-load correction factor is defined as the ration between the vertically integrated suspended-load transport and the transport computed as the simple product of the current velocity times the averaged sediment concentration:

12) \beta_{sk}=\frac{\int _{z_b+a}^{z_s}uc_{k}dz}{UhC_{k}}

For fine cohesive sediments, βsk is close to unity while for coarse sediments it is typically between 0.5 and 0.7. There are several options available for computing the suspended-load correction factor. The simplest is a user-specified constant. There are three additional options are based on assuming vertical profiles for the current velocity and sediment concentration.

By assuming logarithmic current velocity and exponential suspended sediment concentration profiles, an explicit expression for the suspended-load correction factor can be obtained as (modified from Sánchez and Wu 2011b)

13) \beta _{sk} = \frac{\mathrm{Ei}(-\phi_{k})-\mathrm{Ei}(-\phi_{k}A) + \mathrm{e}^{-\phi_k} \ln (A)}{\left[\ln(A)+1\right]\left(\mathrm{e}^{-\phi_k} -\mathrm{e}^{-\phi_kA} \right)}

where

\phi_k: non-dimensional vertical mixing parameter [-]
\omega_{sk} :
sediment fall velocity [L/T]
\varepsilon_{svk} :
vertical suspended sediment mixing coefficient [L2/T]
c_{ska} :
suspended sediment concentration reference height a [M/L3]

A = a/h [-]
Z = z_0/h [-]
u(z) : stream-wise current velocity [L/T]
a : reference height [L]
z_b : bed elevation [L]
z_s : water surface elevation [L]
z_a : apparent roughness length [L]
\mathrm{Ei}(x)=\int _{-\infty}^{x}\frac{e^{-t}}{t}dt : exponential integral

An analytical expression can be computed for the suspended-load correction assuming a power-law velocity and exponential concentration profile as

14) \beta _{sk} = \frac{\frac{m+1}{m}\gamma(\frac{m+1}{m},\phi_k)} {\phi_k^{1/m} \left( 1 - \mathrm{e}^{-\phi_k} \right)}

where

\gamma(s,x)=\int _{0}^{x}t^{s-1}e^{-t}dt : lower incomplete gamma function [-]
m = \frac{\kappa R^{1/6}}{n \sqrt g}: roughness parameter [-]
R : hydraulic radius [L]
g : gravity [L/T2]
n : Manning's roughness coefficient [T/L1/3]

Other velocity and suspended sediment concentration profiles can be easily integrated numerically to compute the suspended-load correction factor. The figure below shows a comparison of three different formulations. The different formulations are similar for 0.8 but vary more for smaller values of  \beta _{sk}.


Figure 2. Example comparison of the suspended-load correction factor for different streamwise and concentration profiles using a water depth of 10 m (32.8 ft); a) Power-law velocity and exponential concentration, b) Log-law velocity and exponential concentration, and c) log-law velocity and Rouse concentration.



Suspended-load Horizontal Diffusion Coefficient

The suspended-load horizontal mixing coefficient (εsk) represents the effects of turbulent diffusion. The horizontal sediment mixing coefficient is assumed to be related to the turbulent eddy viscosity as

15) \varepsilon_{svk} = \frac{\nu_{t}}{\sigma _{sk}}

where

σ_{sk}:Schmidt number for kth grain class [-]
ν_{t} = turbulent eddy viscosity [L2/T]

If a turbulent eddy viscosity is not available either because the flow model solves a Diffusion Wave equation or simply because it was ignored in the flow model, then it may be calculated as

16) ν_t = c_M u_* h

where

c_{M}: empirical coefficient (cM ≈ 0.5 − 6) [-]
u_{*}: bed shear velocity [L/T]
h: water depth [L]

It is noted that the horizontal and vertical sediment mixing coefficients have similar formulations but generally the horizontal coefficient will be much larger.