Suspended Load - Hydro-Informatics

Suspended load refers to fine particle ( $\lesssim$ 1-2 mm) displacement in the water column, where particles are maintained in temporary suspension by the action of upward-moving turbulent eddies. The TELEMAC software suite uses the hydrodynamic Telemac2d/3d models to simulate Suspended load by solving the Advection-Diffusion equations with tracer concentrations. This is why suspended load modeling requires an open boundary LICBOR type for tracers (e.g., 4 or 5) as described in the setup of the boundaries-gaia.cli file.

To activate the simulation of suspended load, add the following to the Gaia steering file:

/ continued: gaia-morphodynamics.cas
/ ...
/ SUSPENDED LOAD
SUSPENSION FOR ALL SANDS : YES

Theoretical Background¶

The governing equation for suspended sediment transport is the advection-diffusion equation (ADE), which describes mass conservation of suspended sediment in the water column:

\frac{\partial (hC)}{\partial t} + \frac{\partial (hUC)}{\partial x} + \frac{\partial (hVC)}{\partial y} = \frac{\partial}{\partial x}\left(\varepsilon_s h \frac{\partial C}{\partial x}\right) + \frac{\partial}{\partial y}\left(\varepsilon_s h \frac{\partial C}{\partial y}\right) + E - D

(1)

where $C$ is the depth-averaged suspended sediment concentration (Gaia expresses it in g/l, numerically equal to kg m $^{-3}$ ), $h$ is water depth (m), $U$ and $V$ are depth-averaged velocity components (m s $^{-1}$ ), $\varepsilon_s$ is the sediment diffusivity coefficient (m $^2$ s $^{-1}$ ), $E$ is the erosion flux from the bed (kg m $^{-2}$ s $^{-1}$ ), and $D$ is the deposition flux to the bed (kg m $^{-2}$ s $^{-1}$ ).

The sediment diffusivity $\varepsilon_s$ is related to the turbulent eddy viscosity $\nu_t$ by:

\varepsilon_s = \frac{\nu_t}{\sigma_s}

(2)

where $\sigma_s$ is the Schmidt number, which Gaia fixes to $\sigma_s = 1.0$ (i.e., the sediment diffusivity equals the turbulent eddy viscosity). An additional constant diffusivity can be set with the COEFFICIENT FOR DIFFUSION OF SUSPENDED SEDIMENTS keyword (real, default 1.E-6 m $^2$ s $^{-1}$ ).

Additional Sediment Parameters¶

Fine sediment mixtures involving very fine cohesive particles (less than 0.06-0.1 mm) are referred to as mud in Gaia and so do the keywords in the following paragraphs. The distinction between non-cohesive sand and cohesive mud is important because their erosion and deposition behaviors differ fundamentally. More information about mud-related keywords can be found in section 4.2 in the Gaia manual.

Deposition Parameters¶

For suspended load, the definition of additional sediment properties for every sediment class is required (or enabled).

Particle settling velocities $w_{s}$ can be defined with the CLASSES SETTLING VELOCITIES keyword to calculate the deposition flux $D$ . The classical Krone (1962) deposition formula is:

D = w_{s} \cdot C \cdot \left(1 - \frac{\tau}{\tau_{cd}} \right) \quad \text{if } \tau < \tau_{cd}

(3)

where $C$ is the suspended sediment concentration (g/l), $\tau$ is the bed shear stress (N m $^{-2}$ ), and $\tau_{cd}$ is the critical shear stress for deposition (N m $^{-2}$ ). If $\tau \geq \tau_{cd}$ , no deposition occurs because turbulence is too strong to allow particles to settle.

If the CLASSES SETTLING VELOCITIES keyword is omitted (or set to -9), Gaia computes $w_s$ for each sediment class internally, selecting one of three grain-size-dependent formulae:

For very fine particles ( $D_{50} < 10^{-4}$ m), Stokes (1850) law applies:

w_{s} = \frac{(s-1) \cdot g \cdot D_{50}^2}{18 \nu}

(4)

For intermediate sizes ( $10^{-4} \leq D_{50} < 10^{-3}$ m), the Rubey--Zanke (1977) formula is used:

w_{s} = \frac{10\nu}{D_{50}}\left(\sqrt{1 + \frac{(s-1) \cdot g \cdot D_{50}^3}{100\nu^2}} - 1\right)

