Suspended Sediment Transport Equation (Non-cohesive)

The sediment transport module calculates sediment transport induced by both suspended load and bedload. The morphological module calculates the bed evolution based on the sediment continuity equation.

Suspended sediment motion is governed by the depth-averaged sediment concentration equation as follows:

(63)\[(\bar{c} H)_t + \nabla_h \cdot (\bar{c} H ({\bf u}_\alpha + \bar{{\bf u} }_2)) =\nabla_h \cdot (k H (\nabla_h \bar{c})) + P - D \label{ad}\]

where \(\bar{c}\) is the non-dimensional depth-averaged sediment concentration normalized by sediment density. \(H(\bf{u}_\alpha + \bar{\bf{u}}_2) =M\) represents the flow rate per unit width defined in Shi et al. (2012), in which \(H=h+\eta\) is the total water depth. The roller-induced extra undertow can be taken into account as an option (see Wave Breaking, roller and undertow, and Physics (dispersion, breaking, friction)). \(k\) is the horizontal sediment diffusion coefficient calculated by the formula given by Elder (1959):

(64)\[k = 5.93 u_{*c} H\]

where \(u_{*c}\) is the shear velocity and can be calculated by van Rijn (1984):

(65)\[u_{*c} = \frac{\kappa}{-1 + \ln (30 H / k_s)} U_c\]

in which \(U_c\) is the depth-averaged total velocity (m/s), \(k_s = 2.5 d_{50}\) is Nikuradse roughness coefficient, \(d_{50}\) is the median grain diameter (mm), and \(\kappa\) is the von Karman constant.

In the advection-diffusion equation, \(P\) and \(D\) represent the erosion rate and deposition rate, respectively. The erosion rate can be calculated using van Rijn’s (1984) pickup function:

(66)\[P = 0.015 \frac{d_{50}}{a} \left ( \frac{|\tau_b| - \tau_{cr}}{\tau_{cr}}\right )^{1.5} d^{-0.3}_{*} w_f\]

where \(a\) is a reference elevation and is a function of total depth (\(a = 0.01 H\)), \(\tau_b\) is the bed shear stress, \(\tau_{cr}\) is the critical shear stress, and \(w_f\) is the settling velocity. \(P\) has the dimension of velocity (m/s) considering the convection-diffusion equation for non-dimensional sediment concentration. \(d_{*}\) is dimensionless grain size defined as:

(67)\[d_{*} = d_{50} \left( \frac{(s-1)g}{\nu^2} \right)^{1/3}\]

where \(s\) is the specific gravity of the sediment, and \(\nu\) is the kinematic viscosity coefficient. The critical bed shear stress \(\tau_{cr}\) used in (66) is defined as:

(68)\[\tau_{cr} = \rho_w (s-1)gd_{50} \theta_{cr}\]

where \(\theta_{cr}\) is the critical Shields parameter, approximately equal to 0.05. Based on the roughness estimate, the shear stress is expressed as:

(69)\[\tau_b = \rho_w \left(\frac{ 0.4}{1+\ln (k_s/30 h)} \right)^2 U_c^2\]

The deposition rate \(D\) can be calculated using the formula of Cao (1999):

(70)\[D = \gamma c w_f (1-\gamma \bar{c})^{m_o}\]

where \(\gamma = \min [2,(1-(1-n)/\bar{c})]\), \(n\) is the sediment porosity, and \(m_0\) is a constant number given as 2.0.


Cao, Z. (1999) “Equilibrium near-bed concentration of suspended sediment”. J. of Hydraulic Eng., 125(12): 1270-1278. DOI: 10.1061/(ASCE)0733-9429(1999)125:12(1270).

Elder, J.W. (1959). “The dispersion of marked fluid in turbulent shear flow”. J. of Fluid Mech., 5(4): 544-560. DOI: 10.1017/S0022112059000374.

Shi, F., J.T. Kirby, J.C. Harris, J.D. Geiman, and S.T. Grilli, (2012). “A high-order adaptive time-stepping TVD solver for Boussinesq modeling of breaking waves and coastal inundation”. Ocean Modelling, 43-44: 36-51. DOI: 10.1016/j.ocemod.2011.12.004.