Table of Contents

SISYPHE User Manual links from User Manuals, SISYPHE


Introduction


Sisyphe is a sediment transport and morphodynamic simulation module which is part of the hydroinformatic finite elements system Telemac. In this module, sediment transport rates, decomposed into bed-load and suspended load, are calculated at each grid point as a function of various flow (velocity, water depth, wave height, etc.) and sediment (grain diameter, relative density, settling velocity, etc.) parameters. The bedload is calculated by using one of the classical sediment transport formulae from the literature. The suspended load is determined by solving an additional transport equation for the depth-averaged suspended sediment concentration. The bottom evolution equation (the Exner equation) can be solved by using either a finite element or a finite volume formulation. Sisyphe is applicable to non-cohesive sediments (uniform or graded), cohesive sediments as well as to sand-mud mixtures. The sediment composition is represented by a finite number of classes, each characterized by its mean diameter, grain density and settling velocity. Sediment transport processes can also include the effect of bottom slope, rigid beds, secondary currents, slope failure, etc. For cohesive sediments, the effect of bed consolidation can be accounted for. Sisyphe can be applied to a large variety of hydrodynamic conditions including rivers, estuaries and coastal applications, where the effects of waves superimposed to a tidal current can be included. The bed shear stress, decomposed into skin friction and form drag, can be calculated either by imposing a friction coefficient (Strickler, Nikuradse, Manning or Chézy) or by a bed-roughness predictor.

Coupling with hydrodynamics

In Sisyphe , the relevant hydrodynamic variables can be either imposed in the model or calculated by a hydrodynamic computation (“chaining method” or “internal coupling”). It is convenient to use one of the hydrodynamic modules of the Telemac system (Telemac-2d , Telemac-3d or Tomawac) for compatibility reasons (same mesh, same pre- and postprocessor, etc.) and because modules are already coupled, but the user can also choose a different hydrodynamic model. The different methods which can be used to prescribe the hydrodynamics are described next.

Numerical kernel

Sisyphe can be run on the following OS: Unix, Linux and Windows. The latest release of Sisyphe (6.2) uses the version 6.2 of the Bief finite element library (Telemac system library).


Running a sedimentological computation

Input Files

Steering File

The steering file contains the necessary information for running a computation, plus the values of parameters that are different from the defaults. Like in other modules of the Telemac system, the model parameters can be specified in the (obligatory) Sisyphe steering file. The following essential information (input/output) should be specified in the steering file:

Input data files

The geometry file is a binary file that contains the finite element mesh information (e.g., element type, number of elements, coordinates of the nodes and connectivity matrix). This file can also contain bottom topography information and/or friction values at each mesh point.

The geomety file (GEOMETRY FILE) is a binary file in Selafin format issued from the mesh generator. They contain a full description of the (three-node) triangular grid as well as the initial sediment bed level. This initial configuration can also be obtained from the interpolation of bathymetric data file onto each node of the triangular grid.

The FORTRAN FILE contains those subroutines which the user has modified in order to adapt them to their specific needs.

The HYDRODYNAMIC FILE contains the hydrodynamic results of a previous hydrodynamic computation (Telemac-2d or Telemac-3d ). If Telemac is coupled to Sisyphe , the grid must be the same in both cases. When a results file from a finite differences hydrodynamic model needs to be read, the hydrodynamic results file shall be transformed beforehand and interpolated onto the triangular mesh.

The hydrodynamics will be given either by the last time step of that file in steady mode, or by the time steps related to the period of the tide or flood being considered in unsteady mode. For an explanation on steady/unsteady mode, see §XX). That file should contain current data and wave data if the effect of waves is considered. The WAVE FILE contains the hydrodynamic results of a previous wave computation (Tomawac ). Only the last record will be read if the keyword WAVE EFFECTS is activated.

The PREVIOUS SEDIMENTOLOGICAL COMPUTATION FILE is provided for implementing a computational sequence, it contains the results of a previous sedimentological computation made using Sisyphe on the same grid. If that file exists, then the initial sedimentological conditions of the computation will be given by the last time step in the file. This file will be read only if the keyword COMPUTATION CONTINUED is activated.

The REFERENCE FILE is the results file from a previous computation. The last time step is compared to the new last time step. All the files must be either binary (Selafin format) or text (ascii format).

Geometry file

A binary file containing information of the triangular grid. It can be used as input file for all hydrodynamic and wave propagation models in the Telemac Modelling System (Telemac-2d , Telemac-3d or Tomawac ).

Boundary conditions file

The format of the boundary condition file is the same as for Telemac-2d or Telemac-3d. This file can be created or modified using a text editor. Each line in this file is related to a grid boundary. Boundary points are listed in the file in the following way: First the domain outline points, proceeding counterclockwise from the lower left corner, then the islands proceeding clockwise , also from the lower left corner. The boundary conditions for sediment transport computations are imposed at the specific place of the tracer boundary condition as in Telemac-2d .

Obligatory keywords
Names of principal input/output files (other files are optional)

When combining a previous computation:


Coupling hydrodynamics/morphodynamics

We describe here two possibilities for linking the hydrodynamic and the morphodynamic models: chaining (the flow is obtained from a previous hydrodynamic simulation assuming a fixed bed) or internal coupling (both the flow and bed evolution are updated at each time step).

Chaining method

Both models (hydrodynamic and morphodynamic) are run independently: during the first hydrodynamic simulation the bed is assumed to be fixed. Then, in the subsequent morphodynamic run, the flow rate and free surface are read from the previous hydrodynamic results file. This ‘chaining method’ is only justified for relatively simple flows, due to the difference in time-scales between the hydrodynamics and the bed evolution. For unsteady tidal flow, Sisyphe can be used in an unsteady mode: the flow field is linearly interpolated between two time steps of the hydrodynamic file. For steady flow, the last time step of the hydrodynamic file is used and the flow rate and free surface assumed to stay constant while the bed evolves.

At each time step, the flow velocity is updated by assuming simply that both the flow rate and the free surface elevation are conserved, such that, in the case of deposition, the flow velocity is locally increased, whereas in the case of erosion, the flow velocity decreases. This rather schematic updating does not take into account any deviation of the flow. It is only suitable for simple flows (2D processes) and assuming relatively small bed evolutions. However it can be responsible for numerical instabilities (cf. [01b], Hervouet and Villaret, 2004). The morphodynamic computation is stopped when the bed evolution reaches a certain percent of the initial water depth. This simple updating of the flow field is no longer valid when the bed evolution becomes greater than a significant percentage of the water depth, specified by the user. At this point, it is recommended to stop the morphodynamic calculation and to recalculate the hydrodynamic variables.

It should be noted that with this simple method, the sediment mass continuity may not be satisfied because of potential losses due to changes in the flow depth as the bed evolves. When the flow is steady (‘steady case’ = YES), only the last record of the previous result file will be used. Otherwise (‘steady case’ = NO), the ‘Tide Period’ and ‘Number of tides or floods’ will be used to specify the sequence to be read on the hydrodynamic files. Hydrodynamic records are interpolated at each time step of the sedimentological computation. Note: an error may occur when the ‘Tidal Period’ is not a multiple of the graphical time steps of the hydrodynamic file (‘Hydrodynamic file is not long enough’). In an unsteady case, the keyword ‘Starting time of the hydrogram’ gives the first time step to be read. If the starting time is not specified, the last period of the hydrogram will be used for sedimentological computation.

For uncoupled mode, the Sisyphe steering file should specify:

- The time steps, graphical or listing output, duration
- The hydrodynamic file as yielded by Telemac -2D, -3D (‘hydrodynamic file’) or even by user-subroutine (condim_sisyphe.f).
- For waves only: the wave parameters can be either calculated by a wave propagation code (TOMAWAC), or defined directly in Sisyphe (condim_sisyphe.f). The effect of waves on bed forms and associated bed roughness coefficient can be accounted with keyword: Effect of waves =’YES’.
- A restart from a previous Sisyphe model run, by setting ‘Computation Continued = Yes’ and specification of sedimentological results in ‘Previous sedimentological computation ’
- Flow options: steady or unsteady options, flow period

Keywords
For time step, duration and output:

For hydrodynamics (imposed flow and updated):

For waves:

Internal coupling

For a more realistic simulation, Sisyphe can be automatically coupled (internally) with the hydrodynamic model, Telemac-2d or – 3D. Sisyphe is called inside the hydrodynamic model without any exchange of data files. The data to be exchanged between the two programs is now directly shared in the memory, instead of being written and read in a file. At each time step, the hydrodynamics variables (velocity field, water depth, bed shear stress,…) are transferred to the morphodynamic model, which sends back the updated bed elevation to the hydrodynamics model.

The internal coupling method is more CPU time consuming than the chaining method. Various techniques can be set up to reduce the CPU time (e.g. parallel processors, see §9). In certain cases, the use of a coupling period >1 allows the bed load transport rates and resulting bed evolution not to be re-calculated at every time step. For the suspended load, a diffusion/transport equation needs to be solved. This transport equation obeys the same Courant number criteria on the time step than the hydrodynamics, and therefore needs to be solved at each time-step (‘coupling period = 1’).

The time step of Sisyphe is equal to the time step of Telemac-2d or – 3D multiplied by the ‘coupling period’. The graphic/listing printout periods are the same as in the Telemac computation.

The Telemac-2d /3D steering file must specify the type of coupling, the name of the Sisyphe steering file, and the coupling period. In addition, the Fortran file of Sisyphe must be sometimes specified in the Telemac steering file (if there is no Fortran file for Telemac ). Some of the keywords of the Sisyphe steering file become obsolete.

Keywords
For internal coupling, the following keywords need to be specified in the Telemac-2d (or –3D) steering files:

All computational parameters (time step, duration, printout, option for friction) need to be specified in the Telemac steering file, but are no longer used by Sisyphe . The values of time step, bottom shear stress… are transferred directly from Telemac to Sisyphe .