(5)

For coarse particles ( $D_{50} \geq 10^{-3}$ m), a constant drag-coefficient relation is used:

w_{s} = 1.1\sqrt{(s-1) \cdot g \cdot D_{50}}

(6)

where $s$ is the relative density of sediment (typically 2.65), $g$ is gravitational acceleration, $D_{50}$ is the grain diameter, and $\nu$ is the kinematic viscosity of water ( $\approx$ 10 $^{-6}$ m $^{2}$ s $^{-1}$ ). The three regimes transition from a viscous ( $Re_p \ll 1$ , Stokes) to a fully turbulent ( $Re_p \gg 1$ , constant drag) settling behavior Dey, 2014.

To take advantage of Gaia’s integrated routines for calculating $w_{s}$ , either do not use the CLASSES SETTLING VELOCITIES keyword in the Gaia steering file, or set its per-class values to -9 (which triggers automatic calculation). Detailed information on the calculation of settling velocities for particular cases (e.g., suspended load calculation for other suspended material than mineral sediment) can be found, for example, in Dey (2014) (book section 1.7). Gaia’s settling velocity algorithm is located in the file settling_vel.f in the /telemac/sources/gaia/ directory.

The critical shear stress $\tau_{cd}$ for mud deposition can be defined with the CLASSES CRITICAL SHEAR STRESS FOR MUD DEPOSITION keyword (default is 1000. N m $^{-2}$ , which effectively disables the deposition threshold; Gaia converts it internally to the critical shear velocity $u_{*,cd} = \sqrt{\tau_{cd}/\rho_w}$ ).

/ continued: gaia-morphodynamics.cas
/ ...
CLASSES SETTLING VELOCITIES : -9;-9;-9
CLASSES CRITICAL SHEAR STRESS FOR MUD DEPOSITION : 1000;1000;1000 / N per m2

Erosion Parameters¶

For cohesive (mud) sediments, Gaia calculates erosion fluxes $E$ using the Partheniades (1965) formula, which is the classical approach for cohesive sediments:

E = \begin{cases} M\cdot \left(\frac{\tau}{\tau_{ce}} - 1\right) & \mbox{ if } \tau > \tau_{ce} \\ 0 & \mbox{ if } \tau \leq \tau_{ce}\end{cases}

(8)

where $M$ denotes the Krone (1962)--Partheniades (1965) erosion constant (in kg m $^{-2}$ s $^{-1}$ ), which can be defined in Gaia with the LAYERS PARTHENIADES CONSTANT keyword (default value: 1.E-03). Moreover, $\tau_{ce}$ (critical shear stress for erosion) can be defined with the LAYERS CRITICAL EROSION SHEAR STRESS OF THE MUD keyword (default is 0.01;0.02;0.03;... for successive layers) in N m $^{-2}$ .

Non-cohesive sand uses an equilibrium-concentration closure

The Partheniades formula above applies to cohesive mud. For non-cohesive sand (the case used in this tutorial), Gaia does not use the Partheniades constant. Instead, the net bed exchange flux is computed from the equilibrium near-bed concentration $C_{eq}$ obtained from the chosen suspension formula following the Celik & Rodi (1988) approach: $E - D = w_s \, (C_{eq} - C_{z_{ref}})$ , where $C_{z_{ref}}$ is the actual near-bed concentration derived from the depth-averaged concentration assuming a Rouse (1939) profile. Erosion ( $E = w_s C_{eq}$ ) dominates when the bed is under-saturated, and deposition ( $D = w_s C_{z_{ref}}$ ) dominates when it is over-saturated.

/ continued: gaia-morphodynamics.cas
/ ...
LAYERS PARTHENIADES CONSTANT : 1.E-03 / in kg per m2 per s
/ LAYERS CRITICAL EROSION SHEAR STRESS OF THE MUD : 0.01;0.1;0.1 / in N per m2

Suspended Load Formulae¶

The sediment transport formulae for suspended load modeling can be defined with the SUSPENSION TRANSPORT FORMULA FOR ALL SANDS keyword, which accepts an integer number defining a formula for calculating the equilibrium near-bed concentration $C_{eq}$ in g/l (the unit Gaia uses internally for all suspended sediment concentrations). The equilibrium concentration represents the sediment concentration at a reference level near the bed under equilibrium conditions (i.e., when erosion equals deposition). The calculated $C_{eq}$ values align with the later definition of initial and boundary conditions for suspended load.

