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 bedload 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 depthaveraged 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 noncohesive sediments (uniform or graded),
cohesive sediments as well as to sandmud 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 bedroughness 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 (Telemac2d , Telemac3d 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
 Minimum set of input files: CLI, SLF,CAS…
 Additional/optional input files: F, REF, RES,…
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 and output files
 Physical parameters (sand diameter, cohesive or not, settling velocity, etc.)
 Main sediment transport processes (transport mechanism and formulae, etc.)
 Additional sediment transport processes (secondary currents, slope effect, etc.)
 Numerical options and parameters (numerical scheme, options for solvers, etc.)
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 (threenode) 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 (Telemac2d or Telemac3d ). 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 (Telemac2d , Telemac3d or Tomawac ).
Boundary conditions file
The format of the boundary condition file is the same as for Telemac2d or Telemac3d. 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 Telemac2d .
Obligatory keywords
Names of principal input/output files (other files are optional)
 BOUNDARY CONDITIONS FILE
 GEOMETRY FILE
 FORTRAN FILE
 RESULTS FILE
When combining a previous computation:
 COMPUTATION CONTINUED (=NO :default option)
 PREVIOUS SEDIMENTOLOGICAL COMPUTATION FILE
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
 Principle
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 timescales 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.
 Flow updating
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.
 Mass continuity
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.
 Steering/fortran files
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 usersubroutine (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:
 TIME STEP, NUMBER OF TIME STEPS
 GRAPHIC PRINTOUT PERIOD
 LISTING VARIABLES FOR GRAPHIC PRINTOUTS
For hydrodynamics (imposed flow and updated):
 HYDRODYNAMIC FILE
 STEADY CASE (=NO, default option)
 TIDE PERIOD (= 44640, default option)
 STARTING TIME OF THE HYDROGRAM ( =0., default option)
 NUMBER OF TIDES OR FLOODS (= 1, default option)
 CRITICAL EVOLUTION RATIO (=0.1, default value)
For waves:
 WAVE FILE, WAVE EFFECTS
Internal coupling
 Principle
For a more realistic simulation, Sisyphe can be automatically coupled (internally) with the hydrodynamic model, Telemac2d 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.
 Time step and coupling period
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 recalculated 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 timestep (‘coupling period = 1’).
The time step of Sisyphe is equal to the time step of Telemac2d or – 3D multiplied by the ‘coupling period’. The graphic/listing printout periods are the same as in the Telemac computation.
The Telemac2d /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 Telemac2d (or –3D) steering files:
 COUPLING WITH = SISYPHE
 COUPLING PERIOD (=1, default value)
 NAME OF SISYPHE STEERING FILE
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:
 TIME STEP
 GRAPHIC PRINTOUT PERIOD
 LISTING PRINTOUT PERIOD
 LAW OF BOTTOM FRICTION
 COEFFICIENT
Flowsediment interactions
Sediment and fluid properties
Fine sediment particles of diameters less than present complex cohesive properties which affect the sediment transport processes. For noncohesive sediments (median diameter ), the grain diameter and grain density are the key parameters which determine its resistance to erosion and sediment transport rate.
For cohesive sediments (), the grain diameter is no longer the key sediment parameter: The settling velocity now depends on the concentration and other physicochemical 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, noncohesive sediments characterized by one single value for the grain size and grain density with 2650 kg\,m which can be transported both as bedload and suspended load.
Keywords
The physical properties of the sediment are always defined in the
Sisyphe steering file using the following keywords:
 COHESIVE SEDIMENTS (= NO , default option)
 NUMBER OF SIZE CLASSES OF BED MATERIAL (NSICLA= 1, default option)
 SEDIMENT DIAMETERS (m for noncohesive sediments)
 SEDIMENT DENSITY ( kgm, default value)
 SETTLING VELOCITIES
The physical properties of the fluid are defined by:
 VISCOSITY OF THE FLUID ( 10 ms, by default)
 DENSITY OF THE FLUID (1000 kgm, by default)
 GRAVITY ACCELERATION ( ms)
Bed shear stress
Hydrodynamic model
The currentgenerated 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 Telemac2d , the bed shear stress t0 is calculated at each time step from the depthaveraged mean flow velocity by use of a quadratic friction coefficient Cd as follows:
where is the fluid density, and the depthaveraged 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 depthaveraged velocity, except if a Nikuradse friction law is applied. In that case, the friction velocity u, defined by
is related to the near bed flow velocity ) by a logarithmic velocity profile:
where is the von K\'arm\'an constant, is the vertical distance from a rough boundary, expressed as a function of the Nikuradse bed roughness (), with the grain roughness height. For flat beds, the roughness height has been shown to be approximately
Role of bed forms
A natural sediment bed is generally covered with bed forms (length and height ). 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 megaripples. The dimensions of dunes scale with the water depth (, ). In most cases, large scale models do not resolve the small to medium scale bed forms (ripples, megaripples) which need therefore to be parameterized by an increased friction coefficient:
The local skin friction component determines the bedload 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 currentinduced 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 nearbed layer, where interparticle interactions occur. The suspended load is defined as the depthintegrated 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 includes both a bedload and suspended load :
Shields parameter
Most sand transport formulae assume a critical value of the bed shear stress, noted , for the onset of erosion: for . The Shields parameter or nondimensional critical shear stress is defined by:
can be either specified in the \sisyphe steering file by use of keyword SHIELDS PARAMETERS, or calculated by the model as a function of nondimensional grain diameter. The Van Rijn formula has been programmed
where the nondimensional diameter is defined by:
Keywords
 BEDLOAD (= YES : default option)
 SUSPENDED LOAD (=NO : default option)
 SHIELDS PARAMETER’ (default option: Van Rijn formula)
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 3Dmodels, the bed shear stress is directly related to near bed velocity assuming a log profile (eq. 2), whereas in 2Dmodels, 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 depthaveraged velocity in Sisyphe (alone or internally coupled with Telemac2d ). When Sisyphe is internally coupled to Telemac3d , 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 TELEMAC2D
The depthaveraged 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:
 deviation correction due to sloping bed effects
 secondary currents due to river meandering
Coupling with TELEMAC2D
When the model is coupled with Telemac2d , the values of the
friction coefficients (and therefore, the bed shear stress) are provided
by Telemac2d .
Coupling with TELEMAC3D
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 Telemac2d 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:
 Chézy coefficient Ch (KFROT = 2)
 Strickler coefficient (KFROT = 3)
 Manning friction (KFROT = 4)
 Nikuradse bed roughness (KFROT = 5)
Keywords
In the Sisyphe steering file, the total bed shear stress is calculated
based on:
 LAW OF BOTTOM FRICTION (KFROT =3, default option)
 BOTTOM FRICTION COEFFICIENT (St=50, default value)
Similar keywords are available in both Telemac2d 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:
Physically, the skin bed roughness should be smaller than the total bed roughness (i.e. ). 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 timevarying 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 can be admitted. Different methods are programmed in \sisyphe in order to calculate the bedform correction factor , according to the keyword `Skin friction correction'(ICR) :
 ICR = 0: no correction: the total friction issued from Telemac is directly used for sand transport calculations (=1).
 ICR = 1: the skin roughness is assumed to be proportional to the sand grain diameter like in the case of flat beds (). The proportionality coefficient is specified by the keyword ‘Ratio between skin friction and mean diameter’ (KSPRATIO).
where et are both quadratic friction coefficients related respectively to total friction and skin friction: Cd is obtained from Telemac2d or Telemac3d and is calculated from :
 ICR = 2: we use a bedform predictor to calculate the bedform roughness . Both and should influence the transport rates. Here we assume:
where the quadratic friction due to bedforms is calculated as a function of .
Keywords The option selected for the skin friction correction is based on keywords:
 SKIN FRICTION CORRECTION (ICR=0, default option)
 RATIO BETWEEN SKIN FRICTION AND MEAN DIAMETER (KSPRATIO=3,default value)
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 KFROT = 0 : the bed is assumed to be flat
 For KFROT = 1: the bed roughness is decomposed into different components (skin, bedforms, etc.) and each component is calculated by Sisyphe depending on the flow regime.
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 waveinduced bedform bed roughness is calculated as a function of the waveinduced bedform height :
For currents only, the van Rijn predictor has been implemented. The total bed roughness can be decomposed into a grain roughness , a smallscale ripple roughness , a megaripple component , and a dune roughness :
Both small scale ripples and grain roughness have an influence on the sediment transport laws, while the megaripples 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 :
where
Keywords
 LAW OF BOTTOM FRICTION (KFROT =0 or 1, flat bed or rippled bed for total friction)
 SKIN FRICTION CORRECTION (ICR=1 or 2 flat bed or rippled bed , for skin friction correction)
Bedload transport
Exner equation
Equilibrium conditions
Bedload occurs in a thin nearbed highconcentrated 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 bedload 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:
where is the bed porosity ( for non cohesive sediment), the bottom elevation, and (ms) the solid volume transport (bedload) 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 bedload only. However, it can be extended to total load (including the suspended load) assuming quasisteady and uniform flow conditions. Otherwise, a different treatment should be applied to both bed load and suspended load, in order to take nonequilibrium flow conditions into account properly.
Keywords
 BEDLOAD = yes (default option)
 SUSPENDED LOAD = No (default option)
 COEFFICIENT FUNCTION OF THE POROSITY (XKV=1/(1 n)=1.6, default value)
Bedload transport formulae
For currents only (no wave effects), a large number of semiempirical formulae can be found in the literature to calculate the bedload transport rate. Sisyphe offers the choice between 10 different bedload formulae including the MeyerPeter and Müller, EngelundHansen and Einstein Brown formulae. Most sediment transport formulae assume threshold conditions for the onset of erosion (e.g. MeyerPeter 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 nondimensional currentinduced sand transport rate, noted , can be expressed as a function of the nondimensional skin friction or Shields parameter, noted , as defined by:
, the correction factor for skin friction
Keywords
 BEDLOAD TRANSPORT FORMULA (default option is MP formula, ICF=1)
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(usersubroutine), 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.
 ICF= 1: MeyerPeter formula This classical bedload formula has been validated for coarse sediments in the range ( mm mm). It is based on the concept of initial entrainment:
where is the critical Shields parameter ( by default).
 ICF = 2: EinsteinBrown formula. This bedload formula is recommended for gravel ( mm) and large bed shear stress . The solid transport rate is expressed as:
with
and
 ICF= 3 or 30: EngelundHansen formula. The EngelundHansen formula predicts the total load (bed load plus suspended load). It is recommended for fine sediments, in the range ( mm mm but beware that the use of a total load formula is only suitable under equilibrium conditions (quasi steady and uniform flow). Two different forms of the same equation are programmed in Sisyphe:
ICF = 30 corresponds to the original formula, where the transport rate is related to the skin friction without threshold:
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:
the dimensionless bed shear stress is calculated as a function of the dimensionless skin friction :
Notice
In the case of coupling between bedload and suspended load (SUSPENSION = YES), a total load formula (N3 and 30) should not be used, in order to avoid the suspended load to be calculated twice.
 ICF = 7: Van Rijn formula. This formula was proposed by Van Rijn (1984a) in order to calculate the bedload transport rate for particles of size mm mm:
Validity range
Most sediment transport formulae are based on experiments performed under
steadyflow conditions. They demonstrate rapid variation of the bedload
transport rates as a function of the mean flow intensity: an increase of the
current velocity by will result, depending on the formula being used,
in an increase of the transport rate by over (MeyerPeter),
(EngelundHansen) or almost (EinsteinBrown). 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  MeyerPeter  EinsteinBrown  EngelundHansen  van Rijn 

IFC  1  2  3 or 30  7 
Mode of transport  Bedload  Bedload  Totalload  Bedload 
Sand diameter  > 1 mm  0.23 mm  0.24 mm  0.22 mm 
Rippled bottoms  Yes  Yes  Yes  XXX 
Low bedload flows  Yes  Yes  No  XXX 
High bedload 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 bedevolution 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
 For SLOPEFF=1. The first correction method is based on the Koch and Flokstra’s formula (1981). The solid transport rate intensity Qb0 is multiplied by a factor:
where is the coordinate in the current direction and 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 bedevolution equation. It tends to smoothe the results and reduces unstabilities.
 For SLOPEFF=2, the method of Soulsby (1997) is programmed: The threshold bed shear stress is modified as a function of bed slope and depends on the ‘FRICTION ANGLE OF THE SEDIMENT'
where is the angle of the bed to the horizontal plane and 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:
where is the direction of solid transport in relation to the flow direction, is the direction of bottom stress in relation to the flow direction, and 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:
  DEVIA = 2, the deviation is calculated based on Talmon et al. (1995) and depends on the empirical coefficient .
This default value can be modified by changing the keyword: `PARAMETER FOR DEVIATION' = .
Keywords
 SLOPE EFFECT (=YES, default option)
 FORMULA FOR SLOPE EFFECT (SLOPEFF=1, default option)
 FORMULA FOR DEVIATION (DEVIA=1, default option)
 FRICTION ANGLE OF THE SEDIMENT (, default value)
 BETA ((, default value)
 PARAMETER FOR DEVIATION ((, default value)
Sediment slide
An iterative algorithm prevents the bed slope from becoming greater than the maximum friction angle (). 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.
Keywords
 SEDIMENT SLIDE (= NO, default option)
 FRICTION ANGLE OF SEDIMENT (, default value)
Secondary currents
In meandering streams, centrifugal inertial forces are induced by the curvature of the streamlines. The depthaveraged mean pressure gradient due to the crosssectional variation of the free surface balances the depth averaged centrifugal force, noted F:
where is the crosssectional axis and the crosssectional variation in the free surface. The mass balance in the ydirection gives a zerodepth averaged crosssectional velocity.
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 crosssectional circulation, with a crossvelocity directed near the inner bank, and flow velocity therefore follows a spiral motion which cannot be represented by a depthaveraged 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 bottom shear stress is constant in a crosssection
 the bed roughness is also constant
 the mean water depth h is constant
The near bed angular deviation of the main flow can be shown to be proportional to , as explained in Yalin and Ferreira da Silva (2001). We use here . The bend radius is generally unknown, but can be inferred from the cross sectional variation in the free surface:
Keywords
 SECONDARY CURRENTS (=NO, default option)
Notice
This option should only be activated when the flow is calculated by
a depthaveraged model.
Treatment of rigid beds
Different methods
Nonerodible 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 shallowwater continuity equation.
With finite volumes, two methods are implemented in Sisyphe , depending on the ‘OPTION FOR THE TREATMENT OF NON ERODABLE BEDS’ :
 Option (0, default option), the bed is movable everywhere. The new algorithm for minimum water depth is used
 Options (4) are massconservative and avoid erosion of the rigid beds.
Keywords
 OPTION FOR THE TREATMENT OF NON ERODABLE BEDS (=0, default option)
 MASSLUMPING (=YES, default option)
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:
 In the default option (OPTBAN = 1), equations are solved everywhere. The water depth is set to a minimum value HMIN (keyword MINIMAL VALUE OF THE WATER HEIGHT, default value 1.D3) in order to avoid divisions by zero. However, in case of coupling with Telemac2d or Telemac3d , the original depth is kept for suspension, because a strict compatibility with SaintVenant continuity equation is required which would be spoiled by clipping. This is the recommended option.
 Another option (OPTBAN = 2), removes from the computation elements with points that have a depth less than a minimum HMIN. This option must be considered obsolete, it is actually not exactly massconservative and does not work properly in parallel mode.
Keyword
 TIDAL FLATS (=YES , default option)
 MINIMAL VALUE OF THE WATER HEIGHT (HMIN = 10−3m ,default value)
 OPTION FOR THE TREATMENT OF TIDAL FLATS (OPTBAN = 1, default option)
Suspended load
Suspended load transport equation
Passive scalar hypothesis
For nonequilibrium flow conditions, the bedload and the suspended load are dealt with separately treatment. The interface between the bedload and suspended load is located at :
 In the thin highconcentrated bedload layer () interparticle interactions and flowturbulent interactions strongly modify the flow structure. Equilibrium conditions are however a reliable assumption to relate the bedload to the current induced bed shear stress.
 In the upper part of the flow, for dilute suspension (), clear flow concepts still apply, and the sediment grains can be regarded as a passive scalar which follows the mean and turbulent flow velocity, with an additional settling velocity term.
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 semiempirical formula. The Van Rijn formula which is valid for noncohesive spherical particles and dilute suspensions, has been programmed in Sisyphe (cf. Van Rijn, 1993, cf. [31], p. 3.13):
where is the relative density and is gravity.
Outside the bedload layer sediment particles are assumed to follow the mean flow velocity, noted , with an additional vertical settling velocity, noted .
3D transport equation
The 3D advectiondiffusion equation for the suspended sediment concentration is:
where is the sediment diffusivity coefficient. The following boundary conditions must be satisfied both at the free surface () and at the interface between the bedload and the suspended load ():
where is the free surface elevation the position of the interface between the bedload layer and the suspension. 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 Telemac3d.
2D transport equation
The 2D transport equation for the depthaveraged mean concentration C is obtained by depthintegration :
with:
where is the water depth, assuming the bedload layer thickness to be small. After simplification of the advection terms and using the continuity equation, the following approximate depthaveraged transport equation can be solved in its nonconservative form:
where and are the depthaveraged convective flow velocities in the x and y directions.
 Convection: The convective velocity can be corrected from the depthaveraged mean velocity in order to account for the fact that most sediment is transported near the bed.
 Diffusion The diffusion terms can be set to zero (‘DIFFUSION’ = NO) According to the choice of the parameter ‘OPDTRA’ (key word), the diffusion terms can be simplified, and equation becomes
The value of the dispersion coefficient 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. ( and ).
Different schemes are available for solving the nonlinear 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
 SUSPENSION (= NO by default)
 DIFFUSION (= YES, by default)
 OPTION FOR THE DIFFUSION OF TRACER (OPDTRA=1,the default option corresponds to the simplified expression of the diffusion term)
 OPTION FOR THE DISPERSION (=1, default option)
 DISPERSION ALONG THE FLOW ( ms, default value)
 DISPERSION ACROSS THE FLOW ( ms, default value)
 TYPE OF ADVECTION (RESOL=1 , method of characteristics by default)
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 depthaveraged velocity and no longer satisfies the shallowwater continuity equation (‘correction on convection velocity’ =YES ).
Bed evolution
Bed evolution due to the suspension
The bottomevolution equation due to the suspension accounts for the vertical fluxes at the interface between the bedload and the suspended load. The bedevolution equation becomes:
where is the bottom elevation, the interface between the bedload and suspendedload, and the bed porosity. Erosion and deposition fluxes are calculated at the bedload/suspended load interface ().
Units consistency
The net sediment flux ED is in m s. For consistency with the bedload units (ms), 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 : (kgm) = .
Keywords
 MASS CONCENTRATION (NO, default option)
 BEDLOAD =NO (YES, default option)
 SUSPENDED LOAD =YES (NO, default option)
 ‘coefficient function of the porosity’
Equilibrium concentrations
Erosion/Deposition rate
For noncohesive sediments, the net erosion minus deposition flux is determined based on the concept of equilibrium concentration (Celik and Rodi, 1988):
where is the equilibrium nearbed concentration and is the nearbed concentration, calculated at the interface between the bedload and the suspended load, . Two different semiempirical formula from the literature are programmed in Sisyphe :
 ICQ= 1: Zyserman and Fredsoe formulation:
where is the critical Shields parameter and ,the nondimensional skin friction which is related to the Shields parameter, by use of the ripple correction factor.
 ICQ= 2: Bijker (1992) formula: The equilibrium concentration corresponds to the volume concentration at the top of the bedload layer. It can related to the bed load transport rate:
with an empirical factor b = 6.34.
For ICQ= 2…
For ICQ= 2, both bedload and suspended load should be calculated at each time step. The ‘PERIOD OF COUPLING’ should be set to 1(coupled
mode) and ‘BEDLOAD’ = YES.
Reference level
The reference elevation corresponds to the interface between bedload and suspendedload. For flat beds, the bedload layer thickness is proportional to the grain size, whereas when the bed is rippled, the bedload layer thickness scales with the equilibrium bed roughness ().
The definition of the reference elevation needs also to be consistent with the choice of equilibrium bedconcentration formula.
 ICQ= 1: According to Zyserman and Fredsoe, the reference elevation should be set to . In \sisyphe we take , where is the skin bed roughness for flat bed (), the default value of proportionality factor is KSPRATIO =3).
 ICQ= 2: According to Bijker, .
Vertical concentration profile
We assume here a Rouse profile for the vertical concentration distribution, which is theoretically valid in uniform steady flow conditions:
where is the Rouse number defined by:
and the von Karman constant (), and the friction velocity corresponding to the total bed shear stress as commented before, and is the reference elevation above the bed elevation.
Ratio between the reference and depthaveraged concentration By depthintegration of the Rouse profile, the following relation can be established between the mean (depthaveraged) concentration and the reference concentration:
where:
In Sisyphe, we use the following approximate expression for :
with .
Keywords
 REFERENCE CONCENTRATION FORMULA (ICQ)
Convection velocity
A straightforward treatment of the advection terms would imply the definition of an advection velocity and replacement of the depthaveraged velocity U by:
A correction factor is introduced in Sisyphe , defined by:
The convection velocity should be smaller than the mean flow velocity () 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 :
with and
see Huybrechts et al., 2010.
Keywords
 CORRECTION ON CONVECTION VELOCITY’ (= NO, default option)
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 bedevolution 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 (depthaveraged) concentration is then calculated assuming equilibrium concentration at the bed and a Rouse profile correction for the F factor.
Keywords
 INITIAL SUSPENSION CONCENTRATIONS (CS0 = 0, default value)
 EQUILIBRIUM INFLOW CONCENTRATION (=NO, default option)
 CONCENTRATION PER CLASS AT BOUNDARIES (=0, default value)
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 wavegenerated 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:
where is the wave pulsation, , where is the wave length.
The wave number is calculated from the wave dispersion relation:
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 nonsteady 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 waveparameters, required as input parameters for the morphodynamic computation, are the significant wave height and the peak wave period . The waveangle 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 , which means that the direction of propagation of the wave is parallel to the axis).
If the wave field is calculated by TOMAWAC, the results hydrodynamic file should contain the following variables:
 HM0: significant wave height (corresponds to )
 TPR5: peak period (corresponds to )
 DMOY: mean wave direction, relative to axis (corresponds to )
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 (timevarying 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
 WAVE EFFECTS (=NO, by default)
 WAVE FILE
Usersubroutine condim_sisyphe.f
Waveinduced bottom friction
Waveinduced bottom friction
The maximum stress due to waves is calculated at each time step as a function of the waveorbital velocity by use of a quadratic friction coefficient due to waves:
The wave friction factor is calculated as a function of relative density:
where is the orbital amplitude on the bed (), and the bed roughness. In Sisyphe, we use the Swart (1974) formula:
Wavecurrent interaction
For combined waves and currents, the waveinduced (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 , as a function of the bottom shear stress due to currents only, , and the maximum shear stress due to waves only, . In Bijker :
The mean () and maximum () shear stress can also be calculated following Soulsby's method. Nondimensional variables are defined:
They can be parameterized based on different turbulence models:
The various model coefficients () are empirical functions of friction coefficients , and also of the wavecurrent angle :
and similar expressions for , and .
The 26 fitted constants (, , , , , , and ) depend on the turbulence model selected. The coefficients which have been used in \sisyphe correspond to the model of HuynhThanh et Temperville (1991) [24b], , :
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 shearstress to represent wave current interactions, noted , is related to the time averaged meanvelocity:
where
According to Camenen (2002, page 60), is taken as a weighted average between and :
which is equivalent to:
The final expression for the wavecurrent interaction factor is:
This expression of the wave and current friction factor is used in Bailard and in Dibajnia and Watanabe sand transport formulae.
Waveinduced ripples
Equilibrium ripples are generally observed outside the surf zone. Their dimensions can be predicted as a function of waves (orbital velocity and wave period ), for a given uniform sediment diameter . 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, , where . The ripples dimensions, namely their wave length and height , can be calculated as follows:
In the orbital ripples regime, under moderate wave conditions, ripples dimensions are simply proportional to :
For larger waves, in the anorbital ripple regime, ripples dimensions scale with the sand grain diameter:
The effect of ripples is to increase the bed roughness. As in Bijker's formula, we assume in Sisyphe:
The effect of the mean current superimposed to the waves is accounted for by introducing a correction factor, noted , to the orbital excursion, following Tanaka and Dang (1996):
with
Keywords
 BOTTOM FRICTION LAW kfrot= 1 (3, default option)
Waveinduced 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 nonlinear effects, is not accounted for.
 Bijker formula
 The Bijker formula (1968) can be used for determining the total transport rate Qs (bedload + suspension), with its two components, the bedload Qsc and the suspendedload 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.
  Bedload transport
 Bijker extended the steady bedload formula proposed by Frijlink (1952), adding the effect of the wave into it. In nondimensional variables, the bedload transport rate is:

 where is the nondimensional shear stress due to currents alone, the nondimensional shear stress due to wavecurrent interaction, and 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 depthintegration 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

 with

 and

 Keywords
 * BEDLOAD 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:

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

 where is the depthaveraged current velocity, is the RMS orbital velocity of waves, and is the quadratic drag coefficient due to current alone. This formula has been validated assuming a rippled bed roughness: m
 The critical entrainment velocity, is given by:

 Validity range
 m, ms, mm.
 Keywords
 * BEDLOAD TRANSPORT FORMULA: ICF= 5
 * WAVE EFFECTS’ = YES
 In the VanRijnSoulsby formula, the diameter , characteristic of the coarser grains, should also be specified by use of the keyword: `D90'.
 Bailard’s formula
Bailard’s formula [13] The Bailard’s formula is based on an energetic approach. The bedload and the suspendedload components of the sand transport rate are expressed respectively as the third and fourthorder momentum of the nearbed timevarying velocity field
with the friction coefficient which accounts for wavecurrent interactions, , empirical factors (, ), sediment friction angle ( = 0.63) and timeaveraged over a waveperiod.
Under combined wave and currents, the timevarying velocity vector is composed of a mean component and an oscillatory component of amplitude , assuming linear waves; is the angle between the current and the wave direction, is the angle between the mean current direction and the xaxis, is the angle between the wave direction of propagation and the xaxis.
The timevarying velocity field verifies:
For the thirdorder term, one can derive an analytical expression in the general case, whereas for the fourthorder momentum, we have to assume the waves and currents to be colinear (), in order to find an analytical expression. (For the general case of a nonlinear noncolinear wave, we would have to integrate numerically the fourthorder velocity term (cf. Camenen, 2002 [14a]). This method is however not efficient.
 Thirdorder term:
 Fourthorder term:
Keywords
 BEDLOAD TRANSPORT FORMULA: ICF= 8
 WAVE EFFECTS’ = YES
 Dibajnia and Watanabe formula
The Dibajnia and Watanabe (1992) formula is an unsteady formula, which accounts for inertia effects. The sand transport rate is calculated by:
with and , empirical model coefficients. 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.
Timevarying velocity field
In the wave direction, the timevarying velocity is given by:
is defined here as the asymmetry parameter:
The wave cycle is decomposed into two parts:
1. During the first part of the wavecycle (), the velocities are directed onshore (). 2. During the second part (), the velocities are negatives ().
The sand transport rate in the wave direction is:
where and represent the amount of sand transported respectively in the onshore and offshore directions which will be deposited before flow reversal; et represent the remaining part which stays in suspension after flow eversal. The inertia of sediment depends on the ratio , and represented by parameter defined by:
with
1. all sediment is deposited before flow reversal
2. part of the sediment stays in suspension after flow reversal
The critical value of is calculated as a function of the wavecurrent nondimensional friction:
According to Temperville et Guiza (2000):
Keywords
 BEDLOAD TRANSPORT FORMULA: ICF= 9
 WAVE EFFECTS’ = YES
Sand grading effects
Sediment bed composition
Granulometry distribution
The number of classes of bed material is specified in the steering file and limited to . For uniform sediment, the granulometry distribution is represented by one single value () 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 () is defined by its mean diameter and volume fraction in the mixture at every node, noted . The characteristics of each class of sediment (mean diameter, Shields number , settling velocities 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:
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 (, noted ) is calculated based :
Bed structure
The bed is stratified into a number of layers () that enables vertical variations of the sediment bed composition. The percentage of each class of material (noted as well as the thickness of each layer are defined at each point and for each layer .
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 , composition of each layer are specified in usersubroutine 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):
In general, the thickness of the bed is very large (by default, m), so the bottom layer thickness need to be increased.
Assuming the porosity of each class to be identical (parameter is a constant), the total mass of sediments per unit area is:
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 :
Keywords
Initial sediment bed composition is defined by:
 NUMBER OF SIZE CLASSES OF BED MATERIAL (nsicla = 1, by default)
 SEDIMENT DIAMETERS (01;.01; . . . by default)
 SETTLING VELOCITIES (if not specified calculated by SISYPHE)
 SHIELDS PARAMETERS
 INITIAL FRACTION FOR PARTICULAR SIZE CLASS (ava0 = 1; 0. ;0. . . . , by default)
For more than one sizeclasses, the user should specify NSICLA values for the mean diameter, initial fraction. . . For example, for a mixture of two sands
 NUMBER OF SIZECLASSES OF BED MATERIAL = 2
 INITIAL FRACTION FOR PARTICULAR SIZE CLASS =.5;.5
 MEAN DIAMETER OF THE SEDIMENT =.02;.01
 SHIELDS PARAMETERS =.045;.01
 SETTLING VELOCITIES =.10;.05
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 (CVSM) 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:
where is the percentage of class 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 ().
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 ) 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 timeand spatially varying active layer thickness during the simulation. In Sisyphe, we simply assume ( 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 thinsized 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 CVSM overcomes many limitations of the classic layer implementation of Hirano (HRVSM). 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 CVSM output can be found in the new Selafin files VSPRES & VSPHYD in the tmpfolders. As the higher resolution of the CVSM 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
 VERTICAL GRAIN SORTING MODEL = 1
 0 = Layer = HRVSM (HIRANO + RIBBERINK as until SISYPHE v6p1)
 1 = CVSM
 CVSM MAXIMUM SECTIONS = 100
 Should be at least 4 + 4x Number of fractions,
 better > 100, tested up to 10000
 CVSM FULL PRINTOUT PERIOD = 1000
 0 ⇒ GRAPHIC PRINTOUT PERIOD
 Anything greater 0 ⇒ Sets an own printout period for the CVSP
 useful to save disk space!!!
 CVSM PRINTOUT SELECTION = 0251351411181750210433161212128021863187135630271535485
 Add any 2D Mesh Point numbers for .CSVAscii output of the CVSP
 Add 0 for full CVSP output as Selafin3D files (called VSPRES + VSPHYD)
 All files are saved to your working folder and in /VSP & /LAY folders below
 CVSM DYNAMIC ALT MODEL = 5
 Model for dynamic active layer thickness approximation
 0 = constant (Uses Keyword: ACTIVE LAYER THICKNESS)
 1 = Hunziker & Guenther
 2 = Fredsoe & Deigaard (1992)
 3 = van RIJN (1993)
 4 = Wong (2006)
 5 = Malcharek (2003)
 6 = 3*d50 within last time steps ALT'
General formulation
Bedload 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 bedload.
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).
Hidingexposure
In order to calculate the bedload transport rate for each class of sediment, it is necessary to take into account the effect of hidingexposure: 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 “hidingexposure” correction factor (HIDFAC, as defined below in §6.4.
Keywords
Sand grading effects are defined by:
 HIDING FACTOR FORMULA (Hidfac=0, by default)
 ACTIVE LAYER THICKNESS (Elay0= , by default)
 CONSTANT ACTIVE LAYER THICKNESS (= YES, by default)
 NUMBER OF BED LOAD MODEL LAYERS (NOMBAY= 2, by default)
 The number of layers has to be less than 10. The user has to compare the bed layers thickness in comparison to the total bed thickness (given in noerod.f) to ensure that there will be no more than 10 layers during the computation. The composition of the sediment AVAIL is defined at each point and for each layer: AVAIL(NPOIN,10,NSICLA). The sum of AVAIL on the classes must be equal to one everywhere.
 The number of layers is allowed to vary spatially and with time as the sediment bed undergoes erosion (empty sediment layers are suppressed) or deposition (sediment layers can be added).
 A new procedure will be developed in release 6.1 with a constant number of layers which will be compatible with the consolidation algorithm.
Sediment transport formulae
The formula of Hunziker ([24], 1995) is an adaptation of the MeyerPeter Mu”ller formula for fractional transport. The volumetric sediment transport per class is given by:
with the fraction of class in the active layer, the Shields parameter, the corrected critical Shields parameter, the critical Shields parameter (), the hiding / exposure factor, the grain size of class , the mean grain size of the surface layer (m), and 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 . It should be noted that according to Hunziker, stability problems may occur outside the parameter range:
Keywords
Sand grading effects are defined by:
 BEDLOAD TRANSPORT FORMULA (ICF = 6)
 HIDING FACTOR FORMULA (not used if ICF = 6)
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 bedload 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 MeyerPeter and Müller formula. The third formula (Karim et al)can be associated with any bedload predictor. The first two formulae modify the critical Shields parameter that will be used in the MeyerPeter formula, whereas the third formula modify directly the bedload transport rate.
 Formulation of Egiazaroff (HIDFAC=1)
 Formulation of Ashida & Michiue (HIDFAC=2)
 Formulation of Karim et al. (HIDFAC=4)