Sisyphe User Manual En V6p3
sisyphe_user_manual_en_v6p3
User Manual: Pdf
Open the PDF directly: View PDF 
.
Page Count: 76
| Download | |
| Open PDF In Browser | View PDF |
EDF R&D
NATIONAL HYDRAULICS AND ENVIRONMENT LABORATORY
NUMERICAL AND PHYSICAL MODELLING IN RIVER AND COASTAL HYDRODYNAMICS
6 quai Watier - 78401 CHATOU CEDEX, +33 (1) 30 87 79 46
January 13 2014
Front Page
Sisyphe v6.3 User's Manual
Pablo TASSI
Catherine VILLARET
H-P74-2012-02004-EN
1.0
SISYPHE is the state of the art sediment transport and bed evolution module of the TELEMACMASCARET modelling system.
SISYPHE can be used to model complex morphodynamics processes in diverse environments, such
as coastal, rivers, lakes and estuaries, for different flow rates, sediment size classes and sediment
transport modes.
In SISYPHE, sediment transport processes are grouped as bed-load, suspended-load or total-load,
with an extensive library of bed-load transport relations.
SISYPHE is applicable to non-cohesive sediments that can be uniform (single-sized) or non-uniform
(multiple-sized), cohesive sediments (multi-layer consolidation models), as well as sand-mud mixtures.
A number of physically-based processes are incorporated into SISYPHE, such as the influence of
secondary currents to precisely capture the complex flow field induced by channel curvature, the effect
of bed slope associated with the influence of gravity, bed roughness predictors, and areas of
inerodible bed, among others.
For currents only, SISYPHE can be tighly coupled to the depth-averaged shallow water module
TELEMAC-2D or to the three-dimensional Reynolds-averaged Navier-Stokes module TELEMAC-3D.
In order to account for the effect of waves or combined waves and currents, SISYPHE can be
internally coupled to the waves module TOMAWAC.
SISYPHE can be easily expanded and customized to particular requirements by modifying friendly,
easy to read fortran files. To help the community of users and developers, SISYPHE includes a large
number of examples, verification and validation tests for a range of applications.
Accessibility: Internal : EDF SA
Front page
Special mention:
Page I of III
Declassification:
EDF SA 2014
EDF R&D
H-P74-2012-02004-EN
Sisyphe v6.3 User's Manual
Version 1.0
Validation workflow
Author
Pablo TASSI
12/09/13
Reviewer
Nicole GOUTAL
12/19/13
Approval
Charles BODEL
01/05/14
Pre-diffusion
Recipient
Business code
Accessibility: Internal : EDF SA
Page II of III
EDF SA 2014
EDF R&D
H-P74-2012-02004-EN
Sisyphe v6.3 User's Manual
Version 1.0
Diffusion list
Group
P7-LNHE Chefs
P7-LNHE Chefs groupe
P73-MAHTRHYS
P74-SIMPHY
Recipient
Accessibility: Internal : EDF SA
Entity / Structure
Page III of III
Diffusion
EDF SA 2014
EDF R&D
Reference manual
Sisyphe 6.3
H-P74-2012-02004-EN
Version 1.0
AVERTISSEMENT / CAUTION
Optionnal, EDF publications only
L’accès à ce document, ainsi que son utilisation, sont strictement limités aux personnes expressément
habilitées par EDF.
EDF ne pourra être tenu responsable, au titre d’une action en responsabilité contractuelle, en responsabilité délictuelle ou de tout autre action, de tout dommage direct ou indirect, ou de quelque
nature qu’il soit, ou de tout préjudice, notamment, de nature financier ou commercial, résultant de
l’utilisation d’une quelconque information contenue dans ce document.
Les données et informations contenues dans ce document sont fournies ”en l’état” sans aucune
garantie expresse ou tacite de quelque nature que ce soit.
Toute modification, reproduction, extraction d’éléments, réutilisation de tout ou partie de ce
document sans autorisation préalable écrite d’EDF ainsi que toute diffusion externe à EDF du présent
document ou des informations qu’il contient est strictement interdite sous peine de sanctions.
-------
The access to this document and its use are strictly limited to the persons expressly authorized to
do so by EDF.
EDF shall not be deemed liable as a consequence of any action, for any direct or indirect damage, including, among others, commercial or financial loss arising from the use of any information
contained in this document.
This document and the information contained therein are provided ”as are” without any warranty of any kind, either expressed or implied.
Any total or partial modification, reproduction, new use, distribution or extraction of elements
of this document or its content, without the express and prior written consent of EDF is strictly
forbidden. Failure to comply to the above provisions will expose to sanctions.
Accessibility : Interne : EDF SA
Page 1 sur 73
c EDF SA 2013
EDF R&D
Reference manual
Sisyphe 6.3
H-P74-2012-02004-EN
Version 1.0
Optionnal, EDF publications only
Executive Summary
The present document is part of the openTELEMAC documentation (http://www.opentelemac.
org).
Accessibility : Interne : EDF SA
Page 2 sur 73
c EDF SA 2013
EDF R&D
Reference manual
Sisyphe 6.3
H-P74-2012-02004-EN
Version 1.0
Table des matières
1
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6
2
Running a sedimentological computation . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
8
2.1
Input files . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
8
2.2
Steering file . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
8
2.3
Coupling hydrodynamics and morphodynamics . . . . . . . . . . . . . . . . . . . . . . .
8
2.4
Boundary conditions file . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.5
Fortran files . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.6
Parallel computing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
3
Hydrodynamics parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
3.1
Current induced bed shear stresses . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
3.2
Bed roughness . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
3.3
Skin friction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
4
Sediment parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
4.1
Sediment properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
4.2
Shields parameter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
4.3
Settling velocities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
4.4
Transport rate . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
5
Bedload transport . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
5.1
Bedload transport formulae . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
5.2
Sediment transport formulae . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
5.3
Bed slope effect . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
5.4
Sediment Slide . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
5.5
Bedload transport in curved channels . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
6
Suspended load . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
6.1
Main assumptions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
6.2
Two-dimensional sediment transport equation . . . . . . . . . . . . . . . . . . . . . . . . . . 26
6.3
Convection velocity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
6.4
Numerical treatments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
6.5
Equilibrium concentrations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
6.6
Initial and boundary conditions for sediment concentrations . . . . . . . . . . . . . . . . . 30
7
Bed evolution equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
7.1
Bedload . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
7.2
Boundary conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
7.3
Treatment of rigid beds and tidal flats . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
7.4
Tidal flats . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
Accessibility : Interne : EDF SA
Page 3 sur 73
c EDF SA 2013
EDF R&D
7.5
8
Reference manual
Sisyphe 6.3
H-P74-2012-02004-EN
Version 1.0
Bed evolution for suspension . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
Wave effects . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
8.1
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
8.2
Wave-induced bottom friction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
8.3
Wave-induced ripples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
8.4
Wave-induced sediment transport . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
8.5
Internal coupling waves-currents and sediment transport . . . . . . . . . . . . . . . . . . . 40
9
Sand grading effects . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
9.1
Sediment bed composition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
9.2
Sediment transport of sediment mixtures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
10 Cohesive sediment properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
10.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
10.2 Effect of floculation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
10.3 Consolidation process . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
11 Cohesive bed structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
11.1 Multi-layer model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
11.2 Initialization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
12 Erosion and deposition rates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
12.1 Erosion flux . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
12.2 Deposition flux . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
13 Bed Evolution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
14 Mass conservation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
15 Consolidation algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
15.1 Multi-layer empirical algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
15.2 Multi-layer iso-pycnal Gibson’s model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
15.3 Vertical grid Gibson’s model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
16 Closure equations for permeability and effective stress . . . . . . . . . . . . . . . . . . . . . 59
17 Mixed sediments processes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
17.1 Mixte sediment bed composition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
17.2 Erosion/deposition of sand-mud sediment bed . . . . . . . . . . . . . . . . . . . . . . . . 63
17.3 Critical erosion shear stress of sand mud mixture . . . . . . . . . . . . . . . . . . . . 63
17.4 Iterative procedure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
17.5 Definition scheme of sand-mud mixture in Sisyphe . . . . . . . . . . . . . . . . . . . . . 65
17.6 Suspended sediment transport equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
17.7 Fortran subroutines for mixed sediment transport . . . . . . . . . . . . . . . . . . . . . . . . 68
Accessibility : Interne : EDF SA
Page 4 sur 73
c EDF SA 2013
EDF R&D
Reference manual
Sisyphe 6.3
H-P74-2012-02004-EN
Version 1.0
Conventions
In this document, the following conventions are used :
– The names of the different modules are written in caps and small caps
– File names are written in sans serif font
– Keywords, variables, subroutines, etc. are written in monospaced ‘‘typewriter’’ font
Accessibility : Interne : EDF SA
Page 5 sur 73
c EDF SA 2013
EDF R&D
1
Reference manual
Sisyphe 6.3
H-P74-2012-02004-EN
Version 1.0
Introduction
Sisyphe is a sediment transport and morphodynamic model which is part of the hydroinformatics finite element and finite volume system Telemac-Mascaret. In Sisyphe , sediment transport
rates, split into bedload and suspended load, are calculated at each node 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 a classical sediment transport formula
from the literature. The suspended load is determined by solving an additional transport equation for the depth-averaged suspended sediment concentration. The bed evolution equation (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 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 and slope failure. For
cohesive sediments, the effect of bed consolidation can be accounted for.
Sisyphe can be applied to a large variety of hydrodynamic flow conditions including rivers,
estuaries and coastal applications. For the later, 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, Chézy or user defined) or
by a bed-roughness predictor.
In Sisyphe , the relevant hydrodynamic variables can be either imposed in the model (chaining
method) or calculated by a hydrodynamic computation (internal coupling) by using one of the
hydrodynamic modules of the Telemac system (modules Telemac-2d , Telemac-3d or Tomawac )
or an external hydrodynamic model.
Sisyphe can be run on Unix, Linux or Windows. The latest release of Sisyphe (release 6.3) uses
the version 6.3 of the Bief finite element library (Telemac system library). The reader can refer to
the specific documentation for pre- and post-processing tools, e.g. [10].
Outline of the manual
This document is organized in three parts, as follows. In the first part (Part I), the main steps
to run a sedimentological computation are given in Section 2. In Sections 3 and 4, a description of
the main sediment and hydrodynamics parameters such that total bed shear stress, skin friction,
etc. is presented. In the second part (Part II), specific aspects of non-cohesive sediment transport
modelling are give. In Section 5 the bedload is presented. In Section 6, the suspended load transport
is introduced. Details on bed evolution for bedload and suspended-load are given in Section 7.
Section 8 presents the influence of the effect of waves on the sediment transport. In Section 9, sand
grading effects are introduced. Finally, in Part III numerical aspects of cohesive and mixed (sandmud) sediment transport modelling are introduced.
Accessibility : Interne : EDF SA
Page 6 sur 73
c EDF SA 2013
EDF R&D
Reference manual
Sisyphe 6.3
H-P74-2012-02004-EN
Version 1.0
Part I
Basis
by Catherine Villaret & Pablo Tassi
Accessibility : Interne : EDF SA
Page 7 sur 73
c EDF SA 2013
EDF R&D
2
Reference manual
Sisyphe 6.3
H-P74-2012-02004-EN
Version 1.0
Running a sedimentological computation
2.1
Input files
The minimum set of input files to run a Sisyphe simulation includes the steering file (format ascii
file cas), the geometry file (format selafin slf), and the boundary conditions file (format ascii file cli).
Additional or optional input files include the fortran file, the reference file, the results file, etc.
2.2
Steering file
The steering file contains the necessary information for running a computation, including the
physical and numerical parameters that are different from the default values. Like any other modules
of the Telemac system, the model parameters can be specified in the (obligatory) Sisyphe steering
file. The following essential information (input/output) must 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.)
2.3
Coupling hydrodynamics and morphodynamics
We describe here two methods for linking the hydrodynamic and the morphodynamic models :
by chaining (the flow is obtained from a previous hydrodynamic simulation assuming a fixed bed)
or by internal coupling (both the flow and bed evolution are updated at each time step).
2.3.1
Chaining method
– Principle
For this method, both models (hydrodynamic and morphodynamic) are running independently. During the first hydrodynamic simulation, the bed is assumed to be fixed. Then, in
the subsequent morphodynamic step, the flow rate and free surface are read from the previous
hydrodynamic results file. The 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 unsteady mode, that means that the flow field is
linearly interpolated between two time steps of the hydrodynamic file. For steady flows, the
last time step of the hydrodynamic file is used and the flow rate and free surface is assumed
to be constant while the bed evolves.
– Flow updating
At each time step, the flow velocity is updated by assuming that both the flow rate and free
surface elevation are conserved. In the case of deposition, the flow velocity is locally increased,
whereas in the case of erosion, the flow velocity decreases.
This schematic updating does not take into account any deviation of the flow. It is only suitable
for simple flows and assuming relatively small bed evolutions. The morphodynamic computation is stopped when the bed evolution reaches a certain percent of the initial water depth.
This methodology is not 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
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. Furthermore, it can also be
Accessibility : Interne : EDF SA
Page 8 sur 73
c EDF SA 2013
EDF R&D
Reference manual
Sisyphe 6.3
H-P74-2012-02004-EN
Version 1.0
responsible of numerical instabilities.
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 file. Hydrodynamics records are interpolated at each time step of the sedimentological computation.
Important
An error may occur when the TIDE PERIOD is not a multiple of the graphic printout
period 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 the sedimentological computation.
– Steering and fortran files
For uncoupled mode, the Sisyphe steering file should specify :
– The time steps, graphical or listing output, final time
– The hydrodynamic file as yielded by Telemac-2d or Telemac-3d (HYDRODYNAMIC FILE) or
by the 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
– Friction options : Selects the type of formulation used for the bottom friction.
Keywords
For time step, final time 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
– LAW OF BOTTOM FRICTION = 3, default value
– FRICTION COEFFICIENT = 50, default value
For waves :
– WAVE FILE, WAVE EFFECTS
Important
In the case of internal coupling, the friction coefficient is selected in the Telemac-2d or
Telemac-3d steering file except when BED ROUGHNESS PREDICTION is set to YES.
2.3.2
Internal coupling
– Principle
Sisyphe can be internally coupled to the hydrodynamic model, Telemac-2d or Telemac-3d .Sisyphe is internally coupled with 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,
Accessibility : Interne : EDF SA
Page 9 sur 73
c EDF SA 2013
EDF R&D
Reference manual
Sisyphe 6.3
H-P74-2012-02004-EN
Version 1.0
instead of being written and read in a file.
At each time step, the hydrodynamics variables (velocity field, water depth and bed shear
stress) are transferred into the morphodynamic model, which sends back the updated bed elevation to the hydrodynamic model.
– Time step and coupling period
The internal coupling method is more time consuming than the chaining method. Various
techniques can be set up to reduce the CPU time (e.g. parallel processors). 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 suspended load, a convection-diffusion 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 Telemac-3d multiplied by
the COUPLING PERIOD. The graphic and listing printout periods are the same as in the Telemac
computation.
The Telemac-2d and Telemac-3d steering files 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 for example there is no Fortran file for
Telemac-2d or Telemac-3d steering file).
Keywords
For internal coupling, the following keywords need to be specified in the Telemac-2d or
Telemac-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-2d or Telemac-3d steering file, but are no longer used by Sisyphe .
The values of time step, bottom shear stress, etc. are transferred directly to Sisyphe .
For long term simulations, it is possible to use the keyword MORPHOLOGICAL FACTOR (in order
to speed-up the computational time. This method allows to increase the time step used for the
morphological simulations. Further details can be found in [40].
Keywords
For internal coupling, the following keywords need to be specified in the Sisyphe steering
file :
– MORPHOLOGICAL FACTOR =1, default value
Obsolete keywords
For the new version of Sisyphe , the following keywords are obsolete :
Important
The following keywords are no longer in use is (Sisyphe ) steering file :
– TIME STEP
– GRAPHIC PRINTOUT PERIOD
– LISTING PRINTOUT PERIOD
– LAW OF BOTTOM FRICTION
– FRICTION COEFFICIENT
– VARIABLE TIME STEP (DTVAR)
– COMPUTATION METHOD (lengthening of the tide, remplaced by MORPHOLOGICAL FACTOR
– FILTERING COEFFICIENT
– NON EQUILIBRIUM BED LOAD (NOEQBED)
– GRAIN FEEDING (LGRAFED)
Accessibility : Interne : EDF SA
Page 10 sur 73
c EDF SA 2013
EDF R&D
2.4
H-P74-2012-02004-EN
Reference manual
Sisyphe 6.3
Version 1.0
Boundary conditions file
The format of the boundary condition file is the same as for Telemac-2d or Telemac-3d . This
ascii file can be created by a mesh generator (for example, BlueKenue [10]) and modified using a
text editor. Each line is related to a point along the edge of the mesh. 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 edge points are numbered like the file lines ; the numbering first describes the domain outline
in the counterclockwise direction from lower left point (X+Y minimum), then the islands in the
clockwise direction.
For Telemac-2d or Telemac-3d , the information contained in the hydrodynamics boundary
condition file includes :
(1)
LIHBOR
(2)
LIUBOR
(3)
LIVBOR
(4)
QBOR
...
...
(8)
LITBOR
(9)
TBOR
Further details for the boundary condition file for Sisyphe are given in Chapter 5 (bedload) and
Chapter 6 (suspended load).
2.5
Fortran files
Programming can be necessary for particular applications by modifying a user’s subroutine from
the sources. When a subroutine(s) is necessary for a particular case, a Fortran file (keyword FORTRAN
FILE) needs to be specified in the Telemac-2d or Telemac-3d steering file. Some typically used
subroutines are briefly described below :
– Coordinates modification : corrxy.f was designed for carry-out, at the beginning of a computation, a modification of the nodes coordinates, e.g. for a scaling-up (shifting from a small-scale
model to a full size), a rotation or a translation. This subroutine is empty by default and provides, in the form of a comment, an example of programmed change of scale and origin.
– Time variable or space variable friction : strche.f is used for defining a space-variable or
time variable friction coefficient. By default, the value given for the FRICTION COEFFICIENT
keyword in the steering file is imposed throughout the domain. The corstr subroutine is for
imposing a time-variable friction coefficient, which can be computed as a function of liquid
flow, water depth, etc.
– Definition of rigid areas : noerod.f is used for specifying the rigid areas. The depths of nonerodable areas (array ZR) are imposed in this subroutine.
– New sand transport formula : qsform.f can be used to program a sediment transport formula that is different from those already implemented in Sisyphe . The keyword BED-LOAD
TRANSPORT FORMULA should be set to 0.
– Declaration of private variables : subroutines nomvar.f and predes.f. In the standard version
of Sisyphe , it is possible to output a number of computed variables. In some cases, the user
may want to compute variables and write them into the results file (so far, the number of
these variables is limited to four). Likewise, each variable is identified by one or several letters
in the keyword VARIABLES FOR GRAPHIC PRINTOUTS. The new variables are then identified by
letters A, G, L and O which, respectively, correspond to numbers 23, 24, 25 and 26. In the data
structure, these four variables correspond to the PRIVE%ADR(1)%P%R(X), PRIVE%ADR(2)%P%R(X),
PRIVE%ADR(3)%P%R(X) and PRIVE%ADR(4)%P%R(X) arrays (where X is the number of grid points)
that can be used in several programming stages. The new variables are programmed in two
steps :
– Firstly, the names of these new variables should be set by completing the nomvar.f subroutine, which consists of two equivalent structures for English and French languages. Each
structure sets the names of variables of the results file to be generated, then the names of
variables to be read again in the previous computation file in case of a reactivation. This su-
Accessibility : Interne : EDF SA
Page 11 sur 73
c EDF SA 2013
EDF R&D
Reference manual
Sisyphe 6.3
H-P74-2012-02004-EN
Version 1.0
broutine can also be amended when, for instance, the French version of Sisyphe being used,
one wants to make a reactivation on a file generated by the English version. If so, the TEXTPR
array of the subroutine French part shall contain the English names of the variables.
– Secondly, the PREDES subroutine is to be amended for entering the computation(s) of the
new variable(s) thereinto. Besides, the LEO and IMP variables are provided for determining
whether, at the time step being considered, the variable should be printed in the monitoring
listing or in the results file.
Important
In the case of internal coupling, the Sisyphe fortran file can be specified in the Telemac-2d
or Telemac-3d steering file instead. This is necessary if the fortran file is not present for the
hydrodynamics.
2.6
Parallel computing
For running simulations in parallel mode, the user must specify the keyword PARALLEL PROCESSORS
in the Telemac-2d or Telemac-3d steering file (for internal coupling simulations).
Accessibility : Interne : EDF SA
Page 12 sur 73
c EDF SA 2013
EDF R&D
3
Reference manual
Sisyphe 6.3
H-P74-2012-02004-EN
Version 1.0
Hydrodynamics parameters
3.1
3.1.1
Current induced bed shear stresses
Coupled mode
The bed shear stress is one of the most important hydrodynamics parameters regarding sediment
transport applications. The current-generated bed shear stress is used in 2D in the shallow water
momentum equation and in 3D as the bottom boundary condition for the velocity profile. When
Sisyphe is internally coupled with Telemac-2d , the bed shear stress term τ0 = (τx , τy ) is calculated
at each time step from the classical quadratic dependency on the depth–averaged velocity :
τx,y =
1
ρC (u, v)|u|,
2 d
(1)
where (u, v) are the depth-averaged√velocity components along the x − and y− Cartesian directions, respectively ; |u| = |(u, v)| = u2 + v2 the velocity module ; and friction coefficient Cd . The
magnitude of the bed shear stress τ0 can be related to the friction velocity u∗ , defined by :
τ0 = ρu2∗ .
(2)
Sisyphe coupled with Telemac-2d
When the model is coupled with Telemac-2d , the flow variables and the values of the friction
coefficients (and therefore, the bed shear stresses) are provided by Telemac-2d . The 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 (bed slope effect), see § 5.3
– secondary currents, see § 5.5
Sisyphe coupled with Telemac-3d
When Sisyphe is coupled with Telemac-3d , the bed shear stress is aligned with the near bed
velocity in order to account for possible flow deviations. The magnitude of the bed shear stress
is still related to the depth-averaged velocity, except if the Nikuradse friction law is applied.
In this case, the friction velocity is related to the near bed flow velocity u(z1 ) by assuming a
logarithmic velocity profile near the bed :
u∗
z1
u ( z1 ) =
ln
(3)
κ
z0
where z0 is expressed as a function of the Nikuradse bed roughness (z0 = k s /30), with k s the
grain roughness height, z1 is the distance of the first vertical plane from the bed level ; and
κ = 0.4 is the von Kármán constant. For flat beds, the roughness height has been shown to be
approximately k s ≈ 3d50 [29], with d50 the grain size with 50% of the material finer by weight.
3.1.2
Uncoupled mode
The quadratic friction coefficient Cd , which is used to calculate the total bed shear stress, can
be calculated after the selected friction law. Different options, which are consistent with Telemac2d options are available in Sisyphe and depend on the choice of the keywords LAW FOR BOTTOM
FRICTION and on the value of the FRICTION COEFFICIENT. Similar keywords are available in order
to account for the effect of the friction for Telemac-2d , Telemac-3d and Sisyphe .
Accessibility : Interne : EDF SA
Page 13 sur 73
c EDF SA 2013
EDF R&D
H-P74-2012-02004-EN
Reference manual
Sisyphe 6.3
– Chézy coefficient Ch (KFROT = 2)
Version 1.0
2g
Ch2
(4)
Cd =
2g 1
St2 h1/3
(5)
Cd =
2g
M2
h1/3 a
(6)
Cd =
– Strickler coefficient St (KFROT = 3)
– Manning friction Ma (KFROT = 4)
– Nikuradse bed roughness k s (KFROT = 5)
"
Cd = 2
κ
log( 12h
ks )
#2
,
(7)
Keywords
In the Sisyphe steering file, for the uncoupled mode the total bed shear stress is calculated
based on :
– LAW OF BOTTOM FRICTION (KFROT =3, default option)
– BOTTOM FRICTION COEFFICIENT (St=50, default value)
Important
In the case of internal coupling with Telemac-2d or Telemac-3d , the friction keywords provided in the Sisyphe steering file are not longer used.
3.2
3.2.1
Bed roughness
Role of bed forms
A natural sediment bed is generally covered with bed forms, in general characterized by their
length λd and height ηd . The presence of bed forms modifies the boundary layer flow structure, with
the formation of recirculation cells and depressions in the lee of bed forms. Depending on the flow
and sediment transport rates, the size of bed forms ranges from a few centimeters for ripples to
tens of meters for dunes or bars. The dimension of dunes scales with the water depth h, such that
ηd ≈ 0.4h and λd ≈ 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 increading the
friction coefficient. To determine bed roughness, there are two options available in Sisyphe . The
simplest one is to impose the friction coefficient based on friction laws. As in most applications the
bed roughness is unknown, an alternative is to predict the value of the bed roughness as a function
of flow and sediment parameters using a bed roughness predictor.
Important
In coupled mode, the bed roughness predicted by Sisyphe is sent to the hydrodynamics model
in order to avoid inconsistencies between hydrodynamics and morphodynamics simulations.
3.2.2
Bed roughness predictor
Different options are programmed in Sisyphe to predict the total bed roughness through the associated keywords BED ROUGHNESS PREDICTION and BED ROUGHNESS PREDICTOR OPTION. It is recalled
Accessibility : Interne : EDF SA
Page 14 sur 73
c EDF SA 2013
EDF R&D
H-P74-2012-02004-EN
Reference manual
Sisyphe 6.3
Version 1.0
that the bed friction option of Sisyphe is not used in the case of internal coupling with Telemac-2d
or Telemac-3d .
If the keyword BED ROUGHNESS PREDICTION is activated (= YES), the following options are available via the keyword BED ROUGHNESS PREDICTOR OPTION :
– For IKS = 1 : the bed is assumed to be flat k s = k0s = αd50 , with α a constant (assumed to be
equal to 3). The α parameter can be changed according to the keyword RATIO BETWEEN SKIN
FRICTION AND MEAN DIAMETER.
– IKS = 2 : the bed is assumed to be covered by ripples.
– For currents only, the ripple bed roughness is function of the mobility number Ψ = U 2 /(s −
1) gd50 , see [54] :
d50 (85 − 65 tanh(0.015(Ψ − 150))) for Ψ < 250
kr =
20d50
otherwise
– For waves and combined waves and currents, bedform dimensions are calculated as a function of wave parameters following the method of Wiberg and Harris [60]. The wave-induced
bedform bed roughness kr is calculated as a function of the wave-induced bedform height
ηr :
kr = max(k0s , ηr ).
(8)
Then,
k s = k0s + kr .
(9)
– IKS = 3 : for currents only, the van Rijn’s total bed roughness predictor has been implemented [27, 54]. The total bed roughness can be decomposed into a grain roughness k0s , a small-scale
ripple roughness kr , a mega-ripple component k mr , and a dune roughness k d :
q
(10)
k s = k0s + k2r + k2mr + k2d .
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). In Equation (10), the general expression for megaripple roughness k mr is given by :
k mr = 0.00002 f ts h(1 − exp−0.05Ψ )(550 − Ψ),
with
f ts =
d50 /(1.5dsand )
1
for d50 < 1.5dsand
otherwise
and the general expression for dune roughness k d is estimed by :
k d = 0.00008 f ts h(1 − exp−0.02Ψ )(600 − Ψ).
Keywords
In Sisyphe steering file, the bed roughness prediction is calculated based on :
– BED ROUGHNESS PREDICTION, (= NO default option)
– BED ROUGHNESS PREDICTOR OPTION (IKS = 1 : Flat bed (default option), IKS = 2 : Rippled
bed, IKS = 3 : Dunes and mega ripples (van Rijn’s method))
– RATIO BETWEEN SKIN FRICTION AND MEAN DIAMETER = 3, default option. If = 0 use skin
friction prediction from Van Rijn (2007) for currents or the Wiberg and Harris method for
waves.
Accessibility : Interne : EDF SA
Page 15 sur 73
c EDF SA 2013
EDF R&D
3.3
3.3.1
H-P74-2012-02004-EN
Reference manual
Sisyphe 6.3
Version 1.0
Skin friction
Bed shear stress decomposition
In the presence of bedforms, the total bed shear stress can be expressed as the sum of two
components :
τ0 = τ00 + τ000
(11)
where τ0 is the total bed shear stress, τ00 is the grain (or skin) shear stress, and τ000 is the form shear
stress. The local skin friction component determines the bedload transport rate and the equilibrium
concentration for the suspension. The total friction velocity determines the turbulence eddy viscosity/diffusivity vertical distribution in 3D models, and therefore determines both the velocity vertical
profile and the mean concentration profile.
Bedload transport rates are calculated as a function of the local skin friction component τ00 . The
total bed shear stress issued from the hydrodynamics model needs to be corrected in the morphodynamics model as follows :
τ00 = µτ0 ,
(12)
where µ is a correction factor for skin friction. Physically, the skin bed roughness should be smaller
than the total bed roughness (i.e. µ ≤ 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 µ > 1 can be admitted.
3.3.2
Correction factor for skin friction
Different methods are programmed in Sisyphe in order to calculate the bedform correction factor
µ, according to the keyword SKIN FRICTION CORRECTION :
– 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 (k0s ∼ d50 ). The proportionality coefficient is specified by the keyword
RATIO BETWEEN SKIN FRICTION AND MEAN DIAMETER, defined as :
µ=
Cd0
,
Cd
(13)
where Cd et Cd0 are both quadratic friction coefficients related to total friction and skin friction,respectively. Cd is obtained from Telemac-2d or Telemac-3d and Cd0 is calculated from k0s ,
as follows :
"
#2
κ
0
Cd = 2
(14)
log( 12h
)
k0
s
– ICR = 2 the bedform predictor is used to calculate the bedform roughness kr in order to account for the effect of ripples. Both kr and k0s should influence the transport rates. It is assumed
that :
C 00.75 Cr0.25
µ= d
,
(15)
Cd
where the quadratic friction Cr due to bedforms is calculated as a function of kr .
Accessibility : Interne : EDF SA
Page 16 sur 73
c EDF SA 2013
EDF R&D
Reference manual
Sisyphe 6.3
H-P74-2012-02004-EN
Version 1.0
Keywords
The option selected for the skin friction correction is based on keywords :
– SKIN FRICTION CORRECTION (ICR=0, default option)
Accessibility : Interne : EDF SA
Page 17 sur 73
c EDF SA 2013
EDF R&D
4
Reference manual
Sisyphe 6.3
H-P74-2012-02004-EN
Version 1.0
Sediment parameters
4.1
Sediment properties
Fine sediment particles of grain size d50 < 60µm present complex cohesive properties which
affect the sediment transport processes. For non-cohesive sediments (median diameter d50 > 60µm),
the grain diameter and grain density ρs are the key parameters which determine its resistance to
erosion and sediment transport rate.
In Part II, we consider non-cohesive sediments either characterized by one single value for the
grain size d50 and grain density ρs or represented by a number of classes for sand grading.
In Part III, we deal with cohesive sediments (d50 < 60µm) : the grain diameter is no longer the key
sediment parameter and the settling velocity now depends on the flocculation state of suspension,
whereas the critical bed shear strength depends on the consolidation state of the sediment bed.
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 (d50 > 60µm for non-cohesive sediments)
– SEDIMENT DENSITY (ρs = 2650 kgm−3 , default value)
4.2
Shields parameter
The critical Shields number or dimensionless critical shear stress Θc is defined by :
Θc =
τc
g(ρs − ρ)d50
(16)
where τc is the critical shear stress for sediment incipient motion. Values of Θc 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 D∗ :
D∗ = d[(ρs /ρ − 1) g/ν2 ]1/3 .
(17)
The critical Shields number is implemented in the subroutine init sediment.f as follows :
0.24D∗−1 ,
0.14D∗−0.64 ,
Θc =
0.04D∗−0.10 ,
0.013D∗0.29 ,
0.045,
D∗ ≤ 4
4 < D∗ ≤ 10
10 < D∗ ≤ 20
20 < D∗ ≤ 150
150 ≤ D∗
where τc and d50 are in N m−2 and in m, respectively. If the value of the Shields number is not
specified in the steering file, Sisyphe computes the value according the Equation (4.2) as a function
of the non-dimensional grain diameter D∗ . For multi-grain sizes, the Shields parameter needs to be
specified for each class.
Keywords
– SHIELDS PARAMETERS
Accessibility : Interne : EDF SA
Page 18 sur 73
c EDF SA 2013
EDF R&D
H-P74-2012-02004-EN
Reference manual
Sisyphe 6.3
Version 1.0
Attention
For multiple grain size (NSICLA > 1), Θc is an array which allows to use different values for
each class of material.
4.3
Settling velocities
The settling velocity ws is an important parameter for suspended sediment transport. It can
be either specified or calculated by the model as a function of the grain diameter. The van Rijn
formula [52, 53] which is valid for non-cohesive spherical particles and dilute suspensions, has been
implemented in Sisyphe as follows :
(s − 1) gd250
,
18ν
s
3
(
s
−
1
)
gd
10ν
50
ws =
1 + 0.01
− 1 ,
18ν2
d50
q
1.1 (s − 1) gd ,
50
if d50 ≤ 10−4
if 10−4 ≤ d50 ≤ 10−3
otherwise
with s = ρs /ρ0 is the relative density, ν is the fluid viscosity and g is gravity.
Keywords
– SETTLING VELOCITIES
Attention
The default value is not given. If the user does not give a value, the subroutine vitchu.f
computes their value as a function of the grain size, the relative density and the fluid viscosity.
4.4
Transport rate
When the current-induced bed shear stress exceeds the critical threshold value coarse sediment
particles start to move as bedload, while the finer particles are transported in suspension. The total
sediment load Qt includes a bedload Qb and suspended load Qs components :
Qt = Qb + Qs .
(18)
Bedload occurs in a very thin high concentrated near-bed layer, where inter-particle interactions
develop. In Sisyphe , we assume the classical delimitation between bedload and suspended load.
The interface between the bedload and suspended load is located at z = Zre f :
– In the thin high-concentrated bedload layer (z < Zre f ) inter-particle interactions and flowturbulent interactions strongly modify the flow structure. Equilibrium conditions are however
a reliable assumption to relate the bed-load to the current induced bed shear stress.
– In the upper part of the flow (z > Zre f ), 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.
Accessibility : Interne : EDF SA
Page 19 sur 73
c EDF SA 2013
EDF R&D
Reference manual
Sisyphe 6.3
H-P74-2012-02004-EN
Version 1.0
Part II
Non-cohesive sediments
by Catherine Villaret & Pablo Tassi
Accessibility : Interne : EDF SA
Page 20 sur 73
c EDF SA 2013
EDF R&D
5
H-P74-2012-02004-EN
Reference manual
Sisyphe 6.3
Version 1.0
Bedload transport
5.1
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 among 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 and Mü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 Φs , is expressed as :
Φs = q
Qb
(19)
g(s − 1)d3ch
with ρs the sediment density ; s = ρs /ρ the relative density ; d the characteristic sand grain diameter
(= dch for uniform grains) ; and g the gravity. The characteristic sand grain diameter can be chosen as
d50 initially. As presented next, the non-dimensional sand transport rate Φs is, in general, expressed
as a function of the non-dimensional skin friction or Shields parameter θ 0 , defined by :
θ0 =
µτ0
,
(ρs − ρ) gdch
(20)
with the correction factor for skin friction µ and the bottom shear stress τ0 .
5.2
Sediment transport 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 BED-LOAD TRANSPORT FORMULA. By setting ICF = 0, the
user can program a specific transport formula through the subroutine qsform.f.
Keywords
– BED-LOAD TRANSPORT FORMULA (default option : Meyer-Peter-Müller formula ICF=1)
– ICF = 1 (Meyer-Peter-Müller formula) : this classical bed-load formula has been validated for
coarse sediments in the range (0.4 mm < d50 < 29 mm). It is based on the concept of initial
entrainment :
0
if θ 0 < θc
Φb =
0
αmpm (θ − θc )3/2 otherwise
with αmpm a coefficient (= 8 by default), θc the critical Shields parameter (= 0.047 by default).
The coefficient αmpm can be modified in the steering file by the keyword MPM COEFFICIENT (=
8.0 by default). The characteristic sand grain diameter was originally chosen to d64 by MeyerPeter and Müller.
– ICF = 2 (Einstein-Brown formula) : this bed-load formula is recommended for gravel (d50 > 2
mm) and large bed shear stress θ > θc . The solid transport rate (see [17]) is expressed as :
Φ b = F ( D∗ ) f ( θ 0 ) ,
with
F ( D∗ ) =
Accessibility : Interne : EDF SA
2
36
+
3 D∗
0.5
Page 21 sur 73
−
(21)
36
D∗
0.5
,
(22)
c EDF SA 2013
EDF R&D
H-P74-2012-02004-EN
Reference manual
Sisyphe 6.3
and
f (θ 0 ) =
2.15 exp(−0.391/θ 0 )
0
40 θ 3
Version 1.0
if θ 0 ≤ 0.2
otherwise
where the non-dimensional diameter D∗ is defined according to Equation (17).
– ICF = 3 or 30 (Engelund-Hansen formula) : this formula predicts the total load (bedload +
suspended load). It is recommended for fine sediments, in the range 0.2 mm < d50 < 1 mm but
beware that the use of a total load formula is only suitable under equilibrium conditions (quasisteady and uniform flow). The 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 :
Φs = 0.1 θ 05/2
(23)
– ICF = 3 corresponds to the version modified by Chollet and Cunge [12] to account for the
effects of sand dunes. The transport rate is related to the total bed shear stress as :
Φs =
where the dimensionless bed shear
skin friction θ 0 :
0p
2.5(θ 0 − 0.06)
θ̂ =
1.065θ 00.176
0
θ
0.1 5/2
θ̂ ,
Cd
(24)
stress θ̂ is calculated as a function of the dimensionless
if θ 0 < 0.06 (flat bed regime - no transport)
if 0.06 < θ 0 < 0.384 (dune regime)
if 0.384 < θ 0 < 1.08 (transition regime)
if θ 0 > 1.08 (flat bed regime - upper regime)
Attention
To avoid the suspended load to be calculated twice, in the case of coupling between
bedload and suspended load (SUSPENSION = YES), the total load formula (ICF = 3 or
30) should not be used.
– ICF = 7 (van Rijn’s formula) : this formula was proposed by van Rijn [51] to calculate the
bedload transport rate for particles of size 0.2 mm < d50 < 2 mm :
θ p − θcr 2.1
Φb = 0.053D∗−0.3
.
(25)
θcr
Attention
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). The Bijker’s formula (ICF = 4) is a total load formula (bedload
+ suspended load) that can also be used to account for currents only.
5.2.1
Validity range of sediment transport formulae
Most sediment transport formulae are based on experiments performed under fluvial, unidirectional flows. These formulae shown a rapid variation of the bedload transport prediction, as a
function of the mean flow intensity. Therefore, an increasing of the current velocity by 10% will
result, depending on the formula being used, in an increasing of the transport rate of 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. tidal regime), the average transport will be highly influenced by the stronger currents and will not be directly related to the mean
flow.
Accessibility : Interne : EDF SA
Page 22 sur 73
c EDF SA 2013
EDF R&D
H-P74-2012-02004-EN
Reference manual
Sisyphe 6.3
Version 1.0
The validity range of the different formulae is briefly summarized in Table 1.
Formula
IFC
Mode of transport
Validity range (d50 )
Meyer-Peter& Müller
1
bedload
0.4 − 29 mm
Einstein-Brown
2
bedload
0.25 − 32 mm
Engelund-Hansen
3 or 30
total load
0.19 − 0.93 mm
van Rijn
7
bedload
0.20 − 2.0 mm
Table 1 – Validity range of some of the sand transport formulae programmed in Sisyphe (for currents
only).
5.3
Bed slope 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 bedload direction. In Sisyphe , a correction factor can be applied
to both the magnitude and direction of the solid transport rate, before solving the bed evolution
equation. The bed slope effect is activated if the keyword SLOPE EFFECT is present in the Sisyphe
steering file.
Two different formulations for both effects are available depending on the choice of the keywords
FORMULA FOR SLOPE EFFECT, which chooses the magnitude correction and FORMULA FOR DEVIATION,
which chooses the direction correction. At the open boundaries the slope effect is disarmed due to
stability reasons.
5.3.1
Correction of the magnitude of bedload transport rate (formula for slope
effect)
– SLOPEFF = 1 : this correction method is based on the Koch and Flokstra’s formula [31]. The
solid transport rate intensity Qb0 is multiplied by a factor 1 − β∂Z f /∂s, thus :
∂Z f
,
(26)
Qb = Qb0 1 − β
∂s
with s the coordinate in the current direction and β an empirical factor, which can be specified
with the associated 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 smooth the results and is
often used to reduce unstabilities.
– SLOPEFF=2 : this correction is based on the method of Soulsby (1997), in which the threshold
bed shear stress θco is modified as a function of the bed slope χ, the angle of repose of the
sediment φs , modified in the Sisyphe steering file with the keyword FRICTION ANGLE OF THE
SEDIMENT (= 40◦ by default), and the angle of the current to the upslope direction ψ :
θc
cos ψ sin χ + (cos2 χ tan2 φs − sin2 ψ sin2 χ)0.5
=
,
θco
tan φs
(27)
with θc the modified threshold bed shear stress. Be aware that this formula effected the results
with threshold bedload formulae (like Meyer-Peter and Müller) only.
5.3.2 Correction of the direction of bedload transport rate (formula for deviation)
The change in the direction of solid transport is taken into account by the formula :
tan α = tan δ − T
∂Z f
,
∂n
(28)
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 n is the coordinate along the axis perpendicular to the flow.
Accessibility : Interne : EDF SA
Page 23 sur 73
c EDF SA 2013
EDF R&D
H-P74-2012-02004-EN
Reference manual
Sisyphe 6.3
Version 1.0
– DEVIA = 1 : according to Koch and Flochstra, the coefficient T depends exclusively on the
Shields parameter
4
(29)
T=
6θ
– DEVIA = 2 : the deviation is calculated based on Talmon et al. [44] and depends on the Shields
parameter and an empirical coefficient β 2
T=
1
√
(30)
β2 θ
The empirical coefficient β 2 can be modified by the keyword PARAMETER FOR DEVIATION (=
0.85 by default). For rivers could be achieved satisfying results by using a higher value of
β 2 = 1.5.
Keywords
–
–
–
–
–
–
SLOPE EFFECT (= NO, default option)
FORMULA FOR SLOPE EFFECT (SLOPEFF = 1, default option)
FORMULA FOR DEVIATION (DEVIA = 1, default option)
FRICTION ANGLE OF THE SEDIMENT (φs = 40◦ , default value)
BETA (β = 1.3, default value)
PARAMETER FOR DEVIATION (β 2 = 0.85, default value)
Figure 1 – Sloping bed : general slope (adapted from [42]).
5.4
Sediment Slide
An iterative algorithm prevents the bed slope to become greater than the maximum friction angle
(θs ≈ 32◦ − 40◦ ). This option is activated by use of key word SEDIMENT SLIDE. Further details can be
found in the Die-Moran’s thesis [15].
Keywords
– SEDIMENT SLIDE, = NO, default value
– FRICTION ANGLE OF SEDIMENT = 40◦ , default value
Accessibility : Interne : EDF SA
Page 24 sur 73
c EDF SA 2013
EDF R&D
5.5
Reference manual
Sisyphe 6.3
H-P74-2012-02004-EN
Version 1.0
Bedload transport in curved channels
The bedload movement direction deviates from the main flow direction due to helical flow effect [65]. Different authors proposed empirical formulas for evaluating this deviation. The Engelund
formula [19], based on the assumption that the bottom shear stress, the bed roughness and the mean
water depth are constant in the cross-section, is
tan δ = 7
h
r
(31)
where δ is the angle between the bedload movement and main flow direction, h the mean water
depth and r the local radius of curvature. Yalin and Ferrera da Silva [64] have showed that δ is
proportional to h/r. The local radius r can be computed based on the the cross sectional variation of
the free surface [19] :
U2
r = −ρα0 ∂Z ,
g ∂ys
with α0 a coefficient (α0 ≈ 0.75 for very rough beds and α0 ≈ 1.0 for smooth beds).
Keywords
– SECONDARY CURRENTS (= NO, default option)
– SECONDARY CURRENTS ALPHA COEFFICIENT (= 1.0, default option)
Attention
This option should only be activated when the flow is calculated by a depth-averaged model.
Accessibility : Interne : EDF SA
Page 25 sur 73
c EDF SA 2013
EDF R&D
6
H-P74-2012-02004-EN
Reference manual
Sisyphe 6.3
Version 1.0
Suspended load
6.1
6.1.1
Main assumptions
Vertical concentration profile
We assume a Rouse profile for the vertical concentration distribution, which is theoretically valid
in uniform steady flow conditions :
C (z) = Cre f
z−h a
z a−h
R
,
(32)
where R is the Rouse number defined by :
R=
ws
,
κu∗
(33)
with ws > 0 the vertical-settling sediment velocity, κ the von Karman constant (κ = 0.4), u∗ the
friction velocity corresponding to the total bed shear stress, and a the reference elevation above the
bed elevation.
6.2
Two-dimensional sediment transport equation
The two-dimensional (2D) sediment transport equation for the depth-averaged suspended-load
concentration C = C ( x, y, t) is obtained by integrating the 3D sediment transport equation over the
suspended-load zone. By applying the Leibniz integral rule, adopting suitable boundary conditions
and assuming that the bedload zone is very thin, the depth-integrated suspended-load transport
equation is obtained :
∂
∂
∂C
∂C
∂hC ∂(hUC ) ∂(hVC )
+
+
=
hes
+
hes
+ E − D,
(34)
∂t
∂x
∂y
∂x
∂x
∂y
∂y
with h = Zs − Z f ≈ Zs − Zzre f the water depth, assuming that the bedload layer thickness is very
thin, (U, V ) are the depth-averaged velocity components in the x − and y−directions, respectively.
The net sediment flux E − D is determined based on the concept of equilibrium concentration,
see [11] :
(35)
( E − D ) Zre f = ws Ceq − Cre f
where Ceq is the equilibrium near-bed concentration and Cre f is the near-bed concentration, calculated at the interface between the bed-load and the suspended load, z = Zre f . The reference elevation
Zre f corresponds to the interface between bedload and suspended load.
Equation (34), can be expressed in its non-conservative form as :
∂C
∂C
∂C
∂
∂C
∂C
1 ∂
E−D
+U
+V
=
hes
+
hes
+
,
∂t
∂x
∂y
h ∂x
∂x
∂y
∂y
h
|
{z
} |
{z
}
advection
(36)
diffusion
where U and V are the depth-averaged convective flow velocities in the x − and y− directions,
respectively.
Keywords
– SUSPENSION = NON, default value
Accessibility : Interne : EDF SA
Page 26 sur 73
c EDF SA 2013
EDF R&D
H-P74-2012-02004-EN
Reference manual
Sisyphe 6.3
Version 1.0
Attention
In Equation (34), the net sediment flux E − D is expressed in ms−1 . For consistency with
bedload units, expressed in m2 s−1 ), concentration is dimensionless. However, the user can
choose concentration by mass per unit volume of the mixture (kg m−3 ) by using the keyword
MASS CONCENTRATION. The relation between concentration by volume and concentration by
mass is Cm (kg m−3 ) = ρs C.
6.3
Convection velocity
In Equation (34) it is supposed that Qsx = hUC and Qsy = hVC. A straightforward treatment
of the advection terms would imply the definition of an advection velocity and replacement of the
depth-averaged velocities in Eq. (34). For the x − component :
Qsx = hUconv C
A correction factor is introduced in Sisyphe , defined by :
Fconv =
Uconv
U
The convection velocity should be smaller than the mean flow velocity (Fconv ≤ 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 Fconv :
B
I1
I2 − ln 30
,
Fconv = −
I1 ln eB
30
with B = k s /h = Zre f /h and
I1 =
Z 1
(1 − u ) R
B
u
du,
I2 =
Z 1
B
ln u
(1 − u )
u
R
du.
A similar treatment is done for the y−component. For further details, see [27].
Keywords
– CORRECTION ON CONVECTION VELOCITY (= NO, default option)
6.4
6.4.1
Numerical treatments
Treatment of the diffusion terms
According to the choice of the parameter OPDTRA of the keyword OPTION FOR THE DIFFUSION OF
TRACER, the diffusion terms in (36) can be simplified as follows :
– OPDTRA = 1 : the diffusion term is solved in the form ∇ · (ε s ∇ T )
– OPDTRA = 2 : the diffusion term is solved in the form 1h ∇ · (hε s ∇ T )
The value of the dispersion coefficient es depends on the choice of the parameter OPTDIF of the
keyword OPTION FOR THE DISPERSION :
– OPTDIF = 1 : the values of the constant longitudinal and transversal dispersivity coefficients
(T1 and T2 respectively) are provided with the keywords DISPERSION ALONG THE FLOW and
DISPERSION ACROSS THE FLOW
Accessibility : Interne : EDF SA
Page 27 sur 73
c EDF SA 2013
EDF R&D
Reference manual
Sisyphe 6.3
H-P74-2012-02004-EN
Version 1.0
– OPTDIF = 2 : the values of the longitudinal and transversal dispersion coefficients are computed with the Elder model T1 = αl u∗ h and T2 = αt u∗ h, where the coefficients αl and αt can be
provided with the keywords DISPERSION ALONG THE FLOW and DISPERSION ACROSS THE FLOW.
– OPTDIF = 3 : the values of the dispersion coefficients are provide by Telemac-2d
The diffusion term can also be set to zero with keyword DIFFUSION = NO (YES is the option by
default).
Attention
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 conservative schemes.
6.4.2
Treatment of the advection terms
The choice of the numerical scheme is based on keyword TYPE OF ADVECTION. The following
numerical methods are available in Sisyphe :
– Method of characteristics (1)
– Unconditionally stable and monotonous
– Diffusive for small time steps
– Not mass conservative
– Method Streamline Upwind Petrov Galerkin SUPG (2)
– Courant number criteria
– Not mass conservative
– Less diffusive for small time steps
– Conservative N-scheme (similar to finite volumes) (3, 4)
– Solves the continuity equation under its conservative form
– Recommended when the correction on convection velocity is
accounted for
– Courant number limitation (sub-iterations are included to reduce the time step
– Numerical schemes 13 and 14 are the same as 3 and 4 but here adapted to the presence of tidal
flats based on positive water depth algorithm
– Distributive schemes (PSI)
– like scheme 4, but the fluxes are corrected according to the tracer value : this relaxes the
courant number criteria and it is also less diffusive than scheme 4 and 14
– The CPU time is however increased
– This method should not applied for tidal flats
Attention
It is recommended to use scheme 4 or 14 (if tidal flats are present) as a good compromise
(accuracy/computational time)
Keywords
–
–
–
–
–
–
–
–
–
TETA SUSPENSION (= 1, default option)
TYPE OF ADVECTION (= 1, default option)
SOLVER FOR SUSPENSION (= 3, conjugate gradient by default option)
PRECONDITIONING FOR SUSPENSION (= NO, default option)
SOLVER ACCURACY FOR SUSPENSION (1.0 × 10−8 , default option)
MAXIMUM NUMBER OF ITERATIONS FOR SOLVER FOR SUSPENSION (= 50, default option)
OPTION FOR THE DISPERSION (= 1, constant dispersion coefficient by default option)
DISPERSION ALONG THE FLOW (1.0 × 10−2 , default option)
DISPERSION ACROSS THE FLOW (1.0 × 10−2 , default option)
Accessibility : Interne : EDF SA
Page 28 sur 73
c EDF SA 2013
EDF R&D
6.5
6.5.1
Reference manual
Sisyphe 6.3
H-P74-2012-02004-EN
Version 1.0
Equilibrium concentrations
Erosion and deposition rates
For non-cohesive sediments, the net sediment flux E − D is determined based on the concept of
equilibrium concentration, see [11] :
(37)
( E − D ) Zre f = ws Ceq − Cre f
where Ceq is the equilibrium near-bed concentration and Cre f is the near-bed concentration, calculated at the interface between the bed-load and the suspended load, z = Zre f . The reference elevation
Zre f corresponds to the interface between bedload and suspended load. For flat beds, the bed-load
layer thickness is proportional to the grain size, whereas when the bed is rippled, the bedload layer
thickness scales with the equilibrium bed roughness (kr ). The definition of the reference elevation
needs also to be consistent with the choice of the equilibrium near-bed concentration formula.
The keyword REFERENCE CONCENTRATION FORMULA access to four differents semi-empirical formulas :
– ICQ= 1 : Zyserman and Fredsoe formula
Ceq =
0.331(θ 0 − θc )1.75
,
1 + 0.72(θ 0 − θc )1.75
(38)
where θc is the critical Shields parameter and θ 0 = µθ ,the non-dimensional skin friction which
is related to the Shields parameter, by use of the ripple correction factor.
According to Zyserman and Fredsoe, the reference elevation should be set to Zre f = 2d50 . In
Sisyphe we take Zre f = k0s , where k0s is the skin bed roughness for flat bed (k0s ≈ d50 ) , the
default value of proportionality factor is KSPRATIO =3).
– ICQ= 2 : Bijker (1992) formula The equilibrium concentration corresponds to the volume concentration at the top of the bed-load layer. It can related to the bed load transport rate :
Ceq =
Qs
,
bZre f u∗
(39)
with an empirical factor b = 6.34.
According to Bijker, Zre f = kr .
Attention
The near bed concentration computed with the Bijker formula is related to the bedload.
Therefore, this option cannot be used without bedload transport (BED-LOAD = YES). Furthermore, both bedload and suspended load transport must be calculated at each time
step, with the keyword PERIOD OF COUPLING set equal to 1.
– ICQ= 3 : van Rijn formula
Ceq = 0.015d50
(τp /τc − 1)3/2
,
Zre f D∗0.3
(40)
The concentration is defined at Zre f = max(k s /2; 0.01m)
– ICQ= 4 : Soulsby-van Rijn formula
For details, we refer to the reference [42].
6.5.2
Ratio between the reference and depth-averaged concentration
By depth-integration of the Rouse profile (41), the following relation can be established between
the mean (depth-averaged) concentration and the reference concentration :
Cre f = FC,
Accessibility : Interne : EDF SA
Page 29 sur 73
c EDF SA 2013
EDF R&D
H-P74-2012-02004-EN
Reference manual
Sisyphe 6.3
where :
F
−1
=
Zre f
h
R Z
1
Version 1.0
Zre f /h
1−u
u
R
du.
(41)
In Sisyphe , we use the following approximate expression for F :
F −1 =
1
B R 1 − B (1− R ) ,
(1 − Z )
F −1 = − B log B,
for R 6= 1
for R = 1,
with B = Zre f /h.
Keywords
– REFERENCE CONCENTRATION FORMULA =1, default value
6.6
6.6.1
Initial and boundary conditions for sediment concentrations
Initial conditions
The initial values of volume concentration for each class can be either imposed in the subroutine
condim susp.f or specified in the steering file through the keyword INITIAL SUSPENSION CONCENTRATIONS.
Attention
It will not be used if EQUILIBRIUM INFLOW CONCENTRATION = YES.
6.6.2
Boundary conditions
For boundary conditions, the concentration of each class can be specified in the steering file
through keyword CONCENTRATION PER CLASS AT BOUNDARIES. To avoid unwanted erosion or sedimentation at the entrance of the domain, it may be also convenient to use keyword EQUILIBRIUM
INFLOW CONCENTRATION = YES to specify the value of the concentration at inflow, according to the
choice of the REFERENCE CONCENTRATION FORMULA. The depth-averaged equilibrium concentration is
then calculated assuming equilibrium concentration at the bed and a Rouse profile correction for the
F factor. Input concentrations can be also directly specified in the subroutine conlit.f. In conlit.f,
it is assumed that LIEBOR = LICBOR.
Accessibility : Interne : EDF SA
Page 30 sur 73
c EDF SA 2013
EDF R&D
7
H-P74-2012-02004-EN
Reference manual
Sisyphe 6.3
Version 1.0
Bed evolution equation
7.1
Bedload
To calculate the bed evolution , Sisyphe solves the Exner equation :
(1 − n )
∂Z f
+ ∇ · Qb = 0,
∂t
(42)
where n is the non-cohesive bed porosity (specified by the keyword NON COHESIVE BED POROSITY, =
0.40 by default), Z f the bottom elevation, and Qb (m2 s−1 ) the solid volume transport (bedload) per
unit width.
Equation (42) states that the variation of sediment bed thickness can be derived from a simple
mass balance and it is valid for equilibrium conditions.
7.2
Boundary conditions
From v6.2, as a general case the user has to specify two different boundary conditions files :
a conlim file for the hydrodynamics module (e.g., for Telemac-2d ) and a different conlim file for
Sisyphe . This allows the user to apply different conditions for a tracer in Telemac-2d and for bed
evolution or concentration in Sisyphe . However, in most simple applications, the same boundary
condition can be used.
In the sediment transport boundary condition file, the following boundary condition types can
be imposed on the open boundaries :
Imposed bed evolution
For this case, LIEBOR=5 (KENT) and LIQBOR=4 (KSORT) for boundaries with prescribed variation
of elevation. For example, if the user wants to define a boundary with a fixed bed elevation (zero
evolution), the conlim file should look like :
(not used)
(4)
LIQBOR
4
(not used)
(5)
...
...
LIEBOR
5
EBOR
0.
Imposed sediment transport rates
For this case, LIEBOR = 4 (KSORT) and LIQBOR = 5 (KENT). For example, if the user wants to define
a boundary with an imposed transport rate = 1.0 × 10−5 m2 s−1 , the conlim file should look like :
(not used)
(-)
LIQBOR
5
(not used)
(-)
Q2BOR
1.E-5
...
...
LIEBOR
4
EBOR
0.
For this case, a boundary condition file for Sisyphe is needed. LIHBOR is actually not used but is now
available for users. LIQBOR is the boundary condition on solid discharge. Q2BOR is the prescribed
solid discharge given at boundary points, and expressed in m2 s−1 .
Alternatively, the user can give the width-integrated values of the imposed sediment transport
along the boundaries by using the keyword PRESCRIBED SOLID DISCHARGES, which is in m3 s−1 .
When the keyword PRESCRIBED SOLID DISCHARGES is given it supersedes Q2BOR, which is then taken
as a constant profile.
For non-uniform sediment distribution, the imposed values of transport rate are distributed to
every class with the help of the fractions (array AVAIL), this is done for EBOR and QBOR in conlis.f
and may be changed.
Accessibility : Interne : EDF SA
Page 31 sur 73
c EDF SA 2013
EDF R&D
Reference manual
Sisyphe 6.3
H-P74-2012-02004-EN
Version 1.0
Keywords
– BEDLOAD = YES (default option)
– NON COHESIVE BED POROSITY (default option = 0.4).
7.3
Treatment of rigid beds and tidal flats
Non-erodable beds are treated numerically by limiting bedload erosion and letting incoming
sediment pass over. The problem of rigid beds is conceptually trivial but numerically complex.
For finite elements the minimum water depth algorithm allows a natural treatment of rigid
beds, see [25]. The sediment is managed as a layer with a depth that must remain positive, and
the Exner equation is solved similarly to the shallow water continuity equation in the subroutine
positive depths.f.
For finite volume schemes, the treatment must be explicitly specified and different options are
available with the keyword OPTION FOR THE TREATMENT OF NON ERODABLE BEDS :
Keywords
– OPTION FOR THE TREATMENT OF NON ERODABLE BEDS (= 0, default option)
– = 0 : erodable bottoms everywhere (default option)
– = 1 : minimisation of the solid discharge
– = 2 : nul solid discharge
– = 3 : minimisation of the solid discharge in FINITE ELEMENTS / MASS-LUMPING
– = 4 : minimisation of the solid discharge in FINITE VOLUMES
Attention
For the finite element method, the treatment of rigid bed is implicitly done by the numerical
scheme and previous options are now obsolete.
The space location and position of the rigid bed can be changed in the subroutine noerod.f. By
default, the position of the rigid bed is located at z = −100m.
7.4
Tidal flats
Tidal flats are areas of the computational domain where the water depth can become zero during
the simulation.
For finite elements the minimum water depth algorithm allows a natural treatment of tidal flats,
see [25] and the Exner equation is solved similarly to the shallow water continuity equation in the
subroutine positive depths.f.
For finite volume schemes, tidal flats can be treated with the keyword TIDAL FLATS. The companion keyword OPTION FOR THE TREATMENT OF TIDAL FLATS allows to choose between two different
approaches :
– OPTBAN = 1 : is the default option, and the equations are solved everywhere. The water depth
is set to a minimum value, controlled by the keyword MINIMAL VALUE OF THE WATER HEIGHT,
with a default value = 1.D-3m.
– OPTBAN = 2 : this option removes from the computational domain the elements with points
that present a water depth less than a minimum value. This option should be avoided because
it is not exactly mass-conservative.
Accessibility : Interne : EDF SA
Page 32 sur 73
c EDF SA 2013
EDF R&D
Reference manual
Sisyphe 6.3
H-P74-2012-02004-EN
Version 1.0
Keywords
– TIDAL FLATS (= YES, default option)
– MINIMAL VALUE OF THE WATER HEIGHT (= 1.D-3m, default value)
– OPTION FOR THE TREATMENT OF TIDAL FLATS (OPTBAN = 1, default option)
Attention
For the finite element method, the treatment of tidal flats is implicitly done by the numerical
scheme and previous options are now obsolete.
7.5
Bed evolution for suspension
By considering only suspended load sediment transport, the bed evolution is function of the net
vertical sediment flux at the interface between the bedload and the suspended load given by :
(1 − n )
∂Z f
+ ( E − D )z=Zre f = 0
∂t
(43)
with Z f the bed elevation, Zre f the interface between the bedload and suspended load, and n the
bed porosity. Once the flow variables are determined by solving the hydrodynamics and suspended
sediment transport equations, changes of bed level are computed from Equation (43) for the cell
coincident to the bed, calculating at each time step the net sediment flux E − D at the bedloadsuspended load interface (z = Zre f ), as explained later. Further details on the derivation of the mass
balance equation (43) and coupling strategies can be found in [65].
Accessibility : Interne : EDF SA
Page 33 sur 73
c EDF SA 2013
EDF R&D
8
Reference manual
Sisyphe 6.3
H-P74-2012-02004-EN
Version 1.0
Wave effects
8.1
Introduction
In coastal zones, the effect of waves superimposed to a mean current (wave-induced or tidal),
modifies the bed structure. Due to the reduced thickness of the bed boundary layer, the bottom
shear stress increases largely and the resulting sand transport rate could be in many cases of one
order of magnitude greater than in the case of currents alone. Furthermore, the wave-generated
ripples may also have a strong effect on the bed roughness and sand transport mechanisms. In
Sisyphe , these effects can be incorporated into the numerical simulations when the keyword EFFECT
OF WAVES is activated.
To compute sediment transport rates due to the action of waves, the wave height, period and
direction need to be specified. This information can be provided from a fortran file, by the subroutine
condim sisyphe.f, from a file containing the variables computed previously by the wave module,
e.g. Tomawac , or by internal coupling.
The wave orbital velocity U0 is calculated by Sisyphe by assuming the validity of the linear
theory :
Hs ω
,
U0 =
2 sinh(kh)
where Hs is the significant wave height, ω = 2π/Tp is the intrinsic angular frequency, k = 2π/L
is the wave number, with L the wave length. The wave number is calculated from the dispersion
relation (neglecting surface tension and amplitude effects) as :
ω 2 = gk tanh(kh).
8.1.1
Wave information
The wave information, required as input parameters for the morphodynamic model, are the
significant wave height Hs , the peak wave period Tp and the wave direction θw . If θw is not provided,
Sisyphe assumes the default value θw = 90◦ , which means that the direction of propagation of the
wave is parallel to the x-axis. When using the Bailard’s transport formula, this variable must be
always provided.
If the wave field is given from a file containing the variables computed previously on the same
mesh by the wave module, e.g. Tomawac , the keyword WAVE FILE must be specified in the Sisyphe
steering file. This file must contain the significant wave height Hs (HM0), the wave peak period Tp
(TPR5) and mean wave direction θw relative to y-axis (DMOY).
Attention
For the keyword WAVE FILE, variables are given by the last record of the file.
If the currents and wave field is provided from a file containing the variables computed previously on the same mesh by the currents and wave modules, e.g. Telemac-2d and Tomawac ,
respectively, the keyword HYDRODYNAMIC FILE must be specified in the Sisyphe steering file. This
file must contain the hydrodynamic variables (water height and velocity field) as well as the wave
variables : significant wave height Hs (HM0), wave peak period Tp (TPR5) and mean wave direction
θw relative to y-axis (DMOY). The users has also the choice of provide the wave variables from the
subroutine condim sisyphe.f.
Attention
For the keyword HYDRODYNAMIC FILE, the variables are given by the last record of the file if
the case is steady or by the time steps describing the tide or flood for the unsteady case.
Accessibility : Interne : EDF SA
Page 34 sur 73
c EDF SA 2013
EDF R&D
8.2
8.2.1
H-P74-2012-02004-EN
Reference manual
Sisyphe 6.3
Version 1.0
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 Uw by using the quadratic friction coefficient f w due to waves :
τw =
1
2
ρ f w Uw
.
2
The wave friction factor f w is calculated as a function of the orbital amplitude on the bed A0 = Uw /ω
and the bed roughness k s :
f w = f w ( A0 /k s ) ,
In Sisyphe , the formula proposed by Swart et al. [43] is provided :
−0.19
exp −6.0 + 5.2 Aks0
,
if A0 /k s > 1.57
fw =
0.30,
otherwise
8.2.2
Wave-current interaction
For combined waves and currents, the wave-induced (mean and maximum) bottom stresses can
be of one 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 stresses τcw , as a function of the bottom
shear stress due to currents only τc and the maximum shear stress due to waves only τw . Following
Bijker [7] :
1
τcw = τc + τw .
(44)
2
The mean τmean and maximum τmax shear stresses can be computed following the Soulsby’s
method [42]. Non-dimensional variables are defined :
X=
τ0
;
τ0 + τw
Y=
τmean
;
τ0 + τw
Z=
τmax
τ0 + τw
They can be parameterized as :
Y = X (1 + bX p (1 − X )q ) ,
Z = 1 + aX m (1 − X )n
The various model coefficients (a, b, m, n, p, q) are empirical functions of friction coefficients f w and
CD , and the wave-current angle φ :
2 fw
I
I
,
a = ( a1 + a2 |cos φ| ) + ( a3 + a4 |cos φ| ) log10
CD
2 fw
b = (b1 + b2 |cos φ| J ) + (b3 + b4 |cos φ| J ) log10
,
CD
......
with similar expressions for m, n, p, and q. The fitted constants (ai , bi , mi , ni , pi , qi , I and J) depend
on the turbulence model selected. Table 2 shows the coefficients used in Sisyphe corresponding to
the model of Huynh-Thanh and Temperville [28] for I = 0.82, J = 2.7.
Accessibility : Interne : EDF SA
Page 35 sur 73
c EDF SA 2013
EDF R&D
H-P74-2012-02004-EN
Reference manual
Sisyphe 6.3
i
1
2
3
4
ai
−0.07
1.87
−0.34
−0.12
bi
0.27
0.51
−0.10
−0.24
mi
0.72
−0.33
0.08
0.34
Version 1.0
ni
0.78
−0.23
0.12
−0.12
pi
−0.75
0.13
0.12
0.02
qi
0.89
0.40
0.50
−0.28
Table 2 – Constants issued from the turbulence model of Huynh-Thanh and Temperville [28], see
also Soulsby [42], page 91.
8.2.3
Friction coefficient under combined waves and currents
The characteristic shear stress
−
τ→
cw ( t ) , is related to the time-averaged mean velocity :
−−→ 2
−
→
,
τcw (t) = ρ f cw U (t)
where
−−→
U (t)
2
1
= Uc2 + Uw2
2
According to Camenen [8], the characteristic shear stress is taken as a weighted average between
τmean and τmax :
−
τ→
= Xτmean + (1 − X )τmax ,
cw ( t )
which is equivalent to :
−
τ→
cw ( t )
= Yτc + Zτw .
Finally, the expression for the wave-current interaction factor is :
f cw =
Yτc + Zτw
.
2
Uc2 + 12 Uw
(45)
This expression is used by Bailard [2] and Dibajnia and Watanabe [14] to compute sand transport
rates.
8.3
Wave-induced ripples
Equilibrium ripples are generally observed outside the surf zone. Their dimensions can be predicted as a function of waves (orbital velocity U0 and wave period T = 2π/ω), for a given uniform
sediment diameter d50 . In Sisyphe , the procedure of Wiberg and Harris has been implemented [60].
This formulation is only applicable to 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 D0 /d50 , with D0 = 2U0 /ω. The ripples wave length λ and height
η, can be estimated as follows :
– Under moderate wave conditions (orbital regime), the ripple dimensions are proportional to
D0 :
λ = 0.62D0 , η = 0.17λ.
– For larger waves (anorbital regime), ripple dimensions scale with the sand grain diameter :
λ = 535 d50 ,
D0
η = λ exp −0.095 log
η
Accessibility : Interne : EDF SA
2
!
D0
+ 0.042 log
− 2.28 .
η
Page 36 sur 73
c EDF SA 2013
EDF R&D
H-P74-2012-02004-EN
Reference manual
Sisyphe 6.3
Version 1.0
The effect of ripples is to increase the bed roughness. Based on the Bijker’s formula, in Sisyphe it is
assumed that :
k s = max (η, 3d50 ) .
Following Tanaka and Dang [49], the effect of the mean current superimposed to the waves is
accounted by introducing a correction factor α to the orbital velocity :
2.5 U 1.9
α = 1 + 0.81 tanh(0.3S∗2/3 )
,
Uw
with S∗ = d50
8.4
p
(s − 1) gd50 /4ν.
Wave-induced sediment transport
The effect of waves modifies the sand transport rates when it is superimposed to currents. This
effect can be accounted by using one of the following wave sand transport formula implemented
in Sisyphe , namely the Bijker’s formula [7], the Soulsby-van Rijn’s formula [42], the Bailard’s formula [2] and the Dibajnia and Watanabe’s formula [14].
The keyword BED-LOAD TRANSPORT FORMULA allows the user to choose :
– ICF = 4 : the Bijker’s formula (1968) can be used to determine the total transport rate Qs
(bed-load + suspension), with its two components the bedload Qsc and the suspended load
Qss computed separately. For bedload transport, Bijker extended the steady bed-load formula
proposed by Frijlink (1952), by adding the wave effects. In non-dimensional form, the bedload
transport rate is :
1
,
Φb = bθc0.5 exp −0.27
µθcw
where θc is the non-dimensional shear stress due to currents alone, θcw is the non-dimensional
shear stress due to wave-current interaction, and µ is a correction factor which accounts for the
effect of ripples. The shear stress under combined wave and current is calculated by Equation
(44).
In the original formulation, the coefficient b = 5. Recent studies showed that b = 2 provides
better results when comparing with field and experimental data. By default, in Sisyphe b =
2 but this value can be modified with the keyword B VALUE FOR THE BIJKER FORMULA. The
ripple factor correction µ is calculated in the same way as in the Meyer-Peter-Mueller formula
for currents only.
The suspended load transport is solved in a simplified manner, by assuming the concentration
profile to be in equilibrium. The inertia effects are not modeled and it is assumed that no
exchange takes place with the bed load layer.
After depth-integration and by assuming a Rouse profile for the concentration and a logarithmic velocity profile for the mean velocity profile, the suspended load can be written as :
Qss = Qsc I,
where
B A −1
I = 1.83 × 0.216
(1 − B ) A
with
A=
Accessibility : Interne : EDF SA
Ws
,
κu∗
Z 1
1−y A
r
u∗ =
y
B
τcw
,
ρ
Page 37 sur 73
ln
33y
B
dy,
B = k s /h.
c EDF SA 2013
EDF R&D
H-P74-2012-02004-EN
Reference manual
Sisyphe 6.3
Version 1.0
– ICF = 5 : Soulsby-van Rijn formula, the transport rate due to the combined action of waves
and current is provided by the following equation :
"
Qb,s = Ab,s U
0.018 2
U
U +2
CD 0
2
#2.4
0.5
− Ucr
.
This formula can be applied to estimate the components of the total sand transport rate (bedload Qb and suspended load Qs ), and it is suitable for beds covered by ripples. The bedload
and suspended load coefficients Ab,s are computed as :
Ab =
0.005h (d50 /h)1.2
((s − 1) gd50 )1.2
0.012d50 D∗−0.6
As =
,
((s − 1) gd50 )1.2
where U is the depth-averaged current velocity, U0 is the RMS orbital velocity of waves, and
CD is the quadratic drag coefficient due to current alone. This formula has been validated
assuming a rippled bed roughness with k s = 0.18 m. The critical entrainment velocity Ucr is
given by the expression :
4h
log
,
if 0.1mm ≤ d50 ≤ 0.5mm
0.19d0.1
50
10 D
90
Ucr =
4h
,
if 0.5mm ≤ d50 ≤ 2.0mm.
8.5d0.6
50 log10
D90
The diameter D90 , characteristic of the coarser grains, is specified with the keyword D90. The
validity range for the Soulsby-van Rijn formula is h = [1 − 20] m, U = [0.5 − 5] ms−1 , and
d50 = [0.1 − 2.0] mm.
– ICF = 8 : the Bailard’s formula is based on the 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 U (t), as follows :
−
→ 2−
→
f cw
ec
Qb =
U U
g(s − 1) tan ϕ
→
−
→ 3−
f cw
es
Qs =
U U ,
g(s − 1) Ws
with f cw the friction coefficient which accounts for wave-current interactions, es = 0.02, ec = 0.1
empirical factors, ϕ sediment friction angle (tan ϕ = 0.63) and < · > time-averaged over a
wave-period.
−−→
Under combined wave and currents, the time-varying velocity vector U (t) is composed of a
mean component Uc and an oscillatory component of amplitude U0 , assuming linear waves ; φ
is the angle between the current and the wave direction, φc is the angle between the mean current direction and the x −axis, and φw is the angle between the wave direction of propagation
and the x-axis. The time-varying velocity field verifies :
−−→
U (t) = Uc expiφc +U0 cos ωt expiφw .
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 (φ =
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 [8]. However, this
method has showed to be inefficient. Expressions for the third and fourth-order terms are given
below :
Accessibility : Interne : EDF SA
Page 38 sur 73
c EDF SA 2013
EDF R&D
H-P74-2012-02004-EN
Reference manual
Sisyphe 6.3
Version 1.0
– Third-order term :
→
−
→
−
→ 2−
−
→ −
→ −
→
1
U U = Uc (Uc2 + U02 ) + U0 (Uc · U0 )
2
→
−
→ 2−
1
1
= Uc3 + Uc U02 (1 + cos 2φ) cos φc − Uc U02 sin 2φ sin φc
U U
2
2
x
2
→
−
→ −
1
1
= Uc3 + Uc U02 (1 + cos 2φ) sin φc + Uc U02 sin 2φ cos φc
U U
2
2
y
– Fourth-order term :
→
−
→ 3−
U U
1
(24U02 Uc2 + 8Uc4 + 3U04 ) cos φc
8
x
→
−
→ 3−
1
= (24U02 Uc2 + 8Uc4 + 3U04 ) sin φc
U U
8
y
=
– ICF = 9 : the Dibajnia and Watanabe (1992) formula is an unsteady formula, which accounts
for inertia effects. The sand transport rate is calculated by the expression :
−
→
−
→
Qs
Γ
= α 1− β ,
Ws d50
|Γ|
(46)
−
→
with α = 0.0001 and β = 0.55 empirical model coefficients and Γ 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. In the wave direction, the time-varying velocity is given by :
U (t) = Uc cos φ + U0 sin ωt,
r is defined by the asymmetry parameter :
r=
U0
.
Uc cos φ
The wave cycle is decomposed into two parts :
1. During the first part of the wave-cycle (0 < t < T1 ), the velocities are directed onshore
(U (t) > 0).
2. During the second part (T2 = T − T1 ), the velocities are negatives (U (t) < 0).
The sand transport rate in the wave direction is :
0
ΓX0 =
0
u1 T1 (Ω31 + Ω23 ) − u2 T2 (Ω32 + Ω13 )
( u1 + u2 ) T
where Ω1 and Ω2 represent the amount of sand transported in the onshore and offshore directions which will be deposited before flow reversal, respectively, Ω10 and Ω20 represent the
remaining part which stay in suspension after flow reversal. The inertia of sediment depends
on the ratio d50 /ws , and represented by parameter ωi :
ωi =
u2i
,
2(s − 1) gWs Ti
u2i =
2
Ti
with
Z
u>or<0
u2 dt,
and two different cases :
Accessibility : Interne : EDF SA
Page 39 sur 73
c EDF SA 2013
EDF R&D
H-P74-2012-02004-EN
Reference manual
Sisyphe 6.3
Version 1.0
1. All sediment is deposited before flow reversal ωi ≤ ωcrit
Ω i = ωi
2Ws Ti
,
d50
Ωi0 = 0
2. Part of the sediment stays in suspension after flow reversal ωi ≥ ωcrit
Ωi = ωcrit
2Ws Ti
,
d50
Ωi0 = (ωi − ωcrit )
2Ws Ti
.
d50
The critical value of ωcrit is calculated as function of the wave-current non-dimensional friction :
Θcw =
< τcw >
1
.
ρ
(s − 1) gd50
According to Temperville et Guiza :
ωcrit
0.03,
1 − p1 − ((Θ − 0.2)/0.58)2 ,
cw
p
=
0.8104
Θ
−
0.4225,
cw
p
0.7236 Θcw − 0.3162,
if Θcw ≤ 0.2
if 0.2 ≤ Θcw ≤ 0.4
if 0.4 ≤ Θcw ≤ 1.5
otherwise
Attention
In Sisyphe , the solid discharge 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.
8.5
Internal coupling waves-currents and sediment transport
The internal coupling waves-currents and sediment transport is implemented in the TelemacMascaret system. The keyword COUPLING WITH = ’TOMAWAC, SISYPHE’ (in the Telemac-2d steering
file) is used to internally couple with Tomawac and Sisyphe .
The coupling period with Tomawac and Sisyphe can be controlled with the keywords COUPLING
PERIOD FOR TOMAWAC and COUPLING PERIOD FOR SISYPHE, respectively.
In the Telemac-2d cas file, the keyword WAVE DRIVEN CURRENTS = YES allows to incorporate the
effect of waves in the hydrodynamic simulation.
In the Sisyphe cas file, the keyword EFFECT OF WAVES = YES is used to consider the effect of the
waves on the solid transport formula. In Sisyphe , the implemented bedload transport formulae that
consider the combined effect of currents and waves are Bijker (=4), Soulsby-van Rijn (=5), Bailard
(=8) and Dibajnia and Watanabe (=9). Formulae 4, 5 and 8 consider bedload and suspended load
transport. Formula 9 considers total transport.
Accessibility : Interne : EDF SA
Page 40 sur 73
c EDF SA 2013
EDF R&D
9
H-P74-2012-02004-EN
Reference manual
Sisyphe 6.3
Version 1.0
Sand grading effects
9.1
9.1.1
Sediment bed composition
Granulometry distribution
The number of classes of bed material is specified in the steering file with the keyword NUMBER
OF SIZE-CLASSES OF BED MATERIAL and is limited to NSICLA ≤ 20. For uniform sediments, 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 1 < j < NSICLA is defined by its mean diameter d50 (j) and volume fraction in
the mixture at every node AVAI(j). The characteristics of each class of sediment, such as the Shields
number θc (i) or the settling velocity ws (i) for each class can be specified in the steering file or
calculated by the model as for a single sediment class.
The percent of each class of material need to verify :
∑
AVAI(j) = 1
and
0 ≤ AVAI(j) ≤ 1.
(47)
j=1,NSICLA
The initial bed composition can be defined in the steering file for uniform bed. For complex
bed configurations, e.g. spatial variation, vertical structure, etc., the user subroutine init compo.f is
used to define the initial state. The mean diameter dm , ACLADM, is calculated as :
dm =
∑
AVAI(j) d50 (j).
(48)
j=1,NSICLA
9.1.2
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, AVAI(i,j,k) as well as the
thickness of each layer Es 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 active layer. The initial composition of the bed (number of layers, thickness
of each layer Es , composition of each layer AVAI are specified in user’s subroutine init compo.f.
The subroutine init avail.f verifies if the structure of the initial bed composition is compatible
with the position of the rigid bed, as defined in user’s subroutine noerod.f :
Z f (i) − Zr (i) =
∑
Es (k).
(49)
k=1,NOMBLAY
In general, the thickness of the bed is taken to be large (by default, 100 m), so the bottom layer
thickness of the last layer (nearest to the rigid bed) needs to be increased.
Assuming the porosity n of each class to be identical and constant, the total mass of sediments
per unit area is :
Ms ( i ) =
∑
ρs Es (k)(n − 1).
(50)
k=1,NOMBLAY
In each layer, the percent of each class AVAI, which is defined as the volume of each class of
material per the total volume of material, is considered to be a constant. For each layer and at each
point of the domain, the following constraints need to be satisfied :
Accessibility : Interne : EDF SA
Page 41 sur 73
c EDF SA 2013
EDF R&D
H-P74-2012-02004-EN
Reference manual
Sisyphe 6.3
∑
AVAI(i, k, j) = 1
Version 1.0
0 ≤ AVAI(i, k, j) ≤ 1.
and
j=1,NSICLA
Keywords
The initial sediment bed composition is defined by :
– NUMBER OF SIZE-CLASSES OF BED MATERIAL : NSICLA = 1, by default
– SEDIMENT DIAMETERS : FDM = .01, by default
– SETTLING VELOCITIES : if XWC is not given, the subroutine vitchu sisyphe.f is used to
compute the settling velocity based on the Stokes, Zanke or Van Rijn formulae as a function
of the grain size
– SHIELDS PARAMETERS : for multi grain size, the Shields parameter needs to be specified for
each class. If only one value is specified, the shields parameter will be considered constant
for all classes. The default option, no Shields parameter given in steering file, is to calculate
the Shields parameter as a function of the sand grain diameter, see logical CALAC.
– INITIAL FRACTION FOR PARTICULAR SIZE CLASS : AVA0 = 1.;0.;0.;..., by default
For more than one-size classes, the user should specify NSICLA values for the mean diameter,
initial fraction, etc. For example, for a mixture of two classes :
Example
–
–
–
–
–
NUMBER OF SIZE-CLASSES OF BED MATERIAL = 2
SEDIMENT DIAMETERS = 0.5; 0.5
SETTLING VELOCITIES = 0.10; 0.05
SHIELDS PARAMETERS = 0.045; 0.01
INITIAL FRACTION FOR PARTICULAR SIZE CLASS = 0.5; 0.5
Attention
When using a (manually modified) sedimentological restart pay attention, that the following
restrictions are valid :
∑
AVAI(j) = 1
and
0 ≤ AVAI(j) ≤ 1
(51)
j=1,NSICLA
Z f (i) − Zr (i) =
∑
Es (k)
(52)
k=1,NOMBLAY
9.2
Sediment 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 classical active layer model based on Hirano’s concept [23] and a newly developed continuous
vertical sorting model (C-VSM), implemented by Merkel and Kopmann (BAW) [36].
9.2.1
Active layer and stratum
Since only a relatively thin layer is envolved in the transport, the upper part (active part) of the
bed is therefore divided into a thin active layer and an active stratum. The composition of the active
layer is used to calculate the composition of transported material and the intensity of transport rates
Accessibility : Interne : EDF SA
Page 42 sur 73
c EDF SA 2013
EDF R&D
Reference manual
Sisyphe 6.3
H-P74-2012-02004-EN
Version 1.0
for each class of sediment :
Qs =
∑
AVAI(i, 1, j) Qs (j),
(53)
j=1,NSICLA
where AVAI(i,1,j) is the percentage of the class j in the active layer.
The active stratum is the layer that will exchange sediment with the active layer in order to keep
the active layer to a fixed (either constant or calculated) size. In case of high erosion depth, the active
stratum becomes thinner and thinner, since it will be merged with the stratum lying underneath.
The active layer thickness concept is not well defined and depends on the flow and sediment
transport characteristics [52]. According to Yalin [63], 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 few centimeters in the case of rippled beds. In
the presence of dunes, the active layer thickness should be half the dune height [39].
In Sisyphe , the active layer thickness is an additional parameter of the model, called 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 = NON, it is possible to use
a space and time varying active layer thickness during the simulation. In Sisyphe it is assumed that
ELAY0 = 3dm , with dm the mean diameter of the active layer. Using the continuous vertical sorting
model, more formulae for calculating the active layer thickness are available.
The erosion rate is restricted by the total amount of sediments in the active layer, which therefore acts as a rigid bed. The same method applied for rigid bed is applicable for the active layer
formulation.
Attention
When dealing with graded sediment and thin thickness active layers, it is advised to use finite
elements combined with the positive depth algorithm, as used for the treatment of rigid bed,
in order to avoid numerical problems such as negative sediment fractions. The error message
EROSION GREATER THAN ONE LAYER THICKNESS can also appear in the listing file when the
erosion is greater than the active layer thickness. If the sedimentation is greater than the active
layer thickness the error message DEPOSITION MORE THAN ONE LAYER THICKNESS appear in the
listing as well. In both cases it is strongly recommended to reduce the time step. Otherwise
the calculation of the fractions is not valid anymore.
9.2.2
Continuous vertical sorting model
The C-VSM overcomes many limitations of the classical layer implementation of Hirano (HRVSM). Even though it is a different way to manage the grain sorting, it is just another interpretation
of Hirano’s original idea with fewer simplifications. So still an active layer is defined but the grain
distribution in this layer is computed for 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 does not 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 Sisyphe steering file. The full C-VSM output can
be found in the new Selafin files VSPRES & VSPHYD in the tmp-folders. As the higher resolution
Accessibility : Interne : EDF SA
Page 43 sur 73
c EDF SA 2013
EDF R&D
H-P74-2012-02004-EN
Reference manual
Sisyphe 6.3
Version 1.0
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
The initial sediment bed composition is defined by :
– VERTICAL GRAIN SORTING MODEL = 1
– 0 = Layer = HR-VSM (Hirano + Ribberink as until Sisyphe v6p1)
– 1 = C-VSM
– C-VSM MAXIMUM SECTIONS = 100
– Should be at least 4 + 4× Number of fractions,
– better > 100, tested up to 10000
– C-VSM FULL PRINTOUT PERIOD = 1000
– 0 => GRAPHIC PRINTOUT PERIOD
– Anything greater 0 => Sets an own printout period for the CVSP
– Useful to save disk space
– C-VSM PRINTOUT SELECTION = 0|251|3514|...|...
– Add any 2D Mesh Point numbers for .CSV-Ascii 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
– C-VSM DYNAMIC ALT MODEL = 5
– Model for dynamic active layer thickness approximation
– 0 = constant (Use Keyword : ACTIVE LAYER THICKNESS)
– 1 = Hunziker & Guenther
ALT = 5d MAX
– 2 = Fredsoe & Deigaard (1992)
ALT =
2τB
(1 − n) g(ρS − ρ) tan Φ
– 3 = van RIJN (1993)
ALT = 0.3D∗0.7
τB − τC 0.5
d50
τB
– 4 = Wong (2006)
ALT = 5
τB
− 0.0549)0.56 d50
(ρS − ρ) gd
– 5 = Malcharek (2003)
ALT =
d90
τ
max(1, B )
1−n
τC
– 6 = 3*d50 within last time steps ALT’
ALT = 3d50
9.2.3
General formulation
Bedload transport rates are computed separately for each class using classical sediment transport
formulae, corrected for sand grading effects. The Exner equation is then solved for each class of
sediment. The individual bed evolution due to each class of bed material is then summed to give
the global evolution due to bedload. Similarly, the suspended transport equation is solved for each
class of sediment and the resulting bed evolution for each class is then summed to give the global
Accessibility : Interne : EDF SA
Page 44 sur 73
c EDF SA 2013
EDF R&D
H-P74-2012-02004-EN
Reference manual
Sisyphe 6.3
Version 1.0
evolution due to the suspended load.
At each time step, the bed evolution due to bedload and to suspension transport is used to
compute the new bed structure. The algorithm which determines the new bed composition considers
two possibilities, namely deposit or erosion, and ensures mass conservation for each individual class
of material. The algorithm is explained in [22].
9.2.4
Hiding exposure
In order to calculate bedload sediment transport rates for each class of sediment, it is necessary to
consider the effect of hiding exposure. That means that in a sediment mixture, bigger particles will
be more exposed to the flow than the smaller ones and therefore, prediction of sediment transport
rates with classical sediment transport formulas, need to be corrected by use of a hiding-exposure
correction factor. In Sisyphe , the keyword HIDING FACTOR FOR PARTICULAR SIZE CLASS sets the
value of hiding factor for a particular size class and the keyword HIDING FACTOR FORMULA allows
the user to choose among different formulations. Two formulas, Egiazaroff [16] and Ashida & Michiue [1], have been written based on the Meyer-Peter and Müller formula. Both formulas modify
the critical Shields parameter that will be used in the Meyer-Peter formula. Therefore they only can
be combined useful with threshold formulae. The formula of Karim and Kennedy [30] can be used
in combination with any bedload transport predictor. This formula modifies directly the bedload
transport rate.
– HIDFAC = 1 : formulation of Egiazaroff
θcr = 0.047ζ i ,
with ζ i =
log(19)
log(19Di /Dm )
2
– HIDFAC = 2 : formulation of Ashida & Michiue
– HIDFAC = 4 : formulation of Karim et al.
Keywords
Sand grading effects are defined by :
– HIDING FACTOR FORMULA : if HIDFAC = 0 (by default), the user needs to give HIDING FACTOR
FOR PARTICULAR SIZE CLASS
– ACTIVE LAYER THICKNESS : = 10000, by default
– CONSTANT ACTIVE LAYER THICKNESS : = YES, by default
– NUMBER OF BED LOAD MODEL LAYERS : NOMBLAY = 2, by default
9.2.5
Sediment transport formula
The formula of Hunziker [26] is an adaptation of the Meyer-Peter Müller formula for fractional
transport. The volumetric sediment transport per sediment class is given by :
Φb = 5pi [ζ i (µθi − θcm )]3/2
if
µθi > θcr ,
(54)
with pi the fraction of class i in the active layer, θi the Shields parameter, θcm = θcr ( Dmo /Dm )0.33 the
corrected critical Shields parameter, θcr the critical Shields parameter (θcr = 0.047), ξ i = ( Di /Dm )−α
the hiding factor, Di the grain size of class i, Dm the mean grain size of the surface layer (m),
Dmo the mean grain size of the under layer (m), and α a constant. 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 α ≤ 2.3 and Di /Dm ≥ 0.25.
Hunziker developed an own hiding-exposure formula corresponding to his bed load formula.
Accessibility : Interne : EDF SA
Page 45 sur 73
c EDF SA 2013
EDF R&D
Reference manual
Sisyphe 6.3
H-P74-2012-02004-EN
Version 1.0
Part III
Cohesive sediments
by Catherine Villaret, Lan Anh Van, Damien Pham van Bang & Pablo Tassi
Accessibility : Interne : EDF SA
Page 46 sur 73
c EDF SA 2013
EDF R&D
10
Reference manual
Sisyphe 6.3
H-P74-2012-02004-EN
Version 1.0
Cohesive sediment properties
10.1
Introduction
Fine cohesive sediments are predominant in most estuaries, basins and reservoirs. The prediction
of cohesive sediment transport, often associated with pollutants, is crucial regarding environmental
issues.
Cohesive properties appear for fine particles (silts and clay), with diameter less than a limiting
value of about 60 µm, depending on the physico-chemical properties of the fluid and salinity. Fine
cohesive sediments are transported in suspension (no bed-load) and transport processes strongly
depend on the state of floculation of the suspension and consolidation of the bed.
The erosion rate mainly depends on the degree of consolidation of the sediment bed, while the
settling velocity depends on the state of floculation and aggregates properties.
Attention
The separation value at 60µm to discriminate non-cohesive from cohesive sediment is conventional (UK classification). The value is different depending on the country (63µm in Netherlands, 75µm in USA as pointed by Winterwerp and Van Kesteren [61]). Moreover, aggregation
of flocs can lead to the formation of macro-flocs larger than 100µm.
Keywords
In Sisyphe , the simplest case of uniform and cohesive sediments is characterized by one single
value for the primary grain size D50 and density ρs = 2650 kg/m3 for quartz particles, which
is transported in suspension (no bedload).
In the Sisyphe steering file, the physical properties of the sediment are defined :
– COHESIVE SEDIMENT = YES (= NO, default option)
– SEDIMENT DIAMETERS (D50 > 0.00006 m, for non-cohesive sediment)
– SEDIMENT DENSITY (ρs = 2650 Kg/m3 , default value)
Both key words SUSPENDED LOAD = YES, BED LOAD = NO are set automatically, in order to be
consistent with the selected type of sediment, when COHESIVE SEDIMENT = YES
– SUSPENDED LOAD = YES (= NO, default option)
– BED LOAD = NO (= YES, default option)
10.2
Effect of floculation
For cohesive sediments (D50 < 60µm), the primary grain diameter and density no longer control
the transport rates. Due to electro-static attractive forces (van der Waals), individual particles tend to
form flocs. Their size can be an order of magnitude greater than the initial particle diameter, while
the density decreases : the floc characterictics (size and density, as well as its fractal dimension) are
more relevant to compute the settling velocity.
The settling velocity of cohesive particles therefore depends on the suspended sediment concentration as well as on other physico-chemical properties (pH, salinity, cations for instance) of the
suspension which influence the flocculation process. The effect of turbulence level determines also
the floculation state. The critical bed shear strength depends on the consolidation state or age of the
sediment bed (Migniot, 1968).
Floculation can be modelled using process based models (cf. Mc Anally [34]). The settling velocity
is generally one order of magnitude greater than the settling velocity based on the individual particle
size, as a result of flocculation.
The flocculation process is not yet programmed in Sisyphe and the settling velocity should be
specified by the user as a calibration parameter. More physically based relations will be implemented
Accessibility : Interne : EDF SA
Page 47 sur 73
c EDF SA 2013
EDF R&D
Reference manual
Sisyphe 6.3
H-P74-2012-02004-EN
Version 1.0
in a future release of the code.
Keywords
– SETTLING VELOCITIES
Attention
The settling velocity needs to be specified, if not, it will be calculated by the Stokes law as a
function of the individual particle diameter in subroutine vitchu sisyphe.f. As a result of
flocculation, the settling velocity can be approximately an order of magnitude greater.
10.3
Consolidation process
Once the sedimentation process is achieved, a sediment bed is formed. For non-cohesive bed, no
evolution with time will be observed if no erosion or further sedimentation occurs. For cohesive bed,
the concentration will increase with time as the result of self-weight consolidation or compaction.
Two stages, namely the primary and secondary consolidation are classically considered to explain
the progressive compaction of the cohesive sediment bed with time. As the formed bed resulting
from the sedimentation of mud is very soft, it tends to compact under its self-weight. During this
process, the void ratio tends to diminish so that pore water is expulsed upward as a result of water
incompressibility. The dynamic of this primary stage of consolidation is therefore controlled by
the permeability of the bed. At the end of the primary consolidation, the excess pore pressure is
fully dissipated (Figure 2) but the cohesive bed records further compaction. This secondary stage
of consolidation is no longer related to pore pressure but rather to the solid skeleton property. The
secondary consolidation is considered as a result of time effect (viscous or ageing) in the stress-strain
relationship of the soil. Two effects of time are currently considered in modern soil mechanic : the
creep of the soil (viscous behaviour of cohesive geomaterial) and ageing (thixotropy of mud). Even
through the dynamic of the secondary consolidation is much slower than the primary consolidation,
this stage is important for long term prediction.
Keywords
– MUD CONSOLIDATION = NON, by default
11
Cohesive bed structure
11.1
Multi-layer model
Cohesive sediment beds are generally non-uniform : as a result of self-weight consolidation, the
bed becomes stratified, with density increasing with distance from the surface. The top layer is
generally made of freshly deposited soft mud, while the bottom consolidated layers present higher
resistance to erosion. The vertical stratification of the cohesive bed is therefore a key issue, which
controls the amount of material to be put in suspension.
In Sisyphe , cohesive sediment bed can be represented by a fixed number NCOUCH TASS of layers.
The number of layers is fixed with maximum number up to 20. Each layer is characterized by its
concentration and resitance to erosion. The concentration of each layer (Cs ( j))is generally constant
and specified by keyword MUD CONCENTRATION PER LAYER expressed in kg/m3 . It may also vary in
space and time (see consolidation model 3). The resistance of each layer is specified by keyword
CRITICAL EROSION SHEAR STRESS OF THE MUD expressed in N/m2 .
Each layer j (1 ≤ j ≤ NCOUCH TASS) is defined by its concentration Cs , thickness Es and mass
Accessibility : Interne : EDF SA
Page 48 sur 73
c EDF SA 2013
EDF R&D
Reference manual
Sisyphe 6.3
H-P74-2012-02004-EN
Version 1.0
Figure 2 – Vertical variation of density and effective stress within the cohesive sediment bed.
per surface area Ms , which depends on both layer concentration and thickness, such that
Ms (i, j) = Cs ( j) Es (i, j)
where Ms is the mass per unit surface area (Kg/m2 ), Cs is the concentration (kg/m3 ) and Es the
layer thickness.
For specific applications, the bed concentrations Cs (i, j) can also vary from one point to another.
The initial distribution can be specified in subroutine init compo coh.f.
Keywords
– NUMBER OF LAYERS OF THE CONSOLIDATION MODEL : NCOUCH TASS (default value = 1)
should be less than the maximum number of layers (NLAYMAX = 20)
– MUD CONCENTRATION PER LAYER, in Kg/m3 : Cs (default value = 50)
Attention
The keywords NUMBER OF BED LOAD MODEL LAYERS and NUMBER OF LAYERS OF THE
CONSOLIDATION MODEL are essentially the same except that the default values are different.
For cohesive sediments, it is possible to have only one uniform layer, whereas for sand grading
algorithm we need at least 2 layers (the active layer and the substratum). The subroutine
lecdon.f specifies at the end NOMBLAY = NCOUCH TASS
Attention
Due to the default value of NOMBLAY for sand grading effects (NOMBLAY = 2). For uniform beds both keywords need to be specified :
– NUMBER OF BED LOAD MODEL LAYERS = 1
– NUMBER OF LAYERS OF THE CONSOLIDATION MODEL = 1
Accessibility : Interne : EDF SA
Page 49 sur 73
c EDF SA 2013
EDF R&D
H-P74-2012-02004-EN
Reference manual
Sisyphe 6.3
Version 1.0
Figure 3 – Multi-layer bed structure (NCOUCH TASS is the number of layers)
11.2
Initialization
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.
The initial layer thicknesses are specified in user-subroutine init compo coh.f. It is called by
subroutine init mixte.f which just checks if the sum of all layers is equal to the total bed thickness
as specified by the initial bathymetry Z f and rigid bed Zr :
Z f (i ) − Zr (i ) =
∑
Es (i, j)
(55)
j=1,NCOUCH TASS
where i is the number of point, j the number of layer, Z f the bed level and Zr the non-erodible bed
level.
Accessibility : Interne : EDF SA
Page 50 sur 73
c EDF SA 2013
EDF R&D
12
H-P74-2012-02004-EN
Reference manual
Sisyphe 6.3
Version 1.0
Erosion and deposition rates
The erosion and deposition rates are needed to solve the net sediment flux present in the
advection-diffusion equation 34, see Chapter 6. In this equation, the erosion term is treated explicitly
while the deposition term is treated implicitly.
Attention
The sediment concentration is assumed to be uniform over the depth, which is a reasonable
assumption for fine cohesive sediments (ws << u∗ ).
12.1
Erosion flux
Erosion fluxes can be calculated using the Partheniades classical formula which is programmed
in Sisyphe : erosion flux are related to the excess of applied bed shear stress to the bed shear strength
at the bed surface. When the bed is stratified, it is necessary to take into account the vertical increase
of bed shear strength as the bed gets eroded.
Uniform bed
The classical Partheniades formula is applied for cohesive sediments. Assuming the bed to be
uniform (NCOUCH TASS = 1) :
!
u∗ 2
M
−1 ,
if u∗ > u∗e
E=
u∗e
0,
otherwise
where u∗e is the critical erosion velocity and u∗ is the friction velocity related to skin friction, defined
as :
s
τ0
u∗ =
ρ
with τ 0 the bed shear stress corrected for skin friction and ρ the fluid density.
The Partheniades coefficient M is a dimensional coefficient. It is specified in the steering file with
units given in kg/m2 /s. The erosion rate is expressed in m/s, by converting the M coefficient to
m/s (subroutine lecdon.f).
Non uniform bed
In the multi-layer model, the critical erosion shear stress increases as the mud concentration
increases (see Migniot, 1968).
In the multi-layer model (NCOUCH TASS > 1), each layer is characterized by its density and critical
bed shear stress. The erosion rate of each layer E( j) decreases with distance from the surface.
For each layer j and at each time step, the erosion rate E( j) is calculated as a function of the
difference between the applied bed shear stress and the critical bed shear strength, using (4). For
each E( j) > 0, the layer is erodible.
The erosion flux E is determined as the mean erosion rate, averaged over the eroded depth :
E=
1
∆zer
Z ∆zer
0
E(z)dz =
1
ρs dt
Z dt
0
dMs
where dMs is the eroded mass per unit surface area (Kg/m2 ), ∆zer is the maximum depth to be
eroded, E(z) is the erosion flux (m/s) at distance z from the bed interface, and dt the time step.
Accessibility : Interne : EDF SA
Page 51 sur 73
c EDF SA 2013
EDF R&D
H-P74-2012-02004-EN
Reference manual
Sisyphe 6.3
Version 1.0
Iterative procedure
The potential depth to be eroded is estimated using the iterative procedure described below and
programmed in subroutine suspension erosion coh.f.
At each time step, the top layer is first eroded if E(1) > 0. Once it is empty, the next layer is
eroded etc... For each erodible layer j with erosion rate (E( j) > 0), the time interval dt( j) to erode it
completely is estimated by a simple mass balance from the mass of the layer Ms ( j) :
dt( j) =
Ms ( j )
Es ( j)Cs ( j)
=
ρs E( j)
ρs E( j)
For the first layer (j = 1), there are two possibilities :
– dt(1) > dt : layer 1 is emptied partially, such that −Cs (1)dEs (1) = ρs dtE( j)
– dt(1) < dt : layer 1 is emptied completly, and the next (second) layer can be eroded, such that
dEs (1) = − Es (1).
The last (non-empty) layer (noted jmax ) to be eroded is obtained by emptying all successive top layers
until :
jmax −1
∑
jmax
dt( j) < dt <
j =1
∑ dt( j)
j =1
For jmax , the time left to erode the last layer is
j max −1
∑
dtmax = dt −
dt( j)
j =1
The mass potentially eroded during this process is
j max −1
Ms =
∑
Ms ( j) + ρs E( jmax )dtmax
j =1
The erosion flux (m/s) is therefore estimated as follows :
E=
Ms
ρs dt
Keywords
– COHESIVE SEDIMENTS (= YES, default option)
– PARTHENIADES CONSTANT (M = 10−3 Kg/m2 , default option)
– CRITICAL EROSION SHEAR STRESS OF THE MUD (τce = 0.01 N/m2 , default option)
12.2
Deposition flux
The deposition flux is calculated as a function of the near bed concentration. In the case of
fine cohesive sediments (ws << u∗ where ws is the settling velocity, and u∗ , the turbulent friction
velocity), the bed concentration is approximately equal to the depth-averaged concentration. The
deposition flux can be treated as an implicit term, whereas the erosion flux is explicit.
D=
Accessibility : Interne : EDF SA
0
ws C
1−
u∗
u∗d
2 !
Page 52 sur 73
if u∗ > u∗d
otherwise
c EDF SA 2013
EDF R&D
H-P74-2012-02004-EN
Reference manual
Sisyphe 6.3
Version 1.0
where C is the depth-averaged concentration, ws , the settling velocity, and u∗d , the critical deposition
velocity which represents the limiting shear velocity, above which the sediment flocs are broken and
resuspended.
Keywords
– CRITICAL SHEAR VELOCITY FOR MUD DEPOSITION (u∗d = 1000 m/s, default option)
– SETTLING VELOCITIES (by default, ws is calculated by the model based on the Stokes law
and individual particle diameter)
13
Bed Evolution
The bed evolution for cohesive sediments is solved with the Equation (43), see Chapter 6.
Uniform bed
The uniform bed (NCOUCH TASS = 1) is characterized by a single density Cs = Cs (1). The bed
evolution at each time step is therefore :
dZ f =
ρs ( D − E)dt
Cs
(dZ f < 0 for net erosion and dZ f > 0 for net deposition).
Non-uniform beds
For non-uniform beds (NCOUCH TASS > 1), two different cases occur.
– In the case of net deposition (D − E > 0), sediment is deposited in the first top layer of density
Cs (1) :
dZ f = dEs (1) =
ρs ( D − E)dt
>0
Cs (1)
– In the case of net erosion (E − D > 0), sediment is eroded layer by layer (dEs ( j) < 0), until the
eroded mass balances the net eroded flux. The last layer to be eroded jmax is determined in
order to satisfy :
j max −1
∑
j max
Msn ( j) < ρs ( E − D )dt <
j =1
∑
Msn ( j)
j =1
Top layers are emptied (dEs ( j) = − Es ( j) for j = 1 to jmax − 1) ; the last layer is eroded up to a
the depth dEs ( jmax ) < 0 in order to satisfy mass conservation :
j max −1
∑
Msn ( j) − dEs ( jmax )Cs ( jmax ) = ρs ( E − D )dt
j =1
|
{z
}
eroded mass
The variation of bed elevation is therefore :
j max −1
dZ f = −
∑
!
Esn ( j)
+ dEs ( jmax ) < 0
j =1
The thickness and mass of each layer are updated at the end of each time step (n + 1), from their
initial values (n) :
for j = (1, jmax − 1)
Accessibility : Interne : EDF SA
Esn+1 s( j) = 0
Msn+1 ( j) = 0
Page 53 sur 73
c EDF SA 2013
EDF R&D
H-P74-2012-02004-EN
Reference manual
Sisyphe 6.3
and
f or j = jmax
Version 1.0
Esn+1 ( jmax ) = Esn ( jmax ) + dEs ( jmax )
Msn+1 ( jmax ) = Esn+1 ( jmax )Cs ( jmax )
Underneath layers (from jmax + 1 to NCOUCH TASS) are not eroded and keep their thickness
and mass.
14
Mass conservation
The subroutine suspension bilan coh.f calculates at each time step the mass of sediments in
the computational domain and ensures that the sum of the different components is compensated by
the fluxes at the liquid boundaries. The mass of sediment in suspension M1 is given by :
y
x
M1 = ρs
c( x, y, z)δv = ρs
C ( x, y)hδs
S
with C the depth-averaged volume concentration of sediment in suspension, h the water depth and
S the surface of the computational domain.
In finite elements, each 2D variable is decomposed on basis function φi :
x
M1 = ρs ∑ Ci hi
φi δs = ρs ∑ Ci hi Si
i =1,N poin
i =1,N poin
with Si the integral of basis function. The mass of sediment in the bed M2 is expressed as :
y
x Z Z f
M2 =
Cs δv =
Cs δz δs
Zr
with Cs the sediment bed mass concentration, Z f the bed level and Zr the non erodible bed level.
After discretization of the bed layers into layers of variable thickness Es and concentration Cs
!
!
x
M2 =
∑ Ms ( j) δs = ∑
∑ Ms ( j ) Si
i =1,N poin
j=1,Nomblay
j=1,Nomblay
i
where Msi is the total mass per surface area at node i.
Keywords
– MASS-BALANCE (= NON, default option)
15
Consolidation algorithm
In Sisyphe , the effect of consolidation is modeled through physical as well as empirically-based
models. In a so-called iso-pycnal scheme, the bed is represented by a number of layers of increasing
concentration. The concentration of each layer is fixed, while their thickness varies in time as the
bed undergoes erosion/deposition/consolidation. This scheme is used in a semi-empirical model,
originally developed by Villaret and Walther [58]. A new scheme developed by Van [35] is based on
the Gibson’s theory. An other third scheme has been implemented which is similar to the Gibson
consolidation model developed in Telemac-3d by LeNormand et al. [32]. In this third model, the
bed is discretised in layers of given thicknesses and time-varying concentrations. This third model
has been shown to be unstable and CPU time-consuming (cf. [35]).
Accessibility : Interne : EDF SA
Page 54 sur 73
c EDF SA 2013
EDF R&D
H-P74-2012-02004-EN
Reference manual
Sisyphe 6.3
Version 1.0
Consolidation theroy
1. Most consolidation theories describe the primary consolidation since they are based on the
theory by Terzaghi (1923) [50].
2. The theory by Terzaghi relates the deformation of the bed to the permeability and the effective (or solid) stress which is defined as the difference between the total stress and the pore
pressure (principle of effective stress [50]). It corresponds to the stress that is transmitted directly through the contacts between solid particles. It is also called osmotic pressure by some
authors. Been and Sills (1981) evidenced the progressive increase of effective stress during the
consolidation process.
3. The theory by Terzaghi was originally formulated for infinitesimal strains so that permeability
and compressibility could be assumed constant. The theory was extended for large deformations by Gibson et al. [20, 21] as pointed out by Been and Sills [6].
As an illustration of the process in the water column, Figure 3 proposes a schematic representation of both, the sedimentation and the consolidation. In the right side of Figure 3, we present
the validity range of the Kynch theory of sedimentation and of the Gibson theory of large strain
consolidation. Both theories are unified by Toorman [47, 48] which pointed out the fact that Gibson
model can also describe the sedimentation of particle depending on the choice of closure equations.
The general model (Equations 56, 57 and 58) proposed by Gibson et al. [20, 21] represents the
different stages of consolidation. This equation is based on a two-phase approach by considering
continuity and motion equations for both fluid and solid phase to obtain the general equation that
reads in material coordinate ζ which represents the volume of solids :
∂e
+
∂t
ρs − ρ f
ρf
!
d
de
k
1+e
∂e
∂
+
∂ζ
∂ζ
k
dσ0 ∂e
gρ f (1 + e) de ∂ζ
!
= 0,
(56)
where e is the void ratio, k the hydraulic permeability (in m/s) and σ0 the effective (or solid) stress.
The Gibson equation can be also written in Eulerian framework as :
∂e
+ (1 + e )2
∂t
ρs − ρ f
ρf
!
∂
∂z
k
(1 + e )2
+
(1 + e )2 ∂
gρ f ∂z
k ∂σ0
1 + e ∂z
=0
(57)
This equation is equivalent to :
∂φ
∂
−
∂t
∂z
"
k ∂σ0
k ( s − 1) φ +
γ f ∂z
! #
φ =0
(58)
where φ stands for the sediment volume concentration, k the hydraulic permeability, s the density
ratio between sediment and fluid (=ρs /ρ f ), γ f the unit weight of fluid (=gρ f , g being the acceleration
of gravity), z the vertical coordinate (positive upward), and t the time.
The equation of Gibson has been widely used in various numerical consolidation models [3, 6,
47, 48] as well as compared with the experimental results. It has been implemented in Telemac-3d
by LeNormand et al. [32] or coupled with Telemac-2d by Thiebot [45]. Their method of resolution
has been adapted to Sisyphe by Van [35].
The main difficulty in using the Gibson model is related to the choice for closing the problem.
Two closure equations, for the permeability k and for the effective stress σ0 , is indeed required to
obtain the time evolution of vertical concentration profiles. However, the formulation of these closure
equations remains an open problem and a shared protocol to determine their parameter values are
still lacking as reported by Toorman [47, 48], Bartholomeeusen et al. [3] and Van [35].
Accessibility : Interne : EDF SA
Page 55 sur 73
c EDF SA 2013
EDF R&D
Reference manual
Sisyphe 6.3
H-P74-2012-02004-EN
Version 1.0
Figure 4 – Diagram of different processes involved in the settling transport (left : non-cohesive,
right : cohesive).
15.1
Multi-layer empirical algorithm
Time evolution
The consolidation effect is reproduced by assuming that the vertical flux of sediment from layer
j to underneath layer j + 1 is proportional to the mass of sediments, Ms (kg/m2), contained in the
layer j.
dMs ( j)
= a i Ms ( j )
dt
(59)
The transfer mass coefficients a j (in s−1 ) are specified in the steering file (MASS TRANSFER PER
LAYER). They correspond to a characteristic timescale to transfer mass from one layer to another.
This multilayer consolidation model is not related to Gibson equation so that it should be considered as empirical. Despite its apparent simplicity, it can qualitatively reproduce an increase of mud
bed deposit with time, while ensuring mass conservation (transfer coefficient of the last bottom layer
is zero). The model results are highly sensitive to the specified values of the mass transfer coefficients
ai which are however difficult to calibrate. The main advantage of this iso-pycnal formulation relies
in the fact that no vertical grid is required to compute the concentration evolution of the different
layers.
Discretization
The vertical resolution of semi-empirical equation (59) is based on the multi-layer iso-pycnal
model with fixed concentrations. The consolidation is then reproduced by mass transfer between
layers of the model. The transfer coefficients are fixed for each layer.
The set of mass-transfer coefficients a( j) in (s−1 ) are selected by calibration, in order to find best
agreement of time-varying concentration profiles between model and experiment. Physically, they
represent the inverse of the residence time per layer. The values a( j) are found to decrease from
top to bottom, as the time scale of consolidation (or residence time) increases. The mass transfer of
the last layer is set to zero, in order to insure no mass loss at the rigid bed level (impermeability
condition).
Accessibility : Interne : EDF SA
Page 56 sur 73
c EDF SA 2013
EDF R&D
Reference manual
Sisyphe 6.3
H-P74-2012-02004-EN
Version 1.0
Keywords
–
–
–
–
COHESIVE SEDIMENT = YES (= NO, default option)
MUD CONSOLIDATION = YES (= NO, default option)
CONSOLIDATION MODEL = 1, default option
NUMBER OF LAYERS OF THE CONSOLIDATION MODEL (NOMBLAY=10 by default, maximum value
is 20 )
– MUD CONCENTRATION PER LAYER (in kg/m3 , E = 50.; 100.; ... by default)
– MASS TRANSFER PER LAYER (= 5.0E − 5, ..., 0., by default)
15.2
15.2.1
Multi-layer iso-pycnal Gibson’s model
Theoretical background
The previous equation presents the advantage of not using a vertical mesh. However it is not
physically based model since it neglects the Gibson theory. Improvement of this category of model
was proposed by Sanchez [41] and Thiebot [45]. In these formulations, the sediment flux from one
layer to the other is given not in term of empirical mass transfer coefficient (ai ) but in term of
theoretical flux which is given by the Gibson theory.
15.2.2
Numerical discretization
The model originally developed by Sanchez [41] and Thiebot [45] is a 1DV sedimentationconsolidation multi-layer model, based on an original technique to solve Gibson equation. The advantage of this representation (Equation 60) if compared with previous one (Equation 59) relies on
the right hand side term which relates the mass of sediments, Msi (kg), to the net flux entering or
leaving the layer. This equation enables therefore to consider the flux of sedimentation and consolidation as provided by the Gibson theory.
dMsi
= ( Fi (t) − Fi+1 (t))∆tπr2
dt
(60)
In Equation 60, a circular shape (of radius r) is assumed for the surface. Equivalently, Equation 61
is expressed in term of layer thickness Esi (m), recalling Msi = Csi πr2 Esi. . This last expression
becomes independent on the settling column geometry.
dEsi
F (t) − Fi+1 (t)
= i
dt
Csi
(61)
The concentration of the different bed layers Cs(i ) is fixed. As the sedimentation and consolidation progress, the sediment is transferred to the more concentrated layers, and the thickness of these
layers increase as well.
Accessibility : Interne : EDF SA
Page 57 sur 73
c EDF SA 2013
EDF R&D
H-P74-2012-02004-EN
Reference manual
Sisyphe 6.3
Version 1.0
The mass conservation is ensured by requiring at each moment, in each layer, the equality between the mass contained in a layer at time t + ∆t and the mass present in this layer at time t in which
the outgoing mass was removed and the incoming mass was added (means the mass that crossed
the upper and lower sections respectively during the time ∆t). The outgoing and incoming masses
are taken into account by sediment flux noted Fi (t).
Fi (t) =
(Vs,i (t) − Vs,i−1 (t))Ci−1 Ci
Ci−1 − Ci
(62)
where Vs,i is the falling velocity of the layer i, and can be defined as :
Vs,i (Ci ) =
k(Ci )Ci
k(Ci )Ci
1
1
−
ρs
ρf
!
1
1
−
ρs
ρf
!
if Ci ≤ Cgel
,
+ k(Ci )
σ0 (Ci−1 ) − σ0 (Ci )
1
( Epi−1 (t) + Epi (t))
2
(63)
,
otherwise
where k is the permeability, σ0 is the effective stress, Cgel the transition concentration between
sedimentation and consolidation schemes (see Camenen & Pham Van Bang [9]).
15.3
Vertical grid Gibson’s model
Vertical-grid models propose a natural resolution of the Gibson equation as a vertical grid is used
to compute the concentration at each point. This category of model uses common techniques (finite
difference or finite volume or finite element methods) for resolving partial differential equations.
As they use a vertical grid, the connection with Sisyphe which is depth integrated model is not
straightforward.
In the new version of Sisyphe , such a category of model is also available. Strictly speaking,
the original Gibson model (in material coordinate) used in Telemac-3d which was developed by
Lenormant et al. [32] is connected to Sisyphe .
The finite difference method is used. The implicit scheme leads to a tridiagonal matrix that is
solved by using a classical double sweep algorithm. More details are given in Telemac-3d user
Accessibility : Interne : EDF SA
Page 58 sur 73
c EDF SA 2013
EDF R&D
Reference manual
Sisyphe 6.3
H-P74-2012-02004-EN
Version 1.0
guide, LeNormand et al. [32] and Van [35].
16
Closure equations for permeability and effective stress
All of the Gibson based models, i.e. Multi-layer iso-pycknal Gibson’s model (section 15.2) and
vertical grid Gibson model (section 15.3), requires two closure equations. The Multi-layer empirical
model (section 15.1) only needs closure on ai, the mass transfer coefficients.
In this section, some available empirical formula are detailed. In general, a decreasing (or increasing) function of solid concentration is considered for permeability K (or effective solid stress). These
expressions are logarithmic or exponential or power law functions (see Bartholomeeusen et al. [3]).
16.0.1
Closure equations
– Bed permeability k
Bartholomeeusen et al. [3] introduced typical functions, in the form of either power or exponential, to relate the permeability k with the void ratio e :
A 1 e A2
k=
A φ − A2
1 s
exp(− A1 + A2 e)
The value of these coefficients A1 and A2 depend on grain size distribution, organic content,
activity and pore size distribution.
– Effective stress σ0
The similar way can be opted in the determination of σ0 (e) (σ0 (C ), σ0 (φ)), which gives :
(
0
e = − B1 σ B2 + B3
B
σ0 = B1 φs 2
or
e = B1 (σ0 + B2 )− B3
σ0 = exp( B4 + B5 e)
According to Winterwerp & van Kesteren (2004), the validity range of the power-type functions (e.g. k = A1 e A2 , σ0 = B1 φsB2 ) is much larger than that of the exponential-type functions.
Moreover, there is a physical insight in the formulation of the power-type functions. Therefore,
it is recommended to use power-type functions. However, two different power functions for
permeability may be necessary to represent the two separate processes : sedimentation and
consolidation.
16.0.2
Determination of parameters
No standard is found in the literature to recommend a specific methodology for calibrating the
empirical functions for both permeability and effective stress. Most of study reported some fitting
exercice on settling curve, i.e. the position of supernatant/suspension interface. They considered
mostly the least square technique for the adjustment to experimental results. However as stated
by Toorman [48] a more relevant adjustment is offered when concentration profiles are available
from Gamma-ray techniques (see Been and Sills [6] or Bartholomeeusen et al. [3]), X-ray technique
(Villaret et al. [59]) or MRI techniques (Pham Van Bang et al. [38]). In such a situation, the density or
concentration profile are adjusted on the measurement for different time. A specific procedure has
been recently proposed by Thiebot et al. [46], also based on least square method.
The user should consider the previous method (adjustment on settling curve or on concentration
profile) as conventional even through no shared procedure has been internationality recognized as
the best one (as previously discussed in 10.3). An recent alternative was successfully applied to
Gironde mud (Van [35]).
Accessibility : Interne : EDF SA
Page 59 sur 73
c EDF SA 2013
EDF R&D
16.0.3
H-P74-2012-02004-EN
Reference manual
Sisyphe 6.3
Version 1.0
Space-time based method to determine parameters of closure equations
The new methodology to determine parameters is based on space-time analysis of data. It requires data on time evolution of concentration (vertical) profile. If such an information is not available, the proposed method becomes useless. If only the supernatant/suspension interface position
is detailed by data, the user should consider the classical fitting method.
If the space-time resolved data is available, this new methodology could be applied easily. The
method is based on a theoretical framework which is an advantage to help during the calibration
procedure. From time evolution of concentration profiles which are provided by non-intrusive techniques (here the X-ray technique), the procedure uses self-similar analytical solutions to determine
the closure equation relative to the convection (or sedimentation) part and to the diffusion (or consolidation) part of Gibson equation.
∂
∂φ
∂
∂φ
D (φ)
=0
−
[V ( φ ) φ ] −
∂t
∂z
∂z
∂z
where V (φ) = K (s − 1)φ and D (φ) = Kφ/γ f dσ0 /dφ in order to match with Gibson equation 58 in
Eulerian coordinate.
Considering a separation regime between sedimentation and consolidation (or between convection and diffusion), self-similar solutions are obtained for each regime. The separation between both
problems is justified by the fact that effective stress should be zero for a suspension since there is no
direct contact between particles.
16.0.4
Determination of closure equation for the sedimentation
If the concentration of the suspension is lower than a given threshold (the so-called gelling point
for cohesive sediment), the inter-particle contacts are negligible, i.e. there is no solid (or effective)
stress (Camenen & Pham Van Bang [9]). In such a situation the Gibson equation (58) is simply
reduced to the equation of Kynch :
∂φ
∂
(64)
−
[V ( φ ) φ ] = 0
∂t
∂z
where φ = C/ρs is the volume fraction of solids, V (φ) is the settling velocity of the suspension at
concentration φ that is equal to K (s − 1)φ.
Considering the self-similar variable ζ = z/t and similarity solution U, i.e. φ(z, t) = U (ζ ), equation 64 leads to :
df
−ζ
dU
dU
=0
dζ
where f is the solid (or sedimentation) flux, which is equal to −V (φ)φ.
The method is equivalent to the so-called method of characteristics (Leveque [33]). The isoconcentration pattern presents different straight lines in the space-time (z − t) plot. The slopes of
iso-concentration straight lines in the z − t plane are measured to obtain the first derivative of the
solid flux. The sedimentation flux proposed for the Gironde mud (Villaret et al. [59] ; Van [35]) is
given by equations 65 or 66. The first derivative of this closure equation is straightforward [59] : the
determination and validation of its parameters (Vst , φgel , n) has been presented in details in the test
case of settling column of Gironde mud :
φ
f (φ) = Vst (1 − φ) 1 −
φgel
Accessibility : Interne : EDF SA
!n
Page 60 sur 73
φ
for
φ < φgel
(65)
c EDF SA 2013
EDF R&D
H-P74-2012-02004-EN
Reference manual
Sisyphe 6.3
Version 1.0
and
C
C
f (C ) = Vst (1 − ) 1 −
ρs
Cgel
!n
C
ρs
for
C < Cgel
(66)
where Vst is the Stokes velocity of an equivalent sphere, φgel (Cgel ) is the gelling concentration,
n is an exponent. Both parameters have physical meanings in terms of rheology (Pham Van Bang
et al. [38]). Indeed, φgel (Cgel ) is the concentration value from which the effective viscosity of the
suspension diverges. And n is a parameter describing the transition from suspension to a structured
bed. It is worth noting that the so-called Richardson & Zaki empirical model is recovered by setting
φgel = 1.
16.0.5
Determination of closure equation for the consolidation
Still considering separation regime between sedimentation (convective problem) and consolidation (diffusion problem), for cohesive sediment and concentration larger than the gelling point, the
effective solid stress build up. The diffusion term is no longer negligible and becomes the leading
term in the Gibson equation. For concentration larger than the gel point, the convective part is annihilated in the proposed formulation so that only the diffusion term remains. As this term takes
origin from the competition between the seepage flow through a porous media and the effective
stress of the solid skeleton (Camenen & Pham Van Bang [9]), the diffusion is expressed by an expression with the permability K and the effective stress (dσ0 /dφ). The problem is now related to the
non-linearity of the diffusion :
∂
∂φ
∂φ
−
D (φ)
=0
∂t
∂z
∂z
(67)
where D (φ) is the non linear diffusion term, equal to Kφ(dσ0 /dφ)/γ f in order to match with
(58).
The power law is assumed for the diffusion coefficient. Indeed, if we consider a power law for the
permeability and a power law for the effective stress, the diffusion coefficient will also be a power
function of concentration. Here the possible time dependence of the effective stress is investigated
so that the secondary consolidation is also taken into consideration (cf. 10.3). As a result, the non
linear diffusion coefficient is assumed to depend on both concentration and time by an empirical
power law, i.e. D (φ) = D0 φ a tb . Here the time dependence of the consolidation is introduced to
mimic the thixotropic behaviour of mud. Introducing as a self-similar variable, χ = z/tθ with θ =
(1 + b)/(2 + a) and the similarity solution h, i.e. φ(z, t) = h(χ)/tθ in equation (67) leads to :
∂φ
∂
a b ∂φ
−
D0 φ t
=0
∂t
∂z
∂z
d
dh
a
h (χ)
− θχh(χ) = 0
dχ
dχ
→
The similarity solution h is obtained after integration of the previous ODE equation (see Van [35]).
The parameter M is a constant whose value corresponds to the total mass of sediment in the system.
1/a
a (1 + b ) χ2
M
−
2(2 + a )
h(χ) =
0
for
for
"
#
−2M(2 + a) 1/2
χ ∈ 0,
a (1 + b )
−2M(2 + a) 1/2
χ≥
a (1 + b )
In order to find out the closure equation for effective stress, from the above equations, it is
needed to determine the two variables : a, b. This procedure will be presented in detail in test case
Tassement 2, with a given experimental result of settling column.
Accessibility : Interne : EDF SA
Page 61 sur 73
c EDF SA 2013
EDF R&D
Reference manual
Sisyphe 6.3
H-P74-2012-02004-EN
Version 1.0
Keywords
–
–
–
–
–
–
–
COHESIVE SEDIMENT = YES (= NO, default option)
MUD CONSOLIDATION = YES (= NO, default option)
CONSOLIDATION MODEL = 2
GEL CONCENTRATION = CONC GEL
PERMEABILITY COEFFICIENT = COEF N
NUMBER OF LAYERS OF THE CONSOLIDATION MODEL = NOMBLAY (< 20)
MUD CONCENTRATION PER LAYER
Accessibility : Interne : EDF SA
Page 62 sur 73
c EDF SA 2013
EDF R&D
17
Reference manual
Sisyphe 6.3
H-P74-2012-02004-EN
Version 1.0
Mixed sediments processes
17.1
Mixte sediment bed composition
A sand-mud mixture is composed of two classes of bed sediments : The first class (noted 1)
is non-cohesive and represented by its grain diameter D1 and can be transported as bed-load or
suspended load. The settling velocity ws1 is a function of the relative sediment density (s = 1.65)
and grain diameter D1 . So far we assume only suspended load, which implies that the model is
devoted to mixtures of fine sand grains and mud.
The second class (noted 2) is cohesive, grain diameter D2 less than 60 mm. The settling velocity
ws2 is a function of flocs properties which differs from the individual cohesive particles, and needs
to be specified. The mass concentration is the mass of sediments per volume occupied by the mud
component (solid particles + water).
17.2
Erosion/deposition of sand-mud sediment bed
The erosion/deposition fluxes depend on the mass percentage % of each class in the surface layer
(Panagiotopoulos et al., 1997). If the mass percent of the mud class is greater than 50%, the bed is
considered as pure cohesive and the erosion/deposition laws follow the Partheniades classical law
(see below for the calculation of the bed shear strength in a sand/mud mixture).
If the mass fraction of mud percentage is less than 30%, the bed is considered as non-cohesive,
and has little effect on the bed shear strength. If the mass fraction of mud percent is greater than
30%, the bed is considered as cohesive. The rate of erosion is based on the Partheniades erosion law,
where the critical bed shear strength of the mud class depends on the consolidation state. Different
regimes can be considered depending on the volume content of each class, as pictured below :
Figure 5 – (a) Sand regime ( f 2 < 30%) ; (b) Intermediate ; (c) Mud-dominated regime ( f 2 > 50%).
17.3
Critical erosion shear stress of sand mud mixture
The mean bed shear strength (τce ) of the sand/mud mixtures depends on the mass percentage %
of mud f 2 at the surface and depends on the layer. The proposed framework is based on the work
of Waeles (2005). The bed shear strength is calculated in the subroutine suspension flux mixte.f
Accessibility : Interne : EDF SA
Page 63 sur 73
c EDF SA 2013
EDF R&D
Reference manual
Sisyphe 6.3
H-P74-2012-02004-EN
Version 1.0
(TOCE MIXTE, for τce ). For each layer (K), the expression of TOCE MIXTE depends of the mass percent
of mud in the layer :
– For f 2 < 30%, τce = τce (1) for pure sand, based on
τce = τce (1) = Θc gd50 (s − 1)
(68)
– For f 2 > 50%, τce = τce (2) for pure mud as a function of bed layer concentration
τce = τce (2)
(69)
– In the intermediate range (30% < f 2 < 50%) :
τce = f 1 τce (1) + f 2 τce (2)
(70)
Figure 6 – Mixed sediment bed structure.
The potential erosion flux depends on the composition of the bed layer. The first step is to determine a mean erosion rate per layer E (FLUER LOC) :
– For f 2 < 30% (sand dominant) : we use the equilibrium concentration concept :
E1 = ws1 Ceq , for τb τce
(71)
E=
0,
otherwise
– For f 2 > 50% (mud dominant) : the Partheniades erosion law can be applied :
2
E2 = M ττce0
− 1 for τb > τce
E=
0
otherwise
– For 30% < f 2 < 50% (intermediate range) : the mean erosion rate is an average value of the
erosion classes of both classes :
E = f 1 E1 + f 2 E2
Accessibility : Interne : EDF SA
Page 64 sur 73
c EDF SA 2013
EDF R&D
17.4
H-P74-2012-02004-EN
Reference manual
Sisyphe 6.3
Version 1.0
Iterative procedure
The next step is to calculate the mean erosion flux for both classes of bed material. The underlying
assumption is that the eroded depth is the same for both classes. This is done iteratively, by checking
layer by layer (starting from the top) how much sediment there is in the layer and how much time
it would take to empty that layer entirely : This procedure is necessary because the bed is stratified
and the eroded depth is unknown.
At the start of time step (t0 ), the first layer whose thickness is non-zero is eroded, at a given rate
e( j). For each erodible layer j with E( j) > 0, the time interval dt( j) to erode completely layer j is
estimated based on the mass of the layer :
dt( j) =
dEs ( j)Cs ( j)
dMs ( j)
=
,
ρs E( j)
ρs E( j)
where the potential erosion rate E( j) decreases and tends to zero as j increases. The last (non-empty)
layer (noted jmax) to be eroded is obtained by emptying all successive top layers until :
jmax −1
∑
jmax
dt( j) < Dt <
j =1
∑ dt( j).
j =1
j
For jmax , the time left to erode the last layer is Dtmax = Dt − ∑ jmax
=1
mass potentially eroded during this process is :
−1
Dt( j). The total (mud and sand)
jmax −1
Ms =
∑
Ms( j) + ρs E( jmax ) Dtmax .
j =1
The mean mass erosion flux of the sand-mud mixture Es (Kg/m/s) is therefore estimated as
Es =
Ms
.
Dt
Here we assume that both mud and sand are eroded simultaneously to the same depth, therefore
the eroded volume flux of each phases is differentiated by their mass percent :
E1 =
17.5
17.5.1
Ms1
,
ρs Dt
E2 =
Ms2
.
ρs Dt
Definition scheme of sand-mud mixture in Sisyphe
Assumptions
In this framework, the sand particles are assumed to be fully embedded in the cohesive framework. So this is valid for cohesive dominant regime. The density of individual primary or sand
particles is assumed to be constant and equal to ρs = 2650 kg/m3
Accessibility : Interne : EDF SA
Page 65 sur 73
c EDF SA 2013
EDF R&D
Reference manual
Sisyphe 6.3
H-P74-2012-02004-EN
Version 1.0
Figure 7 – Critical erosion shear stress of sand-mud mixture (from Waeles [62]).
17.5.2
Bed structure
The sand-mud mixture is characterized by the volume percent of each class :
Pi =
Vi
Es
= i,
Vt
Et
where Vi is the volume occupied by the sand or mud particles i = 1, 2, Vt = V1 + V2 , is the total
volume, Esi = pi Es is the volume per unit surface area of phase i and corresponds to the vertical
thickness of each phase (voids are not included), and Et is the total volume (sand + mud) per unit
surface area Et = E1 + E2 . It is also useful to define the mass percent of each class, for sand f 1 and
for mud f 2 :
ρs V1
C2 V2
; f2 =
,
f1 =
Mt
Mt
where C2 is the mud density, which varies with distance for the bed surface and Mt the total mass
of the sand-mud mixture :
Mt = C2 V2 + ρs V1 .
Once sediment particles have been put in suspension, they are transported independently, by
solving for each class a transport equation for the volume solid concentrations, defined as :
Ci =
Vi
,
Vt
with Vt = V1 + V2 + Vwater . The volume and density of suspended cohesive flocs vary as a result of
flocculation, so what is conserved here is the volume of primary constitutive particles of density ρs .
17.6
Suspended sediment transport equation
Both sand and mud are treated separately and independently. We solve here the volume 2D (
depth integrated) concentrations of the primary particles, assuming constant density ρs :
∂Ci
∂C
∂C
∂
∂C
∂
∂C
( E − Di )
+ Uconv i + Vconv i =
es i +
es i
+ i
∂t
∂x
∂y
∂x
∂x
∂y
∂y
h
Accessibility : Interne : EDF SA
Page 66 sur 73
c EDF SA 2013
EDF R&D
Reference manual
Sisyphe 6.3
H-P74-2012-02004-EN
Version 1.0
– For sand : D1 = ws1 CF1 , with F1 the ratio between the near bed concentration and depthaveraged concentration.
– For mud :
(
ws2 C 1 − ττ , for τ ≤ τd
d
D2 =
0
otherwise
with τd the critical bed shear stress for deposition.
Figure 8 – Mixed sediment in suspension.
Attention
In the subroutines, the definitions of Sisyphe variables are : horizontal point index : I ; layer
index : K ; number of classes of bed material (NSICLA = 2) ; J = 1 : sand ; J = 2 : mud ; grain
sizes : FDM(J) and number of bed layers : COUCH TASS.
Attention
For pure mud, the bed shear strength increases as the concentration of the mud increases
while consolidation is taking place. Mud concentration per layer : CONC VASE (K) (Kg/m3) ;
critical erosion shear stress per layer : TOCE VASE(K) (N/m2).
Accessibility : Interne : EDF SA
Page 67 sur 73
c EDF SA 2013
EDF R&D
Reference manual
Sisyphe 6.3
H-P74-2012-02004-EN
Version 1.0
Keywords
To enable the mixed sediment transport in Sisyphe , the following keywords should be specified in the steering file :
– MIXED SEDIMENT = YES (= NO, default option)
– COHESIVE SEDIMENTS = YES (= NO, default option. In Sisyphe , the first sediment class is
always sand and the second is mud)
– NUMBER OF SIZE-CLASSES OF BED MATERIAL = 2 (NSICLA)
– INITIAL FRACTION FOR PARTICULAR SIZE CLASS = f 1 ; f 2 (example : = 0.80; 0.20)
– SEDIMENT DIAMETERS = D1, D2 (for non-cohesive sediments D1>6.0E-5, for cohesive sediments D1<6.0E-5, example with 300 µm for sand and 60 µm for mud : = 0.0003; 0.00006)
– SEDIMENT DENSITY = 2650 (the sediment density is only for sand, while the mud density is
defined in the keyword MUD CONCENTRATION PER LAYER)
– SETTLING VELOCITIES = ws1 , ws2 (example = 4.E − 02; 1.E − 02, if this one is active, only
one value is taken for calculation. It is corrected later in the subroutine)
– INITIAL SUSPENSION CONCENTRATIONS (example = 0.0002; 0.0002)
– NUMBER OF LAYERS OF THE CONSOLIDATION MODEL (NOMBLAY = 10, default option)
– MUD CONCENTRATION PER LAYER (in kg/m3 , CONC VASE= 50.; 100.;..., default option)
– CRITICAL EROSION SHEAR STRESS OF THE MUD (example = 0.5; 10.0)
17.7
Fortran subroutines for mixed sediment transport
– INIT SEDIMENT : This subroutine initializes the bed composition (by calling INIT MIXTE for
mixed sediment transport) and calculates the settling velocity, if not present in the steering
file Moreover, shields parameter and CRITICAL EROSION SHEAR STRESS OF THE SAND known
as TOCE SABLE is calculated in this subroutine as well.
– INIT MIXTE : This subroutine calls INIT COMPO COH to initialize the composition of the sediment
bed and also check the thickness of each bed layer. The initial masses of sand and mud are
also calculated here. It can be seen that XMVS is the sand density defined in the steering file
while CONC(I,J) is the mud density of each bed layer defined in the steering file to known as
MUD CONCENTRATION PER LAYER. For this subroutine, the density of the sand is assumed to be
constant (XMVS) and the following variables are defined :
– Mass per unit surface area of the sand : MS SABLE(I,K) (Kg/m2)
– Mass per unit surface of the mud : MS VASE(I,K) (Kg/m2)
– Total Mass per layer : MS TOTAL(I,K) = MS SABLE(I,K) + MS VASE(I,K)
– Mass percent of each phase per layer : Fs(I,K,1) = MS SABLE(I,K)/MS TOTAL(I,K) and
Fs(I,K,2) = MS VASE(I,K)/MS TOTAL(I,K)
– Initial mass of the sand bed : MASS0
– Initial mass of the mud bed : MASV0
– INIT COMPO COH : This subroutine initializes fraction distribution of sand and mud. The thickness of each bed layer is also defined here. For this simulation, there are two bed layers defined
in this subroutine as the sum of the thickness of mud per layer (EPAI VASE) and the thickness of
sand per layer (EPAI SABLE). The modification is made for reading sand fraction from data file.
We do not take into account the porosity of sand, because we assume the void to be entirely
filled up by the fine particles. This may be wrong when the mud percent % is lower than the
void ratio (in case of bedload layer for example). Up to 6.3, the present the mixte sediment
formulation is only valid assuming the sand grains to be imbedded in the cohesive frame.
– Release 6.3 : the thickness of each sand layer corresponds here to the volume occupied by
the sand particles, per unit surface area.
– Bed layer thickness for sand : EPAI SABLE(K) (m)
– Bed layer thickness for mud (m) : EPAI VASE(K) (m)
– Bed layer thickness (mud and sand) : ES(I,K) = EPAI SABLE (K) + EPAI VASE (K)
– Volume percent of each phase : AVAIL(I,K,1) = EPAI SABLE(K)/ES(I,K) and AVAIL(I,K,2)
= EPAI VASE(K)/ES(I,K)
– Future release : introduce a test on AVAIL (only correct if AVAIL(I,K,2) > void ratio of ≈ 0.4
Accessibility : Interne : EDF SA
Page 68 sur 73
c EDF SA 2013
EDF R&D
Reference manual
Sisyphe 6.3
H-P74-2012-02004-EN
Version 1.0
(=1-XKV))
– SUSPENSION COMPUTATION : This is the main subroutine for the computation of the suspension
load and the bed elevation. It has the following steps :
1. COMPUTES THE REFERENCE ELEVATION --> ZREF
2. ADVECTION VELOCITY --> UCONV, VCONV, TAKING INTO ACCOUNT THE VERTICAL PROFILE
OF CONCENTRATIONS AND VELOCITIES
3. EROSION FLUX : FLUER In this part, the subroutine SUSPENSION FLUX MIXTE is called to
compute the suspension for the sand and mud mixture.
4. DEPOSITION FLUX : FLUDPT =WC*T2 This subroutine calls SUSPENSION DEPOT to compute
the sand and mud deposition.
5. DIFFIN A SPECIFIC TREATMENT IS DONE IF THE ADVECTION METHOD IS THE CHARACTERISTICS:
FREE OUTPUTS ARE TREATED LIKE DIRICHLET. THIS SPECIFIC TREATMENT IS CANCELLED
HERE BY SENDING A MODIFIED VALUE FOR RESOL.
6. BOUNDARY CONDITIONS : CBORCONC V The keyword IMP INFLOW C is found here for controlling the calculation of equilibrium concentration at open boundaries.
7. SOLVING TRANSPORT EQUATION IF METHOD OF CHARACTERISTICS
8. SOURCE AND SINK TERMS (IMPLICIT SOURCE TERM FOR THE DEPOSITION T9 AND EXPLICIT
SOURCE TERM WITHOUT PUNCTUAL SOURCES T11)
9. ADVECTION-DISPERSION STEP CONFIGURATION OF ADVECTION
10. BED EVOLUTION DUE TO NET EROSION/DEPOSITION FLUX The subroutine SUSPENSION EVOL
is called from here to calculate the bed evolution.
Accessibility : Interne : EDF SA
Page 69 sur 73
c EDF SA 2013
EDF R&D
Reference manual
Sisyphe 6.3
H-P74-2012-02004-EN
Version 1.0
Références
[1] Ashida K. and Mishiue M., 1973. Studies on bed load transport rate in alluvial streams, Trans JSCE,
vol.4.
[2] Bailard J., 1981. An energetics total load transport model for a plane sloping beach. Journal of
Geophysical Research, 86, C11, 10938-10954.
[3] Bartholomeeusen, G., Sills, G.C., Znidarcic, D., van Kesteren, W., Merckelbach, L.M., Pyke, R.,
Carrier, W.D., Lin, H., Penumadu, D., Winterwerp, H., Masala, S. and Chan, D., 2002, Sidere :
numerical prediction of large strain consolidation, Géotechnique 52, No. 9, 639-648.
[4] Begnudelli L., Valiani A. and Sanders B.F., 2010. A balanced treatment of secondary currents,
turbulence and dispersion in a depth-integrated hydrodynamic and bed deformation model for
channel bends. Advances in Water Resources, 33(1).
[5] Belleudy P., 2000. Numerical simulation of sediment mixture deposition, part 1 : analysis of a
flume experiment, Journal of Hydraulic Research, Vol.38, N◦ 6.
[6] Been, K. and Sills, G.C., 1981, Self weight consolidation of Soft Soils : an Experimental and Theoretical
Study. Geotechnique, Volume 31, No.4, pp. 519-535.
[7] Bijker E.W., 1968. Mechanics of sediment transport by the combination of waves and current. In
Design and Reliability of Coastal Structures, 23rd Int. Conf. on Coastal Engineering, 147-173.
[8] Camenen B., 2002. Modélisation numérique du transport sédimentaire sur une plage sableuse. Ph.D
thesis, Université Joseph Fourier, Grenoble 1.
[9] Camenen B. & Pham Van Bang D., 2011. Modelling the settling of suspended sediments for concentrations close to the gelling concentration, Cont. Shelf Res., 31, S106-S116.
[10] Blue Kenue : Tutorial Manual for the Mesher in BlueKenue. Canadian Hydraulics Centre (CHC)
of the National Research Council, 2013.
[11] Celik I., Rodi W., 1988. Modelling Suspended Sediment Transport in non- equilibrium situations, Journal of Hydraulic Eng., Vol. 114, N◦ 10.
[12] Chollet J.P. and Cunge J.A., 1980. New interpretation of some headlooss - flow velocity relationship for deformable bed, J. Hydr. Eng., 17 (1).
[13] Davies A.G. and Villaret C., 2004. Modelling the effect of wave-induced ripples on littoral sand transport. EDF-LNHE Report HP-75/03/029/A.
[14] Dibajnia M., Watanabe A., 1992. Sheet flow under non-linear waves and currents. Proc. of the Intl.
Conf. on Coast Eng., 2015-2029.
[15] Die-Moran A., 2012. Physical and Numerical modelling investigation of induced bank erosion as a
sediment transport restoration strategy for trained rivers. The case of the Old Rhine (France). Université
Paris-Est.
[16] Egiazaroff I.V., 1965 : Calculation of non-uniform sediment concentrations, J. of Hydr. Div. ASCE,
vol. 91, N◦ 4, pp. 225-248.
[17] Einstein H.A. (1950). The bed load function for sediment transportation in open channel flow. US Dep.
Of Agriculture, Techn. Bull.
[18] Engelund F. and Hansen E., 1967. A monograph on sediment transport in alluvial streams. Techn.
Univers. Of Denmark, Copenhagen, Denmark.
[19] Engelund F., 1974. Flow and bed topography in channel bends. J. Hydraul. Div. Am. Soc. Civ. Eng.,
100, 1631-1648.
Accessibility : Interne : EDF SA
Page 70 sur 73
c EDF SA 2013
EDF R&D
Reference manual
Sisyphe 6.3
H-P74-2012-02004-EN
Version 1.0
[20] Gibson, R.E., Englund, G. L. and Hussey, M. J. L., 1967, The theory of one dimensional consolidation
of saturated clay, I. Finite Non-Linear Consolidation of Thin Homogeneous Layers. Geotechnique, pp.
261-273.
[21] Gibson, R.E., Schiffman, R.L. and Cargill, K.W., 1981, The theory of One dimensional Consolidation of Saturated Clays, II. Finite Nonlinear Consolidation of Thick Homogeneous Layers. Canadian
Geotechnical Journal, Vol. 18, pp. 280-293.
[22] Gonzales de Linares M., 2002. Graded sediment in Sisyphe. Rapport EDF-LNHE, HP-75/2002/68.
[23] Hirano, M., 1970. On phenomena of river-bed lowering and armouring below reservoirs. In Proc. 14th
Hydraul. Lecture Meeting, 13-14 Feb., Civ. Eng. Ass. Hydr. Committee, Hatsumei Kaikan.
[24] Hervouet J.-M., Hydrodynamics of free surface flows : modelling with finite element methods, John
Wiley & Sons, 2007.
[25] Hervouet J.M., Razafindrakoto E., Villaret C., 2011. Dealing with dry zones in free surface flows, a
new class of advection scheme. IAHR conference, July 2011.
[26] Hunziker R.P., 1995. Fraktionsweiser Geschuebetransport, Ph.D. thesis, Mitteeilungen Nr 138 deer
Versuchanstalt fur Wasserbau, Hydrologie und Glaziologie, ETH Zurich, Switzerland.
[27] Huybrechts N., Villaret C., Hervouet JM, 2010. Comparison between 2D and 3D modelling of sediment transport : application to the dune evolution. Proceedings of the RiverFlow conference.
[28] Huynh-Thanh and Temperville A., 1991. A numerical model of the rough turbulent boundary layer
in combined wave and current interaction. In Sand Transport in Rivers Estuaries and the Sea, eds
RL Soulsby and R. Bettess, pp. 93-100, Balkema, Rotterdam.
[29] Julien P. (2010). Erosion and sedimentation, second edition. Cambridge University Press, Cambridge.
[30] Karim M.F. and Kennedy J.F., 1982. A computer based flow and sediment routing. IIH Report N◦ 250
Modelling for streams and its application to the Missouri River, University of Iowa, Iowa City,
IA.
[31] Koch F.G. and Flokstra C., 1981. Bed level computations for curved alluvial channels, XIXth Congress
of the International Association for Hydraulic Research, New Delhi India.
[32] Lenormant, C., Lepeintre, F., Teisson, C., Malcherek, A., Markofsky, M., and Zielke, W., 1993,
Three dimensional modelling of estuarine processes. In MAST Days and Euromar Market, Project
Reports Volume 1.
[33] Leveque, Randal J., 2002, Finite-Volime Methods for Hyperbolics Problems, Cambridge Texts in
Applied Mathematics, ISBN 0-511-04219-1 eBook.
[34] Mc Anally, W.H., 1999. Aggregation and deposition of fine estuarial sediment Floculation. Ph.D.
Thesis, Univ of Fla, Gainesville, Fla.
[35] Van, L.A., 2012. Modélisation du transport de sédiments mixtes sable-vase et application à la morphodynamique de l’estuaire de la Gironde (France). Ph.D. thesis, Université Paris-Est.
[36] Merkel U.H. and Kopmann R., 2012. Continuous vertical grain sorting for Telemac & Sisyphe v6p2,
Proceedings of the XIX Telemac-Mascaret User Conference. Oxford, UK.
[37] Meyer-Peter E. and Muller R., 1948. Formulae for bed-load transport. Int. IARH Congress, Stockholm, Sweden.
[38] Pham Van Bang D., Ovarlez G., Tocquer L., 2007. Density and structural effects on the rheological
characteristics of mud. La Houille Blanche, 2, 85-93.
Accessibility : Interne : EDF SA
Page 71 sur 73
c EDF SA 2013
EDF R&D
Reference manual
Sisyphe 6.3
H-P74-2012-02004-EN
Version 1.0
[39] Ribberink J., 1987. Mathematical modelling of one-dimensional morphological changes in rivers with
non-uniform sediment. PhD thesis, Delft University.
[40] Roelvink J.A., Coastal morphodynamic evolution techniques, Coastal Engineering, 2006, 53, 277-287.
[41] Sanchez, M., 1992, Modélisation dans un estuaire à mare : Rôle du bouchon vaseux dans la tenue des
sols sous marins. PhD thesis. University of Nantes (232 pages).
[42] Soulsby R., 1997. Dynamics of marine sands. Thomas Thelford Edition.
[43] Swart D.H., 1976. Offshore sediment transport and equilibrium beach profiles. Delft Hydraulics Publication 131, Delft University, The Netherlands.
[44] Talmon, A.M., Striksma, N. and van Mierlo, M.C.L.M. (1995). Laboratory measurements of the
direction of sediment tranport on transverse alluvial-bed slopes. Journal of Hydraulics Research,
33, 495-517.
[45] Thiébot J., 2008, Modélisation numérique des processus gouvernant la formation et la dégradation des
massif vaseux. PhD thesis ENGREF, University of Caen, 130 pp.
[46] Thiebot, J., Guillou, S., Brun-Cottan, J-C., 2011, An optimisation method for determining permeability and effective stress relationships of consolidating cohesive sediment deposits, Continental Shelf
Research 31 (2011) S117-S123.
[47] Toorman E. A., 1996, Sedimentation and self-weight consolidation : general unifying theory,
Géotechnique 46, N◦ 1, pp. 101-113.
[48] Toorman, E.A., 1999. Sedimentation and self-weight consolidation : constitutive equations and numerical modelling. Géotechnique, 49(6) :709-726.
[49] Tanaka, H. and Dang, V.T., 1996. Geometry of sand ripples due to combined wave-current flows.
J. of Waterway, Port, Coastal and Ocean Engineering, ASCE, 122 (6), 298-300.
[50] Terzaghi, K. 1923. Die Berechnung des Durchlassigkeitsziffer des Tones aus des Verlauf des hydrodynamischen Spannungerscheinungen, Sitz. Akad. Wissen. Wien, Math. Naturwiss. Kl., Abt IIa, 132,
125-138.
[51] van Rijn L.C., 1984. Sediment transport - Part I : bed load - Part II : suspended load, J. of
Hydraulic Division, Proc. ASCE, 110, HY10, 1431-56, HY11, 1613-41.
[52] van Rijn L.C., 1987. Mathematical modelling of morphological processes in the case of suspended sediment transport, Doctoral Thesis, Faculty of civil engineering, Delft University of technology).
[53] van Rijn L.C., 1993. Principles of Sediment Transport in Rivers, Estuaries and Coastal Seas. Aqua
Publications.
[54] van Rijn, L.C., 2007 Unified view of sediment transport by currents and waves. 1. Initiation of
motion, bed roughness, and bed-load transport. Journal of Hydraulic Engineering, 133, 6, 649-667.
[55] Villaret C., 2001. Prise en compte des fonds non-érodables dans le logiciel Sisyphe, Rapport EDF-LNH,
HP-75/02/045/B.
[56] Villaret C., 2001. Modélisation du transport solide par la formule de Bijker-Etude de sensibilité et tests
validation. Rapport EDF R&D LNHE HP-75/2001/66/A.
[57] Villaret C. and Davies A.G., 2004. Numerical modeling of littoral sand transport. ICCE Conference,
Portugal, September.
[58] Villaret C., Walther R., 2008. Numerical modeling of the Gironde estuary. Physics of Estuaries and
Coastal Sediments, Liverpool, August 2008.
Accessibility : Interne : EDF SA
Page 72 sur 73
c EDF SA 2013
EDF R&D
Reference manual
Sisyphe 6.3
H-P74-2012-02004-EN
Version 1.0
[59] Villaret C., Van L.A., Huybrechts N., Pham Van Bang D., Boucher O. 2010. Consolidation effects
on morphodynamics modelling : application to the Gironde estuary. La Houille Blanche, N◦ 6-2010,
15-24.
[60] Wiberg P.L. and Harris C.K., 1994. Ripple geometry in wave-dominated environments. Journal
of Geophysical Research 99 (C1), 775-789.
[61] Winterwerp J.C. and van Kesteren W.G.M., 2004. Introduction to the Physics of Cohesive Sediment
Dynamics in the Marine Environment. Elsevier, 576 pages.
[62] Waeles B., 2005. Detachment and transport of slay sand gravel mixtures by channel flows. Ph.D. thesis,
University of Caen.
[63] Yalin, MS, 1977. Mechanics of sediment transport. Pergamon Press, Oxford.
[64] Yalin M.S. and Ferrera da Silva A.M.F., 2001. Fluvial Processes. IAHR monograph Delft Hydraulics.
[65] Wu W., 2008. Computational River Dynamics. Taylor and Francis Group, London, UK.
[66] Zyserman J.A. and Fredsoe J. 1994. Data analysis of bed concentration of suspended sediment,
Journal of Hydraulic Engineering, ASCE, Vol. 120, N◦ 9, pp 1021-1042.
Accessibility : Interne : EDF SA
Page 73 sur 73
c EDF SA 2013
Source Exif Data:
File Type : PDF File Type Extension : pdf MIME Type : application/pdf PDF Version : 1.6 Linearized : No Create Date : 2014:01:13 10:07:58+01:00 Creator : pdfsam-console (Ver. 2.1.1e) Modify Date : 2018:01:15 17:04:30Z Has XFA : No XMP Toolkit : Adobe XMP Core 5.6-c015 84.159810, 2016/09/10-02:41:30 Creator Tool : pdfsam-console (Ver. 2.1.1e) Metadata Date : 2018:01:15 17:04:30Z Producer : iText 2.1.7 by 1T3XT Document ID : uuid:6e3f167e-5233-1846-b3d9-e2ef067c8b71 Instance ID : uuid:1be247a9-0b69-d14b-b354-e003f0a04625 Format : application/pdf Page Count : 76EXIF Metadata provided by EXIF.tools