The following integer numbers can be used for calculating $C_{eq}$ with the SUSPENSION TRANSPORT FORMULA FOR ALL SANDS keyword:

1 for the Zyserman & Fredsøe (1994) formula (default and used in this tutorial):
- Empirical formula based on experimental data from Guy et al. (1966)
- Uses a skin friction correction (cf. bedload corrections) for the Shields parameter
- Applicable for non-cohesive sediments in fluvial environments
- Reference (near-bed) elevation $z_{ref} = \alpha_{k_s} \cdot D_{50}$ (default $3.0 \cdot D_{50}$ , modifiable with RATIO BETWEEN SKIN FRICTION AND MEAN DIAMETER)
- Defined in /telemac/sources/gaia/suspension_fredsoe.f
- Formula: $C_{eq} = \frac{0.331 \cdot (\theta' - \theta_{cr})^{1.75}}{1 + 0.72 \cdot (\theta' - \theta_{cr})^{1.75}}$ where $\theta' = \mu\theta$ is the skin-friction Shields parameter and $\theta_{cr}$ is the critical Shields parameter
2 for the Bijker (1992) formula:
- Calculates suspended load concentration as a function of bedload and a reference skin-friction elevation
- Requires that bedload calculation is activated (BED LOAD FOR ALL SANDS : YES)
- Suitable for combined bedload-suspended load calculations
- Reference elevation $z_{ref} = k_{sr}$ (the rippled bed roughness)
- Defined in /telemac/sources/gaia/suspension_bijker.f
3 for the Van Rijn (1984) formula:
- Counterpart of the van Rijn bedload formula
- Uses a skin friction correction (cf. bedload corrections) for the Shields parameter
- Reference elevation $z_{ref} = 0.5 \cdot k_s$ where $k_s$ is the total roughness (from the hydrodynamics steering file)
- Originally developed for sand transport in rivers and estuaries
- Defined in /telemac/sources/gaia/suspension_vanrijn.f
4 for the Soulsby (1997)-Van Rijn (2007) formula:
- Uses orbital velocity of waves (i.e., suggested application: coastal/marine regions)
- Combines current and wave effects on sediment suspension
- Read more about suspended load and waves in section 5.1 of the Gaia manual
- Defined in /telemac/sources/gaia/suspension_sandflow.f

/ continued: gaia-morphodynamics.cas
/ ...
SUSPENSION TRANSPORT FORMULA FOR ALL SANDS : 1

Initial and Boundary Conditions¶

Gaia enables a class-wise definition of initial concentrations for suspended load following the order of sediment class definitions. The following list definition sets the initial concentration for the 0.5-mm sediment class (recall its definition) to 0.6 g/l and 0.0 g/l for the 0.02-m and 0.1-m sediment size classes. The definition of initial suspended sediment concentrations can be overridden in 2d at boundary nodes by setting the EQUILIBRIUM INFLOW CONCENTRATION keyword to YES (requires that the tracer boundary is set to 5).

/ continued: gaia-morphodynamics.cas
/ ...
INITIAL SUSPENDED SEDIMENTS CONCENTRATION VALUES : 0.6;0.;0.

Read more about the definition of initial conditions in section 2.1.1 in the Gaia manual.

Boundary Prescriptions¶

The per-sediment class suspended load concentrations can be prescribed similar to the initial concentrations with the PRESCRIBED SUSPENDED SEDIMENTS CONCENTRATION VALUES keyword. Alternatively, the EQUILIBRIUM INFLOW CONCENTRATION keyword may be used to automatically compute the inflow concentration based on the equilibrium formula (option 1-4 defined above). None of these keywords is used in this tutorial because the model starts with a defined initial concentration and allows the system to evolve.

/ continued: gaia-morphodynamics.cas
/ ...
/ PRESCRIBED SUSPENDED SEDIMENTS CONCENTRATION VALUES : 0.6;0.;0. / g/l
/ EQUILIBRIUM INFLOW CONCENTRATION : YES / not used in this tutorial

Gaia can be run with liquid boundary files for assigning time-dependent suspended load fluxes (the outflow should be kept in equilibrium). Solid flux time series can be implemented using the already applied 455-5 upstream boundary type, analogous to the descriptions of the Telemac2d unsteady boundary setup. More information about suspended load boundary conditions can be found in section 2.1.2 in the Gaia manual.