keywords no longer in use is (Sisyphe ) steering file:


Flow-sediment interactions


Sediment and fluid properties

Fine sediment particles of diameters less than $60\mu m$ present complex cohesive properties which affect the sediment transport processes. For non-cohesive sediments (median diameter $d_{50} >60 \mum$), the grain diameter and grain density $\rho_s$ are the key parameters which determine its resistance to erosion and sediment transport rate.

For cohesive sediments ($d_{50} < 60 \mum$), the grain diameter is no longer the key sediment parameter: The settling velocity now depends on the concentration and other physico-chemical properties of the suspension due to flocculation, whereas the critical bed shear strength depends on the consolidation state of the sediment bed.

For the simplest case, we consider uniform, non-cohesive sediments characterized by one single value for the grain size $d_{50}$ and grain density $\rho_s=$ with $\rho_s=$2650 kg\,m$^{-3}$ which can be transported both as bed-load and suspended load.

Keywords
The physical properties of the sediment are always defined in the Sisyphe steering file using the following keywords:

The physical properties of the fluid are defined by:

Bed shear stress

Hydrodynamic model

The current-generated bed shear stress is used in both the shallow water momentum equation as well as the bottom boundary condition for the velocity profile. When Sisyphe is coupled with Telemac-2d , the bed shear stress t0 is calculated at each time step from the depth-averaged mean flow velocity by use of a quadratic friction coefficient Cd as follows:

\begin{equation*} \tau _0 =\frac{1}{2} \rho C_d U^{2} \end{equation*}

where $\rho$ is the fluid density, and $U$ the depth-averaged velocity. When \sisyphe is coupled with \telddd, the bed shear stress is aligned with the near bed velocity in order to account for possible vertical flow deviations. The magnitude of the bed shear stress is still related to the depth-averaged velocity, except if a Nikuradse friction law is applied. In that case, the friction velocity u$_{*}$, defined by

\begin{equation*} \tau _0 = \rho u_*^2 \end{equation*}

is related to the near bed flow velocity $u(z_1)$) by a logarithmic velocity profile:

\begin{equation*} u(z_1) = \frac{u_*}{\kappa} \ln\left(\frac{z_1}{z_0}\right) \end{equation*}

where $\kappa = 0.4$ is the von K\'arm\'an constant, $z_0$ is the vertical distance from a rough boundary, expressed as a function of the Nikuradse bed roughness ($z_0 = k_s/30$), with $k_s$ the grain roughness height. For flat beds, the roughness height has been shown to be approximately $k_s\approx 3 d_{50}$

Role of bed forms

A natural sediment bed is generally covered with bed forms (length $\lambda$ and height $\eta$). The presence of bed forms greatly modifies the boundary layer flow structure, with the formation of recirculation cells and depressions in the lee of bed forms. Depending on the flow rate and the sediment transport rate, the size of bed forms ranges from a few centimeters for ripples to a few tens of meter for mega-ripples. The dimensions of dunes scale with the water depth ($\eta_D \approx 0.4 h$, $\lambda_D\approx 6-10h$). In most cases, large scale models do not resolve the small to medium scale bed forms (ripples, mega-ripples) which need therefore to be parameterized by an increased friction coefficient:

\begin{equation*} \tau_0 = \tau' + \tau'' \end{equation*}

The local skin friction component determines the bed-load transport rate and the equilibrium concentration for the suspension. The total friction velocity determines the (spatially averaged along bed forms) turbulence eddy viscosity/diffusivity vertical distribution in 3D models, and therefore determines both velocity vertical profile and mean concentration profile.

Sediment transport

Bedload and suspended load

When the current-induced bed shear stress increases above a critical threshold value, sediment particles start to move as bed load, while the finer particles are transported in suspension. Bedload occurs in a very thin high concentrated near-bed layer, where inter-particle interactions occur. The suspended load is defined as the depth-integrated flux of sediment concentration, from the top of the bedload layer up to the free surface. For dilute suspension ,clear flow concepts (turbulence diffusion, eddy viscosity, logarithmic velocity profile) are considered to be valid.The total sediment load $Q_t$ includes both a bedload $Q_b$ and suspended load $Q_s$:

\begin{equation*} Q_t = Q_b + Q_s. \end{equation*}

Shields parameter

Most sand transport formulae assume a critical value of the bed shear stress, noted $\tau_c$, for the onset of erosion: $Q_s = 0$ for $\tau_0 < \tau_c$. The Shields parameter or non-dimensional critical shear stress $\theta_c$ is defined by:

\begin{equation*} \theta_c = \frac{\tau_c}{g(\rho_s -\rho)d_{50}} \end{equation*}

can be either specified in the \sisyphe steering file by use of keyword SHIELDS PARAMETERS, or calculated by the model as a function of non-dimensional grain diameter. The Van Rijn formula has been programmed

\begin{equation*} \begin{align*} D_* \leq 4,\quad & \theta_c = 0.24 D_*^{-1} \\ 4 < D_* \leq 10,\quad & \theta_c = 0.14 D_*^{-0.64} \\ 10 < D_* \leq 20,\quad & \theta_c = 0.04 D_*^{-0.10} \\ 20 < D_* \leq 150,\quad & \theta_c = 0.013 D_*^{0.29} \\ 150 \leq D_*,\quad & \theta_c = 0.045, \end{align*} \end{equation*}

where the non-dimensional diameter $D_*$ is defined by:

\begin{equation*} D_* =\left(\frac{g(\rho_s/\rho-1)}{\nu^2}\right) ^{1/3} d_{50}. \end{equation*}

Keywords

Friction laws

Differences between 2D and 3D models

The intensity of bed shear stress is calculated by use of same friction laws in both 2D and 3D models: the bed shear stress is related to the depthaveraged flow velocity, except when the Nikuradse friction law is applied. In 3D-models, the bed shear stress is directly related to near bed velocity assuming a log profile (eq. 2), whereas in 2D-models, the log profile is assumed to be valid up to the free surface.

The direction of the bed shear stress and resulting bedload transport rate is assumed to be in the direction of the depth-averaged velocity in Sisyphe (alone or internally coupled with Telemac-2d ). When Sisyphe is internally coupled to Telemac-3d , the bed shear stress and resulting transport rate are assumed to be in the direction of the near bed velocity. The 3D model gives a more accurate estimate of the bottom friction, since it accounts for a possible vertical deviation of the current.

Coupling with TELEMAC-2D

The depth-averaged bed shear stress and resulting bedload transport rates are assumed to be in the direction of the mean flow velocity, except when the sediment transport formulation accounts for:

Coupling with TELEMAC-2D
When the model is coupled with Telemac-2d , the values of the friction coefficients (and therefore, the bed shear stress) are provided by Telemac-2d .

Coupling with TELEMAC-3D

Both bed shear stress and bed load transport rate are aligned with the near bed velocity. This 3D approach is more physical, and takes into account possible recirculation and veering of the flow, such that corrections for secondary currents for example are no longer necessary.

Uncoupled model

The quadratic friction coefficient Cd which is used to calculate the total bed shear stress can be calculated based on the selected friction law. Different options, which are consistent with the Telemac-2d options, are available in Sisyphe and depend on the choice of the ‘LAW FOR BOTTOM FRICTION’ and on the value of the ‘FRICTION COEFFICIENT’. Keywords:

\begin{equation*} C_d = \frac{2g}{C_h^2} \end{equation*}

\begin{equation*} C_d = \frac{2g}{S_t^2}\frac{1}{h^{1/3}} \end{equation*}

\begin{equation*} C_d =\frac{2g}{h^{1/3} } M_a^2 \end{equation*}

\begin{equation*} C_d = 2\left[\frac{\kappa}{\log(\frac{12h}{k_s})}\right]^2 \end{equation*}

Keywords
In the Sisyphe steering file, the total bed shear stress is calculated based on:

Similar keywords are available in both Telemac-2d and -3D in the case of internal coupling.

Skin friction correction

Total bed shear stress decomposition

Bedload transport rates are calculated as a function of the local skin friction component t0. The total bed shear stress issued from the hydrodynamics model needs to be corrected in the morphodynamics model as follows:

\begin{equation*} \tau'= \mu \tau_0, \end{equation*}