Numerical Parameters¶

Most numerical parameters for suspended load modeling depend on hydrodynamic Telemac2d/3d steering file definitions. Additional keywords directly affecting the simulation of suspended load should be declared in the Gaia steering file.

For instance, the SCHEME FOR ADVECTION ... keywords for velocities, tracers, and turbulence modeling are defined with the hydrodynamics (Telemac2d/3d) steering file’s general numerical parameters for finite elements. In addition, the advection scheme for suspended load can be defined in the Gaia steering file with the SCHEME FOR ADVECTION OF SUSPENDED SEDIMENTS keyword that accepts one of the following integer keywords (for 2d only):

1 for the unconditionally stable, non-conservative, but diffusive (for small timesteps) Method of Characteristics scheme.
2 for the non-conservative Streamline Upwind Petrov Galerkin (SUPG) scheme that uses the CFL condition and is less diffusive than the Characteristics (1) scheme.
3 or 4 for the conservative N-scheme (distributive) with timestep reduction based on the CFL condition. Option 4 includes mass-lumping for improved stability. These options should not be used in the presence of tidal flats (use 13 or 14 instead).
5 for the mass-conservative PSI distributive scheme (default), which corrects fluxes according to tracer concentrations and is less diffusive than 4 or 14. Computation time with 5 is longer than with 4 or 14. This option should not be used in the presence of tidal flats.
13 and 14 for the Edge-based N-scheme (NERD), which is similar to 3 and 4, but adapted to tidal flats. Option 14 is used in this tutorial according to the recommendation in the Gaia manual.
15 for the mass-conservative ERIA scheme that works with tidal flats.

The options 4 and 14 can be defined along with the keyword definition CORRECTION ON CONVECTION VELOCITY : YES (logical, default NO) that modifies the depth-averaged convection velocity to account for the vertical gradients of velocity and concentration. This setting avoids overestimating suspended load, especially in deep waters, but it is not used in this tutorial.

The SCHEME OPTION FOR ADVECTION OF SUSPENDED SEDIMENTS can be additionally defined to either use a strong (default of 1) or a weak (2) form for advection. A weak form decreases numerical Diffusion, is more conservative, and increases computation time (read more in the Telemac2d steady section).

/ continued: gaia-morphodynamics.cas
/ ...
SCHEME FOR ADVECTION OF SUSPENDED SEDIMENTS : 14
/ CORRECTION ON CONVECTION VELOCITY : YES / use when SCHEME is 4 or 14 for deep water

Read more about the definition of numerical parameters in section 2.1.5 in the Gaia manual.

Morphological Coupling¶

When suspended load is activated together with bed evolution, the erosion and deposition fluxes contribute to the mass balance of the bed through the Exner equation. The net flux (erosion minus deposition) modifies the bed elevation at each time step.

Morphological factor for suspended load

For long-term simulations where morphological timescales are much longer than hydrodynamic timescales, a MORPHOLOGICAL FACTOR can be applied to accelerate bed evolution. This factor multiplies the net erosion/deposition flux, allowing multi-year morphological simulations with reasonable computation times. However, use with caution: morphological factors greater than 10-20 may introduce unrealistic results. The keyword is defined in the Gaia steering file:

MORPHOLOGICAL FACTOR : 10. / accelerate bed evolution 10x


## Example Applications

Examples for the implementation of suspended load come along with the TELEMAC installation (in the `/telemac/examples/gaia/` directory). The following examples in the `gaia/` folder feature (pure) suspended load calculations:

* 2d model of combined cohesive and non-cohesive suspended transport: **hippodrome-t2d/**
* 2d model of cohesive mud mass conservation: **mud_conservation-t2d/**
* 3d model of combined cohesive and non-cohesive suspended transport: **hippodrome-t3d/**
* 3d model of non-cohesive suspended transport with skin friction correction: **lyn-t3d/**
* 3d model of cohesive suspended transport with Rouse vertical profile (cf. [Gaia manual](https://gitlab.pam-retd.fr/otm/telemac-mascaret/-/raw/v9.0.0/documentation/gaia/user/gaia_user_9.0.pdf), section 2.1.2): **rouse-t3d/**
* 3d model of a tidal flume with cohesive sediment: **tidal_flats-t3d/**
* Coupling with waves: **sandpit-t2d/**

```{admonition} Recommended workflow for suspended load simulations
:class: note
1. **Start with hydrodynamics**: Ensure the hydrodynamic model (Telemac2d/3d) is calibrated and produces reasonable flow fields before coupling with Gaia.
2. **Define sediment classes**: Specify grain sizes appropriate for the site. Fine sediments ($D < 0.063$ mm) are typically cohesive; coarser sediments are non-cohesive.
3. **Select suspension formula**: Choose based on the environment (fluvial: `1` or `3`; coastal with waves: `4`).
4. **Set initial conditions**: Use measured or estimated suspended sediment concentrations.
5. **Choose advection scheme**: Use `14` for robustness with tidal flats, or `5` for better accuracy in deep channels.
6. **Calibrate erosion/deposition**: Adjust Partheniades constant $M$, critical shear stresses, and settling velocities to match observed concentrations.
7. **Validate mass balance**: Enable `MASS-BALANCE : YES` in the hydrodynamics steering file to monitor sediment conservation.

References¶

Rouse, H. (1939). An Analysis of Sediment Transportation in the Light of Fluid Turbulence. United States Department of Agriculture. https://resolver.caltech.edu/CaltechAUTHORS:20140529-132455484
Krone, R. B. (1962). Flume studies of the transport of sediment in estuarial shoaling processes final report [Phdthesis, Hydraulic Engineering Laboratory]. http://catalog.hathitrust.org/api/volumes/oclc/8967084.html
Stokes, G. G. (1850). On the effect of the internal friction of fluids on the motion of pendulums. In Transactions of the Cambridge Philosophical Society: Vol. IX (p. 99). http://mural.uv.es/daroig/documentos/stokes1850.pdf
Zanke, U. (1977). Berechnung Der Sinkgeschwindigkeiten von Sedimenten. Mitteilungen Des Franzius-Instituts Für Wasserbau, Technische Universität Hannover, 46, 243–261.
Dey, S. (2014). Fluvial Processes: Meandering and Braiding. In Fluvial Hydrodynamics (pp. 529–562). Springer-Verlag.
Richardson, J. F., & Zaki, W. N. (1954). Sedimentation and Fluidisation: Part I. Chemical Engineering Research and Design, 32, S82–S100. 10.1016/S0263-8762(97)80006-8
Partheniades, E. (1965). Erosion and Deposition of Cohesive Soils. Journal of the Hydraulics Division, 91(1), 105–139. https://cedb.asce.org/CEDBsearch/record.jsp?dockey=0013640
Celik, I., & Rodi, W. (1988). Modelling Suspended Sediment Transport in Nonequilibrium Situations. Journal of Hydraulic Engineering, 114(10), 1157–1191. 10.1061/(ASCE)0733-9429(1988)114:10(1157)
Zyserman, J. A., & Fredsøe, J. (1994). Data Analysis of Bed Concentration of Suspended Sediment. Journal of Hydraulic Engineering, 120(9), 1021–1042. 10.1061/(ASCE)0733-9429(1994)120:9(1021)
Guy, H. P., Simons, D. B., & Richardson, E. V. (1966). Summary of Alluvial Channel Data from Flume Experiments, 1956-61. US Government Printing Office. https://pubs.usgs.gov/pp/0462i/report.pdf
Bijker, E. W. (1992). Mechanics of Sediment Transport by the Combination of Waves and Current [Techreport]. ICCE. https://repository.tudelft.nl/islandora/object/uuid%253A4574aaa1-9246-4981-8305-c3a932e4ce03
Van Rijn, L. C. (1984). Sediment Transport, Part II: Suspended Load Transport. Journal of Hydraulic Division, 110(11), 1613–1641. 10.1061/(ASCE)0733-9429(1984)110:11(1613)
Soulsby, R. (1997). Dynamics of marine sands. Thomas Telford Publishing. 10.1680/doms.25844
Van Rijn, L. C. (2007). Unified View of Sediment Transport by Currents and Waves. I: Initiation of Motion, Bed Roughness, and Bed-Load Transport. Journal of Hydraulic Engineering, 133(6), 649–667. 10.1061/(ASCE)0733-9429(2007)133:6(649)