Physically, the skin bed roughness should be smaller than the total bed roughness (i.e. $\mu \leq 1$). However, in most cases the hydrodynamic friction does not represent the physical bottom friction: the coefficient is generally used as a calibration coefficient in hydrodynamics models. It is adjusted by comparing simulation results with observations of the time-varying free surface and velocity field. Therefore, its model value integrates various neglected processes (side wall friction, possible errors in the bathymetry and input data). Under those conditions, a correction factor $\mu > 1$ can be admitted. Different methods are programmed in \sisyphe in order to calculate the bedform correction factor $\mu$, according to the keyword `Skin friction correction'(ICR) :

\begin{equation*} \mu =\frac{C_d'}{C_d}, \end{equation*}

where $C_d$ et $C_d'$ are both quadratic friction coefficients related respectively to total friction and skin friction: Cd is obtained from Telemac-2d or Telemac-3d and $C_d'$ is calculated from $k_s'$:

\begin{equation*} C_d' = 2\left[\frac{\kappa}{\log(\frac{12h}{k_s'})}\right]^2 \end{equation*}

\begin{equation*} \mu =\frac{C_d'^{0.75} C_r^{0.25}}{C_d}, \end{equation*}

where the quadratic friction $C_r$ due to bedforms is calculated as a function of $k_r$.

Keywords The option selected for the skin friction correction is based on keywords:

Bed roughness predictor

Total bed roughness (uncoupled mode)

Different options are programmed in Sisyphe to predict the total bed roughness. It is recalled that the bed friction option of Sisyphe is not used in the case of internal coupling with Telemac .

For waves and combined waves and currents, bedform dimensions are calculated as a function of wave parameters following the method of Wiberg and Harris. The wave-induced bedform bed roughness $k_r$ is calculated as a function of the wave-induced bedform height $\eta_r$:

\begin{equation*} k_r = \max(k_s', \eta_r). \end{equation*}

For currents only, the van Rijn predictor has been implemented. The total bed roughness can be decomposed into a grain roughness $k_s'$, a small-scale ripple roughness $k_r$, a mega-ripple component $k_{mr}$, and a dune roughness $k_d$:

\begin{equation*} k_s = k_s' + \sqrt{k_r^2 + k_{mr}^2 + k_d^2}. \end{equation*}

Both small scale ripples and grain roughness have an influence on the sediment transport laws, while the mega-ripples and dune roughness only contribute to the hydrodynamic model (total friction).

Bedform correction factor

In order to account for the effect of ripples on the skin friction, the keyword ‘Skin friction correction’ should be set to ICR = 2. For currents only, the ripple bed roughness is a function of the mobility number, according to Van Rijn :

\begin{equation*} \begin{align*} k_r &= d_{50}(85-65\tanh(0.015(\Psi-150))),\,\,\text{for\,\,} \Psi<250\\ k_r &= 20 d_{50},\,\,\text{for\,\,} \Psi \geq 250, \end{align*} \end{equation*}

where $\Psi =U^2/(s-1)gd_{50}$

Keywords


Bed-load transport


Exner equation

Equilibrium conditions

Bed-load occurs in a thin near-bed high-concentrated region. The bedload layer therefore adapts very rapidly to any changes in the flow conditions, such that equilibrium conditions can be considered to be valid. The bed-load transport rate can then be calculated by use of some equilibrium sediment transport formula, as a function of various flow and sediment parameters, assuming that the transport rate corresponds to saturation conditions.

Bed evolution

To calculate the bed evolution, Sisyphe solves the Exner equation:

\begin{equation*} (1-n)\frac{\partial Z_f}{\partial t} + \grad\cdot Q_b =0, \end{equation*}

where $n$ is the bed porosity ($n \approx 0.4$ for non cohesive sediment), $Z_f$ the bottom elevation, and $Q_b$ (m$^{2}$s$^{-1}$) the solid volume transport (bed-load) per unit width.

The previous Equation states that the variation of sediment bed thickness can be derived from a simple mass balance. It is only valid for equilibrium conditions. The previous Equation is strictly valid for bed-load only. However, it can be extended to total load (including the suspended load) assuming quasi-steady and uniform flow conditions. Otherwise, a different treatment should be applied to both bed load and suspended load, in order to take non-equilibrium flow conditions into account properly.

Keywords

Bedload transport formulae

For currents only (no wave effects), a large number of semi-empirical formulae can be found in the literature to calculate the bedload transport rate. Sisyphe offers the choice between 10 different bedload formulae including the Meyer-Peter and Müller, Engelund-Hansen and Einstein- Brown formulae. Most sediment transport formulae assume threshold conditions for the onset of erosion (e.g. Meyer-Peter andMüller, van Rijn and Hunziker). Other formulae are based on similar energy concept (e.g. Engelund- Hansen) or can be derived from a statistical approach (e.g. Einstein- Brown, Bijker, etc.). The non-dimensional current-induced sand transport rate, noted $\Phi_s$, can be expressed as a function of the non-dimensional skin friction or Shields parameter, noted $\theta'$ , as defined by:

\begin{equation*} \theta'=\frac{\mu\tau_0}{(\rho_s-\rho)gd_{50}}, \end{equation*}

$\mu$ , the correction factor for skin friction

\begin{equation*} \Phi_s = \frac{Q_s}{\sqrt{g(s-1)d_{50}^3}} \end{equation*}

Keywords

Sediment transport formulae

Choice of formulae

The choice of a transport formula depends on the selected value of the model parameter ICF , as defined in the steering file by the keyword ‘bedload transport formula’ (ICF). The user can program a specific transport formula through the Qsform.f(user-subroutine), by setting ICF=0. Other sediment transport formulae, described in chapters 8 and 9, have been programmed in Sisyphe to account for the effects of waves (cf. Bijker, Bailard, Dibajnia and Watanabe, etc.) or sand grading (cf. Hunziker). They can nevertheless be used for currents only and uniform grain size.

\begin{equation*} \begin{align*} \Phi_b &=0\quad\text{if}\,\theta'<\theta_c \\ \Phi_b &=8(\theta'-\theta_c)^{3/2}\quad\text{otherwise} \end{align*} \end{equation*}

where $\theta_c$ is the critical Shields parameter ($\theta_c = 0.047$ by default).

\begin{equation*} \Phi_b = F(D_*)f(\theta'), \end{equation*}

with

\begin{equation*} F(D_*) = \left(\frac{2}{3} +\frac{36}{D_*}\right)^{0.5} - \left( \frac{36}{D_*}\right)^{0.5}, \end{equation*}

and

\begin{equation*} \begin{align*} f(\theta') = 2.15\exp(-0.391/\theta') & \quad\text{if}\,\,\theta' \leq 0.2 \\ f(\theta') = 40\,\theta^{'3} & \quad\text{if}\,\,\theta' >0.2 \end{align*} \end{equation*}

ICF = 30 corresponds to the original formula, where the transport rate is related to the skin friction without threshold:

\begin{equation*} \Phi_s = 0.1\,\theta'^{5/2} \end{equation*}

ICF = 3 corresponds to the version modified by Chollet and Cunge in order to account for the effects of sand dunes formation. The transport rate is here related to the total bed shear stress:

\begin{equation*} \Phi_s = \frac{0.1}{C_d} \hat{\theta}^{5/2}, \end{equation*}

the dimensionless bed shear stress $\theta$ is calculated as a function of the dimensionless skin friction $\theta'$:

\begin{equation*} \begin{align*} \hat{\theta} = 0 &\quad \text{if}\,\theta' < 0.06\,\text{(flat bed regime - no transport)} \\ \hat{\theta} = \sqrt{2.5(\theta'-0.06)} &\quad \text{if}\, 0.06 < \theta' < 0.384\,\text{(dune regime)} \\ \hat{\theta} = 1.065\theta'^{0.176} &\quad \text{if}\, 0.384 < \theta' < 1.08\,\text{(transition regime)} \\ \hat{\theta} = \theta' &\quad \text{if}\,\theta' > 1.08\,\text{(flat bed regime - upper regime)} \end{align*} \end{equation*}

Notice
In the case of coupling between bed-load and suspended load (SUSPENSION = YES), a total load formula (N$^\circ$3 and 30) should not be used, in order to avoid the suspended load to be calculated twice.

\begin{equation*} \Phi_b = 0.053 D_*^{-0.3} \left( \frac{\theta_p-\theta_{cr}}{\theta_{cr}} \right)^{2.1}. \end{equation*}

Validity range
Most sediment transport formulae are based on experiments performed under steady-flow conditions. They demonstrate rapid variation of the bed-load transport rates as a function of the mean flow intensity: an increase of the current velocity by $10\%$ will result, depending on the formula being used, in an increase of the transport rate by over $30\%$ (Meyer-Peter), $60\%$ (Engelund-Hansen) or almost $80\%$ (Einstein-Brown). Therefore, any error made when calculating the hydrodynamics will be significantly amplified by the sediment transport rates estimates. On the other hand, under variable flow conditions (e.g. in tidal regime), the average transport will be highly influenced by the stronger currents and will not be directly related to the mean flow.

The validity range of the different formulae can be summarized in the table below:

Formula Meyer-Peter Einstein-Brown Engelund-Hansen van Rijn
IFC 1 2 3 or 30 7
Mode of transport Bed-load Bed-load Total-load Bed-load
Sand diameter > 1 mm 0.2-3 mm 0.2-4 mm 0.2-2 mm
Rippled bottoms Yes Yes Yes XXX
Low bed-load flows Yes Yes No XXX
High bed-load flows Yes Yes Yes XXX

Sloping bed effect

The effect of a sloping bottom is to increase the bedload transport rate in the downslope direction, and to reduce it in the upslope direction. In Sisyphe , we apply a correction factor to both the magnitude and direction of the solid transport rate, before solving the bed-evolution equation. This sloping bed effect is only activated if the keyword ‘SLOPE EFFECT’ is activated. Two different formulations for both effects are available depending on the choice of ‘FORMULA FOR SLOPE EFFECT’ (SLOPEFF), which chooses the magnitude correction and ‘FORMULA FOR DEVIATION’ (DEVIA), which chooses the direction correction.

Magnitude of bedload transport rate

\begin{equation*} 1-\beta\frac{\partial Z_f}{\partial s}, \end{equation*}

where $s$ is the coordinate in the current direction and $\beta $ is an empirical factor, which can be specified by keyword 'BETA' (default value = 1.3). This effect of bed slope is then similar to adding a diffusion term in the bed-evolution equation. It tends to smoothe the results and reduces unstabilities.

\begin{equation*} \frac{\theta_c}{\theta_{co}} = \frac{\cos\psi \sin\beta + ( \cos^2\beta \tan^2\phi_s - \sin^2\psi \sin^2\beta)^{0.5}}{\tan \phi_s}, \end{equation*}

where $\beta$ is the angle of the bed to the horizontal plane and $\psi$ the angle of the current to the upslope direction.

Direction of bedload transport rate

The change in the direction of solid transport is taken into account by the formula:

\begin{equation*} \tan\alpha = \tan\delta -T\frac{\partial Z_f}{\partial n}, \end{equation*}

where $\alpha$ is the direction of solid transport in relation to the flow direction, $\delta$ is the direction of bottom stress in relation to the flow direction, and $n$ is the coordinate along the axis perpendicular to the flow.

- DEVIA = 1: According to Koch and Flochstra, the coefficient T now depends on the Shields parameter:

\begin{equation*} T=\frac{4}{6\theta} \end{equation*}

- DEVIA = 2, the deviation is calculated based on Talmon et al. (1995) and depends on the empirical coefficient $\beta_2 = 0.85$.

\begin{equation*} T=\frac{1}{\beta_2\sqrt{\theta}} \end{equation*}

This default value can be modified by changing the keyword: `PARAMETER FOR DEVIATION' = $\beta_2$.

Keywords

Sediment slide

An iterative algorithm prevents the bed slope from becoming greater than the maximum friction angle ($\phi_s\approx 32^{\circ} - 40^{\circ}$). This option is activated by use of key word `SEDIMENT SLIDE'.

Further developments are needed to introduce the bed stability, in order to represent mud slides.

\begin{equation*} \left\| \overrightarrow{\grad(Z_f)} \right\| <\tan \phi_s \end{equation*}

Keywords

Secondary currents

In meandering streams, centrifugal inertial forces are induced by the curvature $R^{-1}$ of the streamlines. The depth-averaged mean pressure gradient due to the cross-sectional variation of the free surface balances the depth -averaged centrifugal force, noted F:

\begin{equation*} \begin{align*} F_y &= -\frac{1}{\rho}\frac{\partial P}{\partial y} \\ F_y &= \rho \frac{U^2}{R}-\frac{1}{\rho} \frac{\partial P}{\partial y} =-g\frac{\partial Z_s}{\partial y}, \end{align*} \end{equation*}

where $y$ is the cross-sectional axis and $\partial Z_s/\partial y$ the cross-sectional variation in the free surface. The mass balance in the y-direction gives a zero-depth averaged cross-sectional velocity.

\begin{equation*} \rho \frac{U^2}{R} =-g\frac{\partial Z_s}{\partial y} \end{equation*}

However, in the vertical, the centrifugal force increases with distance from the bed whereas the pressure gradient is constant. P and F are therefore out of balance (P greater in the lower part of the flow, whereas F is greater in the upper part). This results in a cross-sectional circulation, with a cross-velocity directed near the inner bank, and flow velocity therefore follows a spiral motion which cannot be represented by a depth-averaged hydrodynamic model.

The correction method, based on Engelund (? ), has been programmed in Sisyphe , in order to account for the effect of these radial stresses on the sand transport rates. This method is based on the following assumptions:

The near bed angular deviation of the main flow $\delta$ can be shown to be proportional to $h/R$, as explained in Yalin and Ferreira da Silva (2001). We use here $\tan\delta = 7h/R$. The bend radius $R$ is generally unknown, but can be inferred from the cross sectional variation in the free surface:

\begin{equation*} R=-\rho\alpha'\frac{U^2}{g\frac{\partial Z_s}{\partial y}}, \end{equation*}

\begin{equation*} R=-\rho\alpha'\frac{U^2}{g\frac{\partial Z_s}{\partial y}}, \end{equation*}

Keywords

Notice
This option should only be activated when the flow is calculated by a depth-averaged model.

Treatment of rigid beds

Different methods

Non-erodible beds are managed numerically by limiting bedload erosion and letting incoming sediment pass over. The problem of rigid beds is conceptually trivial but numerically complex. In finite elements the minimum water depth algorithm allows a natural treatment of rigid beds (cf. Hervouet et al. 2011). The sediment is managed as a layer with a depth that must remain positive, and the Exner equation is managed similarly to the shallow-water continuity equation.

With finite volumes, two methods are implemented in Sisyphe , depending on the ‘OPTION FOR THE TREATMENT OF NON ERODABLE BEDS’ :

Keywords

Subroutine noerod.f The position of the rigid bed is set by default to a very large negative value and can be changed in subroutine noerod.f.

Effect of tidal flats

Tidal flats are areas of simulation domain where depth can become 0 during simulation can be managed by two approaches:

Keyword


Suspended load


Suspended load transport equation

Passive scalar hypothesis

For non-equilibrium flow conditions, the bedload and the suspended load are dealt with separately treatment. The interface between the bed-load and suspended load is located at $z=Z_{ref}$:

Settling velocity

The settling velocity Ws is an important parameter for the suspension. It can be either specified, or calculated by the model as a function of grain diameter, by use of a semi-empirical formula. The Van Rijn formula which is valid for non-cohesive spherical particles and dilute suspensions, has been programmed in Sisyphe (cf. Van Rijn, 1993, cf. [31], p. 3.13):

\begin{equation*} \begin{align*} W_{s} &= \frac{(s-1)g d_{50}^2}{18\nu}, \quad \text{if } d_{50} \leq 10^{-4} \\ W_{s} &= \frac{10\nu}{d_{50}} \left(\sqrt{1+0.01\frac{(s-1)gd_{50}^3}{18\nu^2}}-1\right), \quad \text{if } 10^{-4} \leq d_{50} \leq 10^{-3}\\ W_{s} &= 1.1 \sqrt{(s-1)gd_{50}}, \quad \text{if } 10^{-3} \leq d_{50}, \end{align*} \end{equation*}

where $s=\rho_{s}/\rho_0$ is the relative density and $g$ is gravity.

Outside the bed-load layer sediment particles are assumed to follow the mean flow velocity, noted $\overrightarrow U$, with an additional vertical settling velocity, noted $W_s$.

3D transport equation

The 3D advection-diffusion equation for the suspended sediment concentration is:

\begin{equation*} \frac{\partial C}{\partial t} + \grad\cdot \left(( \overrightarrow{U} +\overrightarrow{W_s}) C\right) = \grad\cdot (\epsilon_s \grad C), \end{equation*}

where $\epsilon_s$ is the sediment diffusivity coefficient. The following boundary conditions must be satisfied both at the free surface ($Z_s$) and at the interface between the bed-load and the suspended load ($z=Z_{ref}$):

\begin{equation*} \begin{align*} \text{at } z &= Z_s,\quad \epsilon_s \frac{\partial C}{\partial z} + W_s C = 0 \\ \text{at } z &= Z_{ref},\quad -\epsilon_s \frac{\partial C}{\partial z} - W_sC = (E-D)_{Z_{ref}} \end{align*} \end{equation*}

where $Z_s$ is the free surface elevation $Z_{ref}$ the position of the interface between the bed-load layer and the suspension. $E-D$ is the net erosion minus deposition flux of sediment. The 3D transport equation and boundary conditions can be solved in the 3D sediment transport model of Telemac-3d.

2D transport equation

The 2D transport equation for the depth-averaged mean concentration C is obtained by depth-integration :

\begin{equation*} \frac{\partial hC}{\partial t} + \frac{\partial h\overline{uC}}{\partial x} + \frac{\partial h\overline{vC}}{\partial y} = \frac{\partial }{\partial x} \left(h\epsilon_s\frac{\partial C}{\partial x}\right) + \frac{\partial }{\partial y} \left(h\epsilon_s\frac{\partial C}{\partial y}\right) + E - D, \end{equation*}

with:

\begin{equation*} X = \overline{x} =\frac{1}{h} \int_{Z_{ref}}^{Z_s} x(z) \d z, \end{equation*}

where $h=Z_s-Z_f \approx Z_s-Z_{zref}$ is the water depth, assuming the bed-load layer thickness to be small. After simplification of the advection terms and using the continuity equation, the following approximate depth-averaged transport equation can be solved in its non-conservative form:

\begin{equation*} \frac{\partial C}{\partial t} + U_{conv} \frac{\partial C}{\partial x} + V_{conv} \frac{\partial C}{\partial y} = \frac{1}{h} \left(\frac{\partial }{\partial x} \left( h\epsilon_s \frac{\partial C}{\partial x} \right) + \frac{\partial }{\partial y} \left(h\epsilon_s \frac{\partial C}{\partial y} \right) \right) + \frac{(E-D)_{z=Zref}}{h}, \end{equation*}

where $U_{conv}$ and $V_{conv}$ are the depth-averaged convective flow velocities in the x and y directions.

\begin{equation*}\ \frac{\partial C}{\partial t} + U_{conv} \frac{\partial C}{\partial x} + V_{conv} \frac{\partial C}{\partial y} = \left(\frac{\partial }{\partial x} \left( \epsilon_s \frac{\partial C}{\partial x} \right) + \frac{\partial }{\partial y} \left(\epsilon_s \frac{\partial C}{\partial y} \right) \right) + \frac{(E-D)_{z=Zref}}{h}, \end{equation*}

The value of the dispersion coefficient $\epsilon_s$ depends on the choice of the parameter OPTDIF and different possibilities are offered (constant coefficient, by default). The Elder model (OPTDIF=2) allows to calculate both the longitudinal and transversal components as a function of the friction velocity. ($\epsilon_s = 6 u_*h$ and $\epsilon_s = 0.3 u_*h$).

Different schemes are available for solving the non-linear advection terms depending on the choice of the parameter the classical characteristics schemes (to the diffusive schemes SUPG and PSI, it is recommended to use the conservative schemes

Keywords

Tidal flats
It is recommended in the presence of tidal flats to use the conservative scheme based on the calculation of flux per segments. Scheme 14 should be used if the convection velocity differs from the depth-averaged velocity and no longer satisfies the shallow-water continuity equation (‘correction on convection velocity’ =YES ).

Bed evolution

Bed evolution due to the suspension

The bottom-evolution equation due to the suspension accounts for the vertical fluxes at the interface between the bed-load and the suspended load. The bed-evolution equation becomes:

\begin{equation*} (1-n)\frac{\partial Z_f}{\partial t} + (E-D) = 0 \end{equation*}

where $Z_f$ is the bottom elevation, $Z_{ref}$ the interface between the bed-load and suspended-load, and $n$ the bed porosity. Erosion and deposition fluxes are calculated at the bed-load/suspended load interface ($z = Z_{ref}$).

Units consistency
The net sediment flux E-D is in m s$^{-1}$. For consistency with the bed-load units (m$^{2}$s$^{-1}$), concentrations and bed porosity are dimensionless. However, the user can choose mass concentration for graphic printouts, by use of key word `MASS CONCENTRATION'. The relation between volume concentration and mass concentration is : $C_m$ (kgm$^{-3}$) = $\rho_s C$.

Keywords

Equilibrium concentrations

Erosion/Deposition rate

For non-cohesive sediments, the net erosion minus deposition flux is determined based on the concept of equilibrium concentration (Celik and Rodi, 1988):

\begin{equation*} \left(E-D \right)_{Z_{ref}} = W_s \left(C_{eq} - C_{Z_{ref}}\right) \end{equation*}

where $C_{eq}$ is the equilibrium near-bed concentration and $C_{zref}$ is the near-bed concentration, calculated at the interface between the bed-load and the suspended load, $z=Z_{ref}$. Two different semi-empirical formula from the literature are programmed in Sisyphe :

\begin{equation*} C_{eq} =\frac{0.331(\theta'-\theta_c)^{1.75}}{1+0.72(\theta'-\theta_c)^{1.75}}, \end{equation*}

where $\theta_c$ is the critical Shields parameter and $\theta'= \mu\theta$ ,the non-dimensional skin friction which is related to the Shields parameter, by use of the ripple correction factor.

\begin{equation*} C_{eq} = \frac{Q_sB}{b Z_{ref} u_*} \end{equation*}

with an empirical factor b = 6.34.

For ICQ= 2…
For ICQ= 2, both bed-load and suspended load should be calculated at each time step. The ‘PERIOD OF COUPLING’ should be set to 1(coupled mode) and ‘BED-LOAD’ = YES.

Reference level

The reference elevation $Z_{ref}$ corresponds to the interface between bed-load and suspended-load. For flat beds, the bed-load layer thickness is proportional to the grain size, whereas when the bed is rippled, the bed-load layer thickness scales with the equilibrium bed roughness ($k_r$).

The definition of the reference elevation needs also to be consistent with the choice of equilibrium bed-concentration formula.

Vertical concentration profile

We assume here a Rouse profile for the vertical concentration distribution, which is theoretically valid in uniform steady flow conditions:

\begin{equation*} C(z)=C_{Z_{ref}} \times \left(\frac{z-h}{z}\frac{a}{a-h}\right)^R, \end{equation*}

where $R$ is the Rouse number defined by:

\begin{equation*} R=\frac{W_s}{\kappa u_*}, \end{equation*}

and $\kappa$ the von Karman constant ($\kappa = 0.4$), and $u_*$ the friction velocity corresponding to the total bed shear stress as commented before, and $a$ is the reference elevation above the bed elevation.

Ratio between the reference and depth-averaged concentration By depth-integration of the Rouse profile, the following relation can be established between the mean (depth-averaged) concentration and the reference concentration:

\begin{equation*} C_{Zref} = F C, \end{equation*}

where:

\begin{equation*} F^{-1} = \left(\frac{Z_{ref}}{h}\right)^R\int_{Zref/h}^1\left(\frac{1-u}{u}\right)^R \d u. \end{equation*}

In Sisyphe, we use the following approximate expression for $F$:

\begin{equation*} \begin{align*} F^{-1} &= \frac{1}{\left(1-Z\right)} B^R\left( 1-B^{(1-R)} \right), \quad \text{for } R \neq 1\\ F^{-1} &= -B \log B, \quad \text{for } R = 1, \end{align*} \end{equation*}

with $B = Z_{ref}/h$.

Keywords

Convection velocity

A straightforward treatment of the advection terms would imply the definition of an advection velocity and replacement of the depth-averaged velocity U by:

\begin{equation*} U_{conv} = \overline{UC}/C. \end{equation*}

A correction factor is introduced in Sisyphe , defined by:

\begin{equation*} F_{conv} =\frac{U_{conv}}{U} \end{equation*}

The convection velocity should be smaller than the mean flow velocity ($F_{conv} \leq 1$) since sediment concentrations are mostly transported in the lower part of the water column where velocities are smaller. We further assume an exponential profile concentration profile which is a reasonable approximation of the Rouse profile, and a logarithmic velocity profile, in order to establish the following analytical expression for $F_{conv}$:

\begin{equation*} F_{conv} =-\frac{I_2 - \ln\left(\frac{B}{30}\right) I_1}{I_1 \ln\left( \frac{eB}{30}\right)}, \end{equation*}

with $B=k_s/h = Z_{ref}/h$ and

\begin{equation*} I_1 =\int_B^1\left(\frac{(1-u)}{u}\right)^R \d u,\quad I_2 = \int_B^1 \ln u \left(\frac{(1-u)}{u} \right)^R \d u, \end{equation*}

see Huybrechts et al., 2010.

Keywords

Initial concentrations, boundary conditions

Initial conditions

The initial concentration for the suspended load can be either imposed within condim_susp.f or specified in the steering file through the keyword ‘Initial suspension concentrations ’ initializes the value of the volume concentration for each class.

Boundary conditions

For the boundary conditions, the concentration of each class can be specified in the steering file through keyword: ‘concentration per class at boundaries’. It may be also convenient to use keyword ‘equilibrium inflow concentration =Yes’: the concentration at the entrance of the domain and at t = 0 is set to its equilibrium value, according to the choice of the ‘reference concentration formula’. Input concentrations can be also directly specified (user subroutines: conlit.f).

Equilibrium conditions

The concentrations at the entrance of the domain can be calculated by Sisyphe assuming equilibrium conditions in order to avoid unwanted bed-evolution at the entrance of the domain, and also at the first time step, it is possible to impose the concentration to its equilibrium value, by activating the keyword ‘EQUILIBRIUM INFLOW CONCENTRATION’. The equilibrium (depth-averaged) concentration is then calculated assuming equilibrium concentration at the bed and a Rouse profile correction for the F factor.

Keywords

User subroutines

The subroutine condim_susp.f can be used to specify the initial conditions for the sediment concentration .The subroutine conlit.f can be used to specify the concentration at the entrance of the domain.


Wave effects


Theoretical background

In the coastal zone, the effect of waves, superimposed to a mean (waveinduced or tidal) current modifies the bottom boundary layer structure. Due to the reduced thickness of the bed boundary layer, bottom shear stresses are largely enhanced and resulting sand transport rates are an order of magnitude greater than in the case of currents alone. In addition, the wave-generated ripples also have a strong effect on the bed roughness and sand transport mechanisms. These effects will only be accounted for when the keyword ‘effect of waves’ is activated.

Linear wave theory

The wave orbital velocity is calculated by Sisyphe , assuming linear theory:

\begin{equation*} U_0=\frac{H_s \omega }{2 \sinh (kh)}, \end{equation*}

where $\omega = 2\pi/T_p$ is the wave pulsation, $k = 2\pi/L$ , where $L$ is the wave length.

The wave number is calculated from the wave dispersion relation:

\begin{equation*} \omega^2 = gk\tanh (kh). \end{equation*}

Wave propagation

The wave parameters (wave height, wave period, wave direction) need to be specified. They can be simply imposed in the Fortran file (user subroutine: condim_sisyphe.f ) or calculated by a previous TOMAWAC simulation [10]. The ‘HYDRODYNAMIC FILE’ can contain all the wave parameters (which is recommended for the non-steady case). The waves can also be specified from an independent wave field. The ‘wave file’ gives the name of a previous wave (TOMAWAC) results file. Only the last record will be read.

Wave parameters

The wave-parameters, required as input parameters for the morphodynamic computation, are the significant wave height $H_s$ and the peak wave period $T_p$. The wave-angle $\theta_w$ needs also to be specified in some particular applications, e.g. when using for example the Bailard's formula. (If not, the default value is $90^\circ$, which means that the direction of propagation of the wave is parallel to the $x$-axis).

If the wave field is calculated by TOMAWAC, the results hydrodynamic file should contain the following variables:

The (last time record of the) ‘wave file’ issued from previous TOMAWAC computation is read by Sisyphe . In addition, the mean current issued from previous Telemac computation is read from the ‘hydrodynamic file’. For unsteady flow conditions (time-varying velocity and wave field), the user should run Tomawac in an unsteady mode, and store all the wave and hydrodynamic relevant variables on the same input ‘hydrodynamic file’.

Keywords

User-subroutine condim_sisyphe.f

Wave-induced bottom friction

Wave-induced bottom friction

The maximum stress due to waves is calculated at each time step as a function of the wave-orbital velocity $U_w$ by use of a quadratic friction coefficient $f_w$ due to waves:

\begin{equation*} \tau_w = \frac{\rho}{2} f_w U_w^2. \end{equation*}

The wave friction factor $f_w$ is calculated as a function of relative density:

\begin{equation*} f_w = f_w \left( A_0/k_s \right), \end{equation*}

where $A_0$ is the orbital amplitude on the bed ($A_0 = U_w/\omega$), and $k_s$ the bed roughness. In Sisyphe, we use the Swart (1974) formula:

\begin{equation*} \begin{align*} f_w &= \exp \left(-6.0 + 5.2\left( \frac{A_0}{k_s} \right)^{-0.19}\right), \quad \text{if } \frac{A_0}{k_s} > 1.57\\ f_w &= 0.30, \quad \text{if } \frac{A_0}{k_s} \leq 1.57. \end{align*} \end{equation*}

Wave-current interaction

For combined waves and currents, the wave-induced (mean and maximum) bottom stresses are an order of magnitude larger than in the case of currents alone. Different models can be found in the literature to calculate the wave and current bottom stress, noted $\tau_{cw}$, as a function of the bottom shear stress due to currents only, $\tau_c$ , and the maximum shear stress due to waves only, $\tau_w$. In Bijker :

\begin{equation*} \tau_{cw} = \tau_c + \frac{1}{2} \tau_w \end{equation*}

The mean ($\tau_{mean}$) and maximum ($\tau_{\max}$) shear stress can also be calculated following Soulsby's method. Non-dimensional variables are defined:

\begin{equation*} X=\frac{\tau_0}{\tau_0+\tau_w};\quad Y=\frac{\tau_{moy}}{\tau_0+\tau_w};\quad Z=\frac{\tau_{\max}}{\tau_0+\tau_w} \end{equation*}

They can be parameterized based on different turbulence models:

\begin{equation*} \begin{align*} Y &= X\left(1+bX^p(1-X)^q\right),\\ Z &= 1+aX^m(1-X)^n \end{align*} \end{equation*}

The various model coefficients ($a, b, m, n, p, q$) are empirical functions of friction coefficients $f_w$, $C_D$ and also of the wave-current angle $\phi$:

\begin{equation*} \begin{align*} a &= (a_1 +a_2 \left|\cos \phi \right|^I)+(a_3+a_4 \left|\cos \phi \right|^I)\log_{10}\left(\frac{2f_w}{C_D} \right),\\ b &= (b_1 +b_2 \left|\cos \phi \right|^J)+(b_3+b_4 \left|\cos \phi \right|^J)\log_{10} \left(\frac{2f_w}{C_D } \right), \end{align*} \end{equation*}

and similar expressions for $m, n, p$, and $q$.

The 26 fitted constants ($a_i$, $b_i$, $m_i$, $n_i$, $p_i$, $q_i$, $I$ and $J$) depend on the turbulence model selected. The coefficients which have been used in \sisyphe correspond to the model of Huynh-Thanh et Temperville (1991) [24b], $I = 0.82$, $J = 2.7$:

$i$ $a_i$ $b_i$ $m_i$ $n_i$ $p_i$ $q_i$
1 -0.07 0.27 0.72 0.78 -0.75 0.89
2 1.87 0.51 -0.33 -0.23 0.13 0.40
3 -0.34 -0.10 0.08 0.12 0.12 0.50
4 -0.12 -0.24 0.34 -0.12 0.02 -0.28

Friction coefficient under combined waves and current

The characteristic shear-stress to represent wave current interactions, noted $\left\langle \left| \overrightarrow{\tau_{cw}} (t)\right| \right\rangle$, is related to the time averaged mean-velocity:

\begin{equation*} \left\langle \left| \overrightarrow{\tau_{cw}} (t)\right| \right\rangle =\rho f_{cw} \left\langle \left| \overrightarrow{U(t)} \right|^2 \right\rangle, \end{equation*}

where

\begin{equation*} \left\langle \left| \overrightarrow{U(t)} \right|^2 \right\rangle =U_c^2 +\frac{1}{2} U_w^2 \end{equation*}

According to Camenen (2002, page 60), $\left\langle \left| \overrightarrow{\tau_{cw}} (t)\right| \right\rangle$ is taken as a weighted average between $\tau_{moy}$ and $\tau_{\max}$:

\begin{equation*} \left\langle \left| \overrightarrow{\tau_{cw}} (t)\right| \right\rangle = X\tau_{moy} + (1-X)\tau_{\max}, \end{equation*}

which is equivalent to:

\begin{equation*} \left\langle \left| \overrightarrow{\tau_{cw}} (t)\right| \right\rangle = Y\tau_c + Z\tau_w. \end{equation*}

The final expression for the wave-current interaction factor is:

\begin{equation*} f_{cw} = \frac{Y\tau_c + Z\tau_w}{{U_c^2} + \frac{1}{2}{U_w^2}}. \end{equation*}

This expression of the wave and current friction factor is used in Bailard and in Dibajnia and Watanabe sand transport formulae.

Wave-induced ripples

Equilibrium ripples are generally observed outside the surf zone. Their dimensions can be predicted as a function of waves (orbital velocity $U_0$ and wave period $T=2\pi/\omega$), for a given uniform sediment diameter $d_{50}$. We follow in \sisyphe the procedure of Wiberg and Harris ([34b], 1994). This formulation is only applicable for oscillatory flow conditions and does not account for the effect of a superimposed mean current.

Ripples can be classified into three types depending on the value of the ratio of wave orbital diameter to mean grain diameter, $D_0/d_{50}$, where $D_0 = 2U_0/\omega$. The ripples dimensions, namely their wave length $\lambda$ and height $\eta$, can be calculated as follows:

In the orbital ripples regime, under moderate wave conditions, ripples dimensions are simply proportional to $D_0$:

\begin{equation*} \begin{align*} \lambda &= 0.62D_0 \\ \eta &= 0.17\lambda. \end{align*} \end{equation*}

For larger waves, in the anorbital ripple regime, ripples dimensions scale with the sand grain diameter:

\begin{equation*} \begin{align*} \lambda &= 535 d_{50} \\ \eta &= \lambda \exp \left(-0.095\left(\log\frac{D_0}{\eta} \right)^2 + 0.042\log\frac{D_0}{\eta}-2.28\right). \end{align*} \end{equation*}

The effect of ripples is to increase the bed roughness. As in Bijker's formula, we assume in Sisyphe:

\begin{equation*} k_s = \max \left(\eta ,3d_{50} \right). \end{equation*}

The effect of the mean current superimposed to the waves is accounted for by introducing a correction factor, noted $\alpha$, to the orbital excursion, following Tanaka and Dang (1996):

\begin{equation*} \alpha = 1 + 0.81\left(\tanh (0.3S_*^{2/3})\right)^{2.5} \left(\frac{U}{U_w}\right)^{1.9}, \end{equation*}

with

\begin{equation*} S_* = d_{50}\frac{\sqrt{(s-1)g d_{50}}}{4\nu}. \end{equation*}

Keywords

Wave-induced sand transport

The main effect of waves is to increase the sand transport rates in comparison to the case of currents alone. This effect can be accounted for by using one of the wave sand transport formula which have been programmed in Sisyphe : the Bijker’s, the Soulsby – Van Rijn formula ,Bailard’s and the Dibajnia and Watanabe formulae. In Sisyphe , the solid flow is assumed to be oriented in the direction of the mean current. The possible deviation of the transport in the direction of waves, which can be important in the near shore zone due to non-linear effects, is not accounted for.

The Bijker formula (1968) can be used for determining the total transport rate Qs (bed-load + suspension), with its two components, the bed-load Qsc and the suspended-load Qss being determined separately. The development is made on the scalar variables, since Bijker’s formula is an extension of a steady flow formula to account for the effect of the wave enhanced shear stress.
The suspended load transport is treated in a simplified manner: the concentration profile is assumed to be continuously in equilibrium. The inertia effects related to the water column loading and unloading will then in no way be modelled. In addition, no exchange takes place with the bed load layer, only the continuity of concentration is ensured in that case.
- Bed-load transport
Bijker extended the steady bed-load formula proposed by Frijlink (1952), adding the effect of the wave into it. In nondimensional variables, the bed-load transport rate is:

\begin{equation*} \Phi_b = b\theta_c^{0.5}\exp\left(-0.27\frac{1}{\mu \theta_{cw}}\right), \end{equation*}

where $\theta_c$ is the non-dimensional shear stress due to currents alone, $\theta_{cw}$ the non-dimensional shear stress due to wave-current interaction, and $\mu$ is a correction factor which accounts for the effect of ripples.
In the original formulation, we have b = 5. Recent studies showed that b = 2 provides better results in comparison to 43 field data and other engineering models. Then the formula is programmed in Sisyphe with b = 2. The ripple correction factor ì is calculated in the same way as in the Meyer Peter formula.
- Suspended load transport
The suspended load transport is obtained by vertically integrating, from a reference level corresponding to the thickness of the bed load layer, the product of the volume concentration by the mean horizontal velocity. After depth-integration and assuming a Rouse profile for the concentration and a logarithmic velocity profile for the mean velocity profile, the suspended load can be written as

\begin{equation*} Q_{ss} = Q_{sc} I, \end{equation*}

with

\begin{equation*} I=1.83\times 0.216\frac{B^{A-1}}{(1-B)^A} \int_B^1 \left(\frac{1-y}{y}\right)^A \ln\left(\frac{33y}{B}\right) \d y, \end{equation*}

and

\begin{equation*} \begin{align*} A &= \frac{W_s}{\kappa u_*},\quad u_*=\sqrt{\frac{\tau_{cw}}{\rho}}\\ B &= k_s/h. \end{align*} \end{equation*}

Keywords
* BED-LOAD TRANSPORT FORMULA: ICF= 4
* WAVE EFFECTS’ = YES
– The Soulsby - Van Rijn formula.The transport rate Qs (volume/volume) due to the combined action of waves and current is provided by the following equation:

\begin{equation*} Q_{b,s} = A_{b,s} U\left[ \left( U^2+2\frac{0.018}{C_D} U_0^2\right)^{0.5}-U_{cr}\right]^{2.4}. \end{equation*}

This formula can be applied to estimate both components of the total sand transport rate (bed-load + suspension), and it is suitable for rippled bed regimes.
$A_{b,s}$ are the bed-load/suspended load coefficients, calculated by:

\begin{equation*} \begin{align*} A_b &= \frac{0.005 h \left(d_{50}/h\right)^{1.2}}{\left((s-1)gd_{50}\right)^{1.2}} \\ A_s &= \frac{0.012 d_{50}D_*^{-0.6}}{\left((s-1)gd_{50}\right)^{1.2}}, \end{align*} \end{equation*}

where $U$ is the depth-averaged current velocity, $U_0$ is the RMS orbital velocity of waves, and $C_D$ is the quadratic drag coefficient due to current alone. This formula has been validated assuming a rippled bed roughness: $k_s = 0.18$ m
The critical entrainment velocity, $U_{cr}$ is given by:

\begin{equation*} \begin{align*} U_{cr} &= 0.19 d_{50}^{0.1}\log_{10}\left(\frac{4h}{D_{90}}\right), \quad \text{if } 0.1 mm\leq d_{50}\leq 0.5 mm \\ U_{cr} &= 8.5 d_{50}^{0.6} \log_{10}\left(\frac{4h}{D_{90}} \right), \quad \text{if } 0.5 mm\leq d_{50}\leq 2.0 mm. \end{align*} \end{equation*}

Validity range
$h = 1-20$ m, $U = 0.5-5$ms$^{-1}$, $d_{50}=0.1-2$mm.
Keywords
* BED-LOAD TRANSPORT FORMULA: ICF= 5
* WAVE EFFECTS’ = YES
In the Van-Rijn-Soulsby formula, the diameter $D_{90}$, characteristic of the coarser grains, should also be specified by use of the keyword: `D90'.

Bailard’s formula [13] The Bailard’s formula is based on an energetic approach. The bedload and the suspended-load components of the sand transport rate are expressed respectively as the third- and fourth-order momentum of the near-bed time-varying velocity field $\overrightarrow{U(t)}$

\begin{equation*} \begin{align*} Q_b &= \frac{f_{cw}}{g(s-1)} \frac{\epsilon_c}{\tan \varphi} \left\langle \left| \overrightarrow{U} \right|^2 \overrightarrow{U} \right\rangle \\ Q_s &= \frac{f_{cw}}{g(s-1)} \frac{\epsilon_s}{W_s} \left\langle \left| \overrightarrow{U} \right|^3 \overrightarrow{U} \right\rangle \end{align*} \end{equation*}

with $f_{cw}$ the friction coefficient which accounts for wave-current interactions, $\epsilon_s$, $\epsilon_c$ empirical factors ($\epsilon_s = 0.02$, $\epsilon_c = 0.1$), $\varphi$ sediment friction angle ($\tan \varphi$ = 0.63) and $<\cdot>$ time-averaged over a wave-period.

Under combined wave and currents, the time-varying velocity vector $\overrightarrow{U(t)}$ is composed of a mean component $U_c$ and an oscillatory component of amplitude $U_0$, assuming linear waves; $\phi$ is the angle between the current and the wave direction, $\phi_c$ is the angle between the mean current direction and the x-axis, $\phi_w$ is the angle between the wave direction of propagation and the x-axis.

The time-varying velocity field verifies:

\begin{equation*} \overrightarrow{U(t)} =U_c \exp^{i\phi_c} + U_0 \cos \omega t \exp^{i\phi_w}. \end{equation*}

For the third-order term, one can derive an analytical expression in the general case, whereas for the fourth-order momentum, we have to assume the waves and currents to be co-linear ($\phi = 0$), in order to find an analytical expression. (For the general case of a non-linear non-colinear wave, we would have to integrate numerically the fourth-order velocity term (cf. Camenen, 2002 [14a]). This method is however not efficient.

\begin{equation*} \begin{align*} \left\langle\left| \overrightarrow{U} \right|^2 \overrightarrow{U} \right\rangle &= \overrightarrow{U_c}(U_c^2+\frac{1}{2} U_0^2)+ \overrightarrow{U_0} (\overrightarrow{U_c}\cdot\overrightarrow{U_0})\\ \left\langle \left| \overrightarrow{U} \right|^2 \overrightarrow{U} \right\rangle_x &= \left[ U_c^3 + U_c U_0^2 (1+\frac{1}{2} \cos 2\phi )\right] \cos \phi_c -\frac{1}{2} U_c U_0^2\sin 2\phi \sin \phi_c\\ \left\langle \left| \overrightarrow{U} \right|^2 \overrightarrow{U} \right\rangle_y &= \left[ U_c^3 + U_c U_0^2 (1+\frac{1}{2} \cos 2\phi )\right] \sin \phi_c +\frac{1}{2} U_c U_0^2\sin 2\phi \cos \phi_c \end{align*} \end{equation*}

\begin{equation*} \begin{align*} \left\langle \left| \overrightarrow{U} \right|^3 \overrightarrow{U} \right\rangle_x = \frac{1}{8} (24 U_0^2 U_c^2 + 8 U_c^4 +3U_0^4)\cos \phi_c\\ \left\langle \left| \overrightarrow{U} \right|^3 \overrightarrow{U} \right\rangle_y = \frac{1}{8} (24 U_0^2 U_c^2 + 8 U_c^4 +3U_0^4)\sin \phi_c \end{align*} \end{equation*}

Keywords

The Dibajnia and Watanabe (1992) formula is an unsteady formula, which accounts for inertia effects. The sand transport rate is calculated by:

\begin{equation*} \frac{\overrightarrow{Q_s} }{W_s d_{50}} =\alpha \frac{\overrightarrow{ \Gamma} }{\left| \Gamma \right|^{1-\beta}}, \end{equation*}

with $\alpha=0.0001$ and $\beta=0.55$, empirical model coefficients. $\overrightarrow{\Gamma}$ represents here the difference between the amounts of sediments transported in the onshore and offshore directions. This formula is used to estimate the intensity of the solid transport rate, the direction is assumed to follow the mean current direction.

Time-varying velocity field

In the wave direction, the time-varying velocity is given by:

\begin{equation*} U(t)=U_c \cos \phi + U_0 \sin \omega t, \end{equation*}

$r$ is defined here as the asymmetry parameter:

\begin{equation*} r=\frac{U_0}{U_c\cos \phi}. \end{equation*}

The wave cycle is decomposed into two parts:

1. During the first part of the wave-cycle ($0 < t < T_1$), the velocities are directed onshore ($U(t) > 0$). 2. During the second part ($T_2 = T-T_1$), the velocities are negatives ($U(t) <0$).

The sand transport rate in the wave direction is:

\begin{equation*} \Gamma_{X'} = \frac{u_1T_1(\Omega_1^3+\Omega_2^{'3})-u_2T_2(\Omega_2^3+\Omega_1^{'3})}{(u_1+u_2)T} \end{equation*}

where $\Omega_1$ and $\Omega_2$ represent the amount of sand transported respectively in the onshore and offshore directions which will be deposited before flow reversal; $\Omega'_1$ et $\Omega'_2$ represent the remaining part which stays in suspension after flow eversal. The inertia of sediment depends on the ratio $d_{50}/W_s$, and represented by parameter $\omega_i$ defined by:

\begin{equation*} \omega_i = \frac{u_i^2}{2(s-1)gW_s T_i}, \end{equation*}

with

\begin{equation*} u_i^2 = \frac{2}{T_i} \int_{u>\text{or}<0} u^2 \d t \end{equation*}

1. all sediment is deposited before flow reversal $\omega_i\leq \omega_{crit}$

\begin{equation*} \Omega_i=\omega_i\frac{2W_s T_i}{d_{50}}, \quad \Omega'_i = 0 \end{equation*}

2. part of the sediment stays in suspension after flow reversal $\omega_i \geq \omega_{crit}$

\begin{equation*} \Omega_i = \omega_{crit} \frac{2W_s T_i}{d_{50}}, \quad \Omega'_i=(\omega_i-\omega_{crit})\frac{2W_sT_i}{d_{50}}. \end{equation*}

The critical value of $\omega_{crit}$ is calculated as a function of the wave-current non-dimensional friction:

\begin{equation*} \Theta_{cw} = \frac{<\tau_{cw} >}{\rho}\frac{1}{(s-1)gd_{50}} \end{equation*}

According to Temperville et Guiza (2000):

\begin{equation*} \begin{align*} \omega_{crit} &= 0.03, \quad\text{if }\Theta_{cw}\leq 0.2 \\ \omega_{crit} &= 1-\sqrt{1-((\Theta_{cw}-0.2)/0.58)^2}, \quad\text{if } 0.2\leq\Theta_{cw}\leq 0.4 \\ \omega_{crit} &= 0.8104\sqrt{\Theta_{cw}} - 0.4225, \quad\text{if } 0.4\leq\Theta_{cw}\leq 1.5\\ \omega_{crit} &= 0.7236\sqrt{\Theta_{cw}} - 0.3162, \quad\text{if } 1.5\leq\Theta_{cw}. \end{equation*}

Keywords


Sand grading effects


Sediment bed composition

Granulometry distribution

The number of classes of bed material is specified in the steering file and limited to $NSICLA < xxx$. For uniform sediment, the granulometry distribution is represented by one single value ($NSICLA =1$) for the whole domain. The mean grain diameter and corresponding settling velocity are defined in the steering file.

For sediment mixtures, the granulometry distribution is discretized in a number of classes. Each class of sediment noted $j$ ($1 < j < NSICLA$) is defined by its mean diameter $d_{50}(j)$ and volume fraction in the mixture at every node, noted $AVAI(j)$. The characteristics of each class of sediment (mean diameter, Shields number $\theta_c (i)$, settling velocities $W_s(i)$ of each class can be specified in the steering file or calculated by the model (same Van Rijn formulae than for single class sediment).

The percent of each class of material need to verify:

\begin{equation*} \sum\limits_{j=1,NSICLA}AVAI(j)=1 \quad \text{and} \quad 0\leq AVAI(j) \leq 1. \end{equation*}

The initial bed composition can be defined in the steering file for simple configuration (uniform bed). For a more complex bed configuration (spatial variation, vertical structure), the user subroutine {\texttt init\_compo.f} must be used to define the initial state.

The mean diameter ($D_m$, noted $ACLADM$) is calculated based :

\begin{equation*} D_m = \sum\limits_{j=1,NSICLA} AVAI(j) D(j). \end{equation*}

Bed structure

The bed is stratified into a number of layers ($NOMBLAY$) that enables vertical variations of the sediment bed composition. The percentage of each class $j$ of material (noted $AVAI(i,j,k)$ as well as the thickness of each layer $E_s$ are defined at each point $i$ and for each layer $k$.

The composition of transported material is computed using the composition of the upper layer (see below the definition of an `active layer').

The initial composition of the bed (number of layers, thickness of each layer $E_s$, composition of each layer $AVAI$ are specified in user-subroutine init_compo.f.

The subroutine init_avail.f checks that the structure of the initial bed composition compatible with the position of the rigid bed (as defined in user subroutine noerod.f):

\begin{equation*} Z_f(i)-Z_r(i) = \sum_{k=1,NOMBLAY} E_s(k). \end{equation*}

In general, the thickness of the bed is very large (by default, $100$ m), so the bottom layer thickness need to be increased.

Assuming the porosity of each class to be identical (parameter $XKV =1/(n-1)$ is a constant), the total mass of sediments per unit area is:

\begin{equation*} M_s(i) = \sum_{k=1,NOMBLAY}\rho_s E_s(k)(n-1). \end{equation*}

In each layer, the AVAI (percent of each class) which are defined as the volume of each class of material per the total volume of material, are considered to be constant. In each layer (and at each point) the constraints on the AVAI need to be satisfied :

\begin{equation*} \sum_{j=1,NSICLA} AVAI(i,k,j) = 1\quad \text{and}\quad 0 \leq AVAI(i,k,j) \leq 1. \end{equation*}

Keywords
Initial sediment bed composition is defined by:

For more than one size-classes, the user should specify NSICLA values for the mean diameter, initial fraction. . . For example, for a mixture of two sands

If not specified the settling velocities and Shield parameters are calculated as a function of grain size for each particles.

Transport of sediment mixtures

The method programmed in SISYPHE for the treatment of mixtures of sediment with variable grain sizes is classical and based on previous models from the literature (as for example the 1D model Sedicoup developed at Sogreah). There are two layer concepts implemented in Sisyphe. The active layer model based on Hirano’s concept ([XX], 1971) (default) and a newly developed continues vertical sorting model (C-VSM) from Merkel ([XX], 2012), ([XX], 2012).

Active layer/Active stratum

Since only a thin upper layer will be transported, the upper layer is therefore subdivided into a thin ‘active layer’ (k = 1) and an ‘active stratum’ (k = 2). The composition of the active layer is used to calculate the composition of transported material and the intensity of transport rates for each class of sediments:

\begin{equation*} Q_s =\sum_{j=1,NSICLA} AVA0(j) Q_s(j), \end{equation*}

where $AVA0(j)$ is the percentage of class $j$ in the active layer.

The active layer thickness is not well defined and depend on the flow and sediment transport parameters (cf. Van Rijn, 1987). According to Yalin (1977), it is proportional to the sand diameter of the upper layer ($1.2d_{90}$).

The active later thickness generally scales with the bed roughness, which is typically of the order of a few grain diameters for flat beds up to a few cm in the case of rippled beds. In the presence of dunes, the active layer thickness should be half the dune height (cf. Ribberink, 1987).

In Sisyphe, the active layer thickness is an additional model parameter (noted $ELAY0$) which can be estimated by calibration and specified by the user in the steering file by keyword `ACTIVE LAYER THICKNESS'. With the option `CONSTANT ACTIVE LAYER THICKNESS' = NO, it is possible to use a time-and spatially varying active layer thickness during the simulation. In Sisyphe, we simply assume $ELAY0 = 3 D_m$ ($D_m$ is the mean diameter of the upper layer).

The erosion rate is restricted by the total amount of sediments in the active layer, which therefore acts as a rigid bed. The same methods applied for rigid beds are applicable for active layer formulation.

In order to avoid numerical problems such as negative sediment fractions, it is advised to use finite elements combined with the positive depth algorithm, as explained in \S 4.5 for the treatment of rigid bed, when dealing with graded sediment and a thin-sized active layer. The error message `TIME STEP SHOULD BE REDUCED' also appears in the listing file when the erosion is greater than half of the active layer thickness. It is stongly recommended to follow the message.

Continuous vertical sorting model

The C-VSM overcomes many limitations of the classic layer implementation of Hirano (HR-VSM). Even though it is a different way to manage the grain sorting, it is just another interpretation of Hiranos original idea with fewer simplifications. So still an active layer is defined. But the grain distribution in this layer is computed each time step from a depth dependent bookkeeping model for each grain size fraction with unlimited numerical resolution for each horizontal node. The new model doesn’t overcome the need to carefully calibrate the same input parameters as all other models, but the new interpretation has the following advantages:

* It is possible to keep minor but prominent grain mixture variations even after a high number of time steps. Smearing effects and bookkeeping accuracy is defined by user defined thresholds or the computational resources, rather than through a fix value.

* Various functions for the impact depth of the shear stress can be chosen to the projects demands. The result is a much more accurate vertical grain sorting, which results in better prognoses for bed roughness, bed forms and erosion stability.

A couple of new keywords must be set in your sis.cas file. The full C-VSM output can be found in the new Selafin files VSPRES & VSPHYD in the tmp-folders. As the higher resolution of the C-VSM needs resources, you can reduce the printout period, or suppress the output at all. The common SISYPHE result files only show the averaged values for the active layer. Even more disk space can be saved, if only few points are printed out as .VSP.CSV files in the subfolder /VSP/. We recommend using between 200 and 1000 vertical sections. More will not improve the accuracy much, and less will lead to increasing data management, as the profile compression algorithms are called more often.

Keywords

\begin{equation*} ALT = 5 d_{MAX} \end{equation*}

\begin{equation*} ALT =\frac{2 \tau_B}{(1-n) g (\rho_S - \rho) \tan \Phi} \end{equation*}

\begin{equation*} ALT =0.3 D_*^{0.7}\frac{\tau_B-\tau_C}{\tau_B}^{0.5} d_{50} \end{equation*}

\begin{equation*} ALT =5 \frac{\tau_B}{(\rho_S - \rho) g d} - 0.0549)^{0.56} d_{50} \end{equation*}

\begin{equation*} ALT =\frac{d_{90}}{1-n} \max(1,\frac{\tau_B}{\tau_C}) \end{equation*}

\begin{equation*} ALT =3 d_{50} \end{equation*}

General formulation

Bed-load transport rates are computed separately for each class using classical sediment transport formulae, corrected for sand grading effects (see below). The Exner equation is solved for each class of material. The individual bed evolutions due to each class of bed material are then added to give the global evolution due to bed-load.

Similarly, the suspended transport equation is solved for each class of material and the resulting bed evolutions for each class are added to give the global evolution due to the suspended load. The concentration of each class of sediments is computed, with the corresponding settling velocity, erosion and deposition flux.

At each time step, the bed evolution (due to bed load and to suspension) is used to compute the new structure (composition and layer thickness) of the new sediment bed mixture. The algorithm which determines the new bed composition considers two possibilities (deposit or erosion) and ensures mass conservation for each individual class of material. The algorithm is explained in Gonzales (2002).

Hiding-exposure

In order to calculate the bed-load transport rate for each class of sediment, it is necessary to take into account the effect of hiding-exposure: in a mixture, the big particles will be more exposed to the flow than the smaller ones which will be protected. Classical sand transport formulae (for uniform sand) need to be corrected by use of a “hiding-exposure” correction factor (HIDFAC, as defined below in §6.4.

Keywords
Sand grading effects are defined by:

Remarks

Sediment transport formulae

The formula of Hunziker ([24], 1995) is an adaptation of the Meyer-Peter Mu”ller formula for fractional transport. The volumetric sediment transport per class is given by:

\begin{equation*} \Phi_b = 5 p_i \left[\zeta_i\left( \mu \theta_i -\theta_{cm} \right) \right]^{3/2} \quad\text{if}\quad \mu \theta_i > \theta_{cr}, \end{equation*}

with $p_i$ the fraction of class $i$ in the active layer, $\theta_i$ the Shields parameter,$\theta_{cm} = \theta_{cr}\left(D_{mo}/D_m\right)^{0.33}$ the corrected critical Shields parameter, $\theta_{cr}$ the critical Shields parameter ($\theta_{cr}=0.047$), $\xi_{i} =\left(D_i/D_m\right)^{-\alpha}$ the hiding / exposure factor, $D_i$ the grain size of class $i$, $D_m$ the mean grain size of the surface layer (m), and $D_{mo}$ the mean grain size of the under layer (m).

Here the critical Shields parameter is calculated, as a function of the dimensionless mean grain size $D_*$. It should be noted that according to Hunziker, stability problems may occur outside the parameter range:

\begin{equation*} \alpha \leq 2.3\quad\text{and}\quad \frac{D_i}{D_m} \geq 0.25. \end{equation*}

Keywords
Sand grading effects are defined by:

Hiding/exposure factors

Sediment transport rates of each class can be calculated using classical sediment transport formulae. Those formulation are initially valid for uniform particles and need to be modified to account for sand grading effects: in a mixture, big particles are more exposed to the flow than if they are alone, and small particles more protected. The bed-load formulae can be modified by use of a hiding/exposure factor.

Different formulations have been implemented in SISYPHE, depending on the choice of ‘HIDING FACTOR FORMULA’ (HIDFAC).

Two formulae (Egiazaroff and Ashida & Michiue) have been written in association with the Meyer-Peter and Müller formula. The third formula (Karim et al)can be associated with any bed-load predictor. The first two formulae modify the critical Shields parameter that will be used in the Meyer-Peter formula, whereas the third formula modify directly the bed-load transport rate.

\begin{equation*} \begin{align*} \theta_{cr} &= \zeta_i 0.047, \\ \zeta_i &= \left[\frac{\log (19)}{\log (19 D_i /D_m)} \right]^2 \end{align*} \end{equation*}