SISYPHE User Manual links from [[User Manuals|User Manuals]], [[Manuals of SISYPHE|SISYPHE]] \\ ===== Introduction ===== \\ Sisyphe is a sediment transport and morphodynamic simulation module which is part of the hydroinformatic finite elements system Telemac. In this module, sediment transport rates, decomposed into bed-load and suspended load, are calculated at each grid point as a function of various flow (velocity, water depth, wave height, etc.) and sediment (grain diameter, relative density, settling velocity, etc.) parameters. The bedload is calculated by using one of the classical sediment transport formulae from the literature. The suspended load is determined by solving an additional transport equation for the depth-averaged suspended sediment concentration. The bottom evolution equation (the Exner equation) can be solved by using either a finite element or a finite volume formulation. Sisyphe is applicable to non-cohesive sediments (uniform or graded), cohesive sediments as well as to sand-mud mixtures. The sediment composition is represented by a finite number of classes, each characterized by its mean diameter, grain density and settling velocity. Sediment transport processes can also include the effect of bottom slope, rigid beds, secondary currents, slope failure, etc. For cohesive sediments, the effect of bed consolidation can be accounted for. Sisyphe can be applied to a large variety of hydrodynamic conditions including rivers, estuaries and coastal applications, where the effects of waves superimposed to a tidal current can be included. The bed shear stress, decomposed into skin friction and form drag, can be calculated either by imposing a friction coefficient (Strickler, Nikuradse, Manning or Chézy) or by a bed-roughness predictor. ==== Coupling with hydrodynamics ==== In Sisyphe , the relevant hydrodynamic variables can be either imposed in the model or calculated by a hydrodynamic computation (“chaining method” or “internal coupling”). It is convenient to use one of the hydrodynamic modules of the Telemac system (Telemac-2d , Telemac-3d or Tomawac) for compatibility reasons (same mesh, same pre- and postprocessor, etc.) and because modules are already coupled, but the user can also choose a different hydrodynamic model. The different methods which can be used to prescribe the hydrodynamics are described next. ==== Numerical kernel ==== Sisyphe can be run on the following OS: Unix, Linux and Windows. The latest release of Sisyphe (6.2) uses the version 6.2 of the Bief finite element library (Telemac system library). \\ ===== Running a sedimentological computation ===== ==== Input Files ==== * Minimum set of input files: CLI, SLF,CAS... * Additional/optional input files: F, REF, RES,... ==== Steering File ==== The steering file contains the necessary information for running a computation, plus the values of parameters that are different from the defaults. Like in other modules of the Telemac system, the model parameters can be specified in the (obligatory) Sisyphe steering file. The following essential information (input/output) should be specified in the steering file: * Input and output files * Physical parameters (sand diameter, cohesive or not, settling velocity, etc.) * Main sediment transport processes (transport mechanism and formulae, etc.) * Additional sediment transport processes (secondary currents, slope effect, etc.) * Numerical options and parameters (numerical scheme, options for solvers, etc.) ==== Input data files ==== The geometry file is a binary file that contains the finite element mesh information (e.g., element type, number of elements, coordinates of the nodes and connectivity matrix). This file can also contain bottom topography information and/or friction values at each mesh point. The geomety file (GEOMETRY FILE) is a binary file in Selafin format issued from the mesh generator. They contain a full description of the (three-node) triangular grid as well as the initial sediment bed level. This initial configuration can also be obtained from the interpolation of bathymetric data file onto each node of the triangular grid. The FORTRAN FILE contains those subroutines which the user has modified in order to adapt them to their specific needs. The HYDRODYNAMIC FILE contains the hydrodynamic results of a previous hydrodynamic computation (Telemac-2d or Telemac-3d ). If Telemac is coupled to Sisyphe , the grid must be the same in both cases. When a results file from a finite differences hydrodynamic model needs to be read, the hydrodynamic results file shall be transformed beforehand and interpolated onto the triangular mesh. The hydrodynamics will be given either by the last time step of that file in steady mode, or by the time steps related to the period of the tide or flood being considered in unsteady mode. For an explanation on steady/unsteady mode, see §XX). That file should contain current data and wave data if the effect of waves is considered. The WAVE FILE contains the hydrodynamic results of a previous wave computation (Tomawac ). Only the last record will be read if the keyword WAVE EFFECTS is activated. The PREVIOUS SEDIMENTOLOGICAL COMPUTATION FILE is provided for implementing a computational sequence, it contains the results of a previous sedimentological computation made using Sisyphe on the same grid. If that file exists, then the initial sedimentological conditions of the computation will be given by the last time step in the file. This file will be read only if the keyword COMPUTATION CONTINUED is activated. The REFERENCE FILE is the results file from a previous computation. The last time step is compared to the new last time step. All the files must be either binary (Selafin format) or text (ascii format). === Geometry file === A binary file containing information of the triangular grid. It can be used as input file for all hydrodynamic and wave propagation models in the Telemac Modelling System (Telemac-2d , Telemac-3d or Tomawac ). === Boundary conditions file === The format of the boundary condition file is the same as for Telemac-2d or Telemac-3d. This file can be created or modified using a text editor. Each line in this file is related to a grid boundary. Boundary points are listed in the file in the following way: First the domain outline points, proceeding counterclockwise from the lower left corner, then the islands proceeding clockwise , also from the lower left corner. The boundary conditions for sediment transport computations are imposed at the specific place of the tracer boundary condition as in Telemac-2d . __Obligatory keywords__\\ Names of principal input/output files (other files are optional) * BOUNDARY CONDITIONS FILE * GEOMETRY FILE * FORTRAN FILE * RESULTS FILE When combining a previous computation: * COMPUTATION CONTINUED (=NO :default option) * PREVIOUS SEDIMENTOLOGICAL COMPUTATION FILE \\ ===== Coupling hydrodynamics/morphodynamics ===== We describe here two possibilities for linking the hydrodynamic and the morphodynamic models: chaining (the flow is obtained from a previous hydrodynamic simulation assuming a fixed bed) or internal coupling (both the flow and bed evolution are updated at each time step). === Chaining method === * Principle Both models (hydrodynamic and morphodynamic) are run independently: during the first hydrodynamic simulation the bed is assumed to be fixed. Then, in the subsequent morphodynamic run, the flow rate and free surface are read from the previous hydrodynamic results file. This ‘chaining method’ is only justified for relatively simple flows, due to the difference in time-scales between the hydrodynamics and the bed evolution. For unsteady tidal flow, Sisyphe can be used in an unsteady mode: the flow field is linearly interpolated between two time steps of the hydrodynamic file. For steady flow, the last time step of the hydrodynamic file is used and the flow rate and free surface assumed to stay constant while the bed evolves. * Flow updating At each time step, the flow velocity is updated by assuming simply that both the flow rate and the free surface elevation are conserved, such that, in the case of deposition, the flow velocity is locally increased, whereas in the case of erosion, the flow velocity decreases. This rather schematic updating does not take into account any deviation of the flow. It is only suitable for simple flows (2D processes) and assuming relatively small bed evolutions. However it can be responsible for numerical instabilities (cf. [01b], Hervouet and Villaret, 2004). The morphodynamic computation is stopped when the bed evolution reaches a certain percent of the initial water depth. This simple updating of the flow field is no longer valid when the bed evolution becomes greater than a significant percentage of the water depth, specified by the user. At this point, it is recommended to stop the morphodynamic calculation and to recalculate the hydrodynamic variables. * Mass continuity It should be noted that with this simple method, the sediment mass continuity may not be satisfied because of potential losses due to changes in the flow depth as the bed evolves. When the flow is steady (‘steady case’ = YES), only the last record of the previous result file will be used. Otherwise (‘steady case’ = NO), the ‘Tide Period’ and ‘Number of tides or floods’ will be used to specify the sequence to be read on the hydrodynamic files. Hydrodynamic records are interpolated at each time step of the sedimentological computation. Note: an error may occur when the ‘Tidal Period’ is not a multiple of the graphical time steps of the hydrodynamic file (‘Hydrodynamic file is not long enough’). In an unsteady case, the keyword ‘Starting time of the hydrogram’ gives the first time step to be read. If the starting time is not specified, the last period of the hydrogram will be used for sedimentological computation. * Steering/fortran files For uncoupled mode, the Sisyphe steering file should specify: : - The time steps, graphical or listing output, duration : - The hydrodynamic file as yielded by Telemac -2D, -3D (‘hydrodynamic file’) or even by user-subroutine (condim_sisyphe.f). : - For waves only: the wave parameters can be either calculated by a wave propagation code (TOMAWAC), or defined directly in Sisyphe (condim_sisyphe.f). The effect of waves on bed forms and associated bed roughness coefficient can be accounted with keyword: Effect of waves =’YES’. : - A restart from a previous Sisyphe model run, by setting ‘Computation Continued = Yes’ and specification of sedimentological results in ‘Previous sedimentological computation ’ : - Flow options: steady or unsteady options, flow period __Keywords__\\ For time step, duration and output: * TIME STEP, NUMBER OF TIME STEPS * GRAPHIC PRINTOUT PERIOD * LISTING VARIABLES FOR GRAPHIC PRINTOUTS For hydrodynamics (imposed flow and updated): * HYDRODYNAMIC FILE * STEADY CASE (=NO, default option) * TIDE PERIOD (= 44640, default option) * STARTING TIME OF THE HYDROGRAM ( =0., default option) * NUMBER OF TIDES OR FLOODS (= 1, default option) * CRITICAL EVOLUTION RATIO (=0.1, default value) For waves: * WAVE FILE, WAVE EFFECTS === Internal coupling === * Principle For a more realistic simulation, Sisyphe can be automatically coupled (internally) with the hydrodynamic model, Telemac-2d or – 3D. Sisyphe is called inside the hydrodynamic model without any exchange of data files. The data to be exchanged between the two programs is now directly shared in the memory, instead of being written and read in a file. At each time step, the hydrodynamics variables (velocity field, water depth, bed shear stress,...) are transferred to the morphodynamic model, which sends back the updated bed elevation to the hydrodynamics model. * Time step and coupling period The internal coupling method is more CPU time consuming than the chaining method. Various techniques can be set up to reduce the CPU time (e.g. parallel processors, see §9). In certain cases, the use of a coupling period >1 allows the bed load transport rates and resulting bed evolution not to be re-calculated at every time step. For the suspended load, a diffusion/transport equation needs to be solved. This transport equation obeys the same Courant number criteria on the time step than the hydrodynamics, and therefore needs to be solved at each time-step (‘coupling period = 1’). The time step of Sisyphe is equal to the time step of Telemac-2d or – 3D multiplied by the ‘coupling period’. The graphic/listing printout periods are the same as in the Telemac computation. The Telemac-2d /3D steering file must specify the type of coupling, the name of the Sisyphe steering file, and the coupling period. In addition, the Fortran file of Sisyphe must be sometimes specified in the Telemac steering file (if there is no Fortran file for Telemac ). Some of the keywords of the Sisyphe steering file become obsolete. __Keywords__\\ For internal coupling, the following keywords need to be specified in the Telemac-2d (or –3D) steering files: * COUPLING WITH = SISYPHE * COUPLING PERIOD (=1, default value) * NAME OF SISYPHE STEERING FILE {{:panneau_attention_40.png |}}All computational parameters (time step, duration, printout, option for friction) need to be specified in the Telemac steering file, but are no longer used by Sisyphe . The values of time step, bottom shear stress... are transferred directly from Telemac to Sisyphe . keywords no longer in use is (Sisyphe ) steering file: * TIME STEP * GRAPHIC PRINTOUT PERIOD * LISTING PRINTOUT PERIOD * LAW OF BOTTOM FRICTION * COEFFICIENT \\ ===== Flow-sediment interactions ===== \\ ==== Sediment and fluid properties ==== Fine sediment particles of diameters less than $60\mu m$ present complex cohesive properties which affect the sediment transport processes. For non-cohesive sediments (median diameter $d_{50} >60 \mum$), the grain diameter and grain density $\rho_s$ are the key parameters which determine its resistance to erosion and sediment transport rate. For cohesive sediments ($d_{50} < 60 \mum$), the grain diameter is no longer the key sediment parameter: The settling velocity now depends on the concentration and other physico-chemical properties of the suspension due to flocculation, whereas the critical bed shear strength depends on the consolidation state of the sediment bed. For the simplest case, we consider uniform, non-cohesive sediments characterized by one single value for the grain size $d_{50}$ and grain density $\rho_s=$ with $\rho_s=$2650 kg\,m$^{-3}$ which can be transported both as bed-load and suspended load. __Keywords__\\ The physical properties of the sediment are always defined in the Sisyphe steering file using the following keywords: * COHESIVE SEDIMENTS (= NO , default option) * NUMBER OF SIZE CLASSES OF BED MATERIAL (NSICLA= 1, default option) * SEDIMENT DIAMETERS ($d_{50} > 60\mu$m for non-cohesive sediments) * SEDIMENT DENSITY ($\rho_s = 2650$ kgm$^{-3}$, default value) * SETTLING VELOCITIES The physical properties of the fluid are defined by: * VISCOSITY OF THE FLUID ($\nu=$ 10$^{-6}$ m$^2$s$^{-1}$, by default) * DENSITY OF THE FLUID ($\rho_{0}=$1000 kgm$^{-3}$, by default) * GRAVITY ACCELERATION ($g=9.81$ ms$^{-2}$) ==== Bed shear stress ==== === Hydrodynamic model === The current-generated bed shear stress is used in both the shallow water momentum equation as well as the bottom boundary condition for the velocity profile. When Sisyphe is coupled with Telemac-2d , the bed shear stress t0 is calculated at each time step from the depth-averaged mean flow velocity by use of a quadratic friction coefficient Cd as follows: \begin{equation*} \tau _0 =\frac{1}{2} \rho C_d U^{2} \end{equation*} where $\rho$ is the fluid density, and $U$ the depth-averaged velocity. When \sisyphe is coupled with \telddd, the bed shear stress is aligned with the near bed velocity in order to account for possible vertical flow deviations. The magnitude of the bed shear stress is still related to the depth-averaged velocity, except if a Nikuradse friction law is applied. In that case, the friction velocity u$_{*}$, defined by \begin{equation*} \tau _0 = \rho u_*^2 \end{equation*} is related to the near bed flow velocity $u(z_1)$) by a logarithmic velocity profile: \begin{equation*} u(z_1) = \frac{u_*}{\kappa} \ln\left(\frac{z_1}{z_0}\right) \end{equation*} where $\kappa = 0.4$ is the von K\'arm\'an constant, $z_0$ is the vertical distance from a rough boundary, expressed as a function of the Nikuradse bed roughness ($z_0 = k_s/30$), with $k_s$ the grain roughness height. For flat beds, the roughness height has been shown to be approximately $k_s\approx 3 d_{50}$ === Role of bed forms === A natural sediment bed is generally covered with bed forms (length $\lambda$ and height $\eta$). The presence of bed forms greatly modifies the boundary layer flow structure, with the formation of recirculation cells and depressions in the lee of bed forms. Depending on the flow rate and the sediment transport rate, the size of bed forms ranges from a few centimeters for ripples to a few tens of meter for mega-ripples. The dimensions of dunes scale with the water depth ($\eta_D \approx 0.4 h$, $\lambda_D\approx 6-10h$). In most cases, large scale models do not resolve the small to medium scale bed forms (ripples, mega-ripples) which need therefore to be parameterized by an increased friction coefficient: \begin{equation*} \tau_0 = \tau' + \tau'' \end{equation*} The local skin friction component determines the bed-load transport rate and the equilibrium concentration for the suspension. The total friction velocity determines the (spatially averaged along bed forms) turbulence eddy viscosity/diffusivity vertical distribution in 3D models, and therefore determines both velocity vertical profile and mean concentration profile. ==== Sediment transport ==== === Bedload and suspended load === When the current-induced bed shear stress increases above a critical threshold value, sediment particles start to move as bed load, while the finer particles are transported in suspension. Bedload occurs in a very thin high concentrated near-bed layer, where inter-particle interactions occur. The suspended load is defined as the depth-integrated flux of sediment concentration, from the top of the bedload layer up to the free surface. For dilute suspension ,clear flow concepts (turbulence diffusion, eddy viscosity, logarithmic velocity profile) are considered to be valid.The total sediment load $Q_t$ includes both a bedload $Q_b$ and suspended load $Q_s$: \begin{equation*} Q_t = Q_b + Q_s. \end{equation*} === Shields parameter === Most sand transport formulae assume a critical value of the bed shear stress, noted $\tau_c$, for the onset of erosion: $Q_s = 0$ for $\tau_0 < \tau_c$. The Shields parameter or non-dimensional critical shear stress $\theta_c$ is defined by: \begin{equation*} \theta_c = \frac{\tau_c}{g(\rho_s -\rho)d_{50}} \end{equation*} can be either specified in the \sisyphe steering file by use of keyword SHIELDS PARAMETERS, or calculated by the model as a function of non-dimensional grain diameter. The Van Rijn formula has been programmed \begin{equation*} \begin{align*} D_* \leq 4,\quad & \theta_c = 0.24 D_*^{-1} \\ 4 < D_* \leq 10,\quad & \theta_c = 0.14 D_*^{-0.64} \\ 10 < D_* \leq 20,\quad & \theta_c = 0.04 D_*^{-0.10} \\ 20 < D_* \leq 150,\quad & \theta_c = 0.013 D_*^{0.29} \\ 150 \leq D_*,\quad & \theta_c = 0.045, \end{align*} \end{equation*} where the non-dimensional diameter $D_*$ is defined by: \begin{equation*} D_* =\left(\frac{g(\rho_s/\rho-1)}{\nu^2}\right) ^{1/3} d_{50}. \end{equation*} __Keywords__ * BED-LOAD (= YES : default option) * SUSPENDED LOAD (=NO : default option) * SHIELDS PARAMETER’ (default option: Van Rijn formula) ==== Friction laws ==== === Differences between 2D and 3D models === The intensity of bed shear stress is calculated by use of same friction laws in both 2D and 3D models: the bed shear stress is related to the depthaveraged flow velocity, except when the Nikuradse friction law is applied. In 3D-models, the bed shear stress is directly related to near bed velocity assuming a log profile (eq. 2), whereas in 2D-models, the log profile is assumed to be valid up to the free surface. The direction of the bed shear stress and resulting bedload transport rate is assumed to be in the direction of the depth-averaged velocity in Sisyphe (alone or internally coupled with Telemac-2d ). When Sisyphe is internally coupled to Telemac-3d , the bed shear stress and resulting transport rate are assumed to be in the direction of the near bed velocity. The 3D model gives a more accurate estimate of the bottom friction, since it accounts for a possible vertical deviation of the current. === Coupling with TELEMAC-2D === The depth-averaged bed shear stress and resulting bedload transport rates are assumed to be in the direction of the mean flow velocity, except when the sediment transport formulation accounts for: * deviation correction due to sloping bed effects * secondary currents due to river meandering {{:panneau_attention_40.png |}}**Coupling with TELEMAC-2D**\\ When the model is coupled with Telemac-2d , the values of the friction coefficients (and therefore, the bed shear stress) are provided by Telemac-2d . === Coupling with TELEMAC-3D === Both bed shear stress and bed load transport rate are aligned with the near bed velocity. This 3D approach is more physical, and takes into account possible recirculation and veering of the flow, such that corrections for secondary currents for example are no longer necessary. === Uncoupled model === The quadratic friction coefficient Cd which is used to calculate the total bed shear stress can be calculated based on the selected friction law. Different options, which are consistent with the Telemac-2d options, are available in Sisyphe and depend on the choice of the ‘LAW FOR BOTTOM FRICTION’ and on the value of the ‘FRICTION COEFFICIENT’. Keywords: * Chézy coefficient Ch (KFROT = 2) : \begin{equation*} C_d = \frac{2g}{C_h^2} \end{equation*} * Strickler coefficient $S_t$ (KFROT = 3) : \begin{equation*} C_d = \frac{2g}{S_t^2}\frac{1}{h^{1/3}} \end{equation*} * Manning friction $M_a$ (KFROT = 4) : \begin{equation*} C_d =\frac{2g}{h^{1/3} } M_a^2 \end{equation*} * Nikuradse bed roughness $k_s$ (KFROT = 5) : \begin{equation*} C_d = 2\left[\frac{\kappa}{\log(\frac{12h}{k_s})}\right]^2 \end{equation*} __Keywords__\\ In the Sisyphe steering file, the total bed shear stress is calculated based on: * LAW OF BOTTOM FRICTION (KFROT =3, default option) * BOTTOM FRICTION COEFFICIENT (St=50, default value) Similar keywords are available in both Telemac-2d and -3D in the case of internal coupling. ==== Skin friction correction ==== === Total bed shear stress decomposition === Bedload transport rates are calculated as a function of the local skin friction component t0. The total bed shear stress issued from the hydrodynamics model needs to be corrected in the morphodynamics model as follows: \begin{equation*} \tau'= \mu \tau_0, \end{equation*} Physically, the skin bed roughness should be smaller than the total bed roughness (i.e. $\mu \leq 1$). However, in most cases the hydrodynamic friction does not represent the physical bottom friction: the coefficient is generally used as a calibration coefficient in hydrodynamics models. It is adjusted by comparing simulation results with observations of the time-varying free surface and velocity field. Therefore, its model value integrates various neglected processes (side wall friction, possible errors in the bathymetry and input data). Under those conditions, a correction factor $\mu > 1$ can be admitted. Different methods are programmed in \sisyphe in order to calculate the bedform correction factor $\mu$, according to the keyword `Skin friction correction'(ICR) : * ICR = 0: no correction: the total friction issued from Telemac is directly used for sand transport calculations ($\mu$=1). * ICR = 1: the skin roughness is assumed to be proportional to the sand grain diameter like in the case of flat beds ($k_s' \sim d_{50}$). The proportionality coefficient is specified by the keyword ‘Ratio between skin friction and mean diameter’ (KSPRATIO). \begin{equation*} \mu =\frac{C_d'}{C_d}, \end{equation*} where $C_d$ et $C_d'$ are both quadratic friction coefficients related respectively to total friction and skin friction: Cd is obtained from Telemac-2d or Telemac-3d and $C_d'$ is calculated from $k_s'$: \begin{equation*} C_d' = 2\left[\frac{\kappa}{\log(\frac{12h}{k_s'})}\right]^2 \end{equation*} * ICR = 2: we use a bedform predictor to calculate the bedform roughness $k_r$. Both $k_r$ and $k_s'$ should influence the transport rates. Here we assume: \begin{equation*} \mu =\frac{C_d'^{0.75} C_r^{0.25}}{C_d}, \end{equation*} where the quadratic friction $C_r$ due to bedforms is calculated as a function of $k_r$. __Keywords__ The option selected for the skin friction correction is based on keywords: * SKIN FRICTION CORRECTION (ICR=0, default option) * RATIO BETWEEN SKIN FRICTION AND MEAN DIAMETER (KSPRATIO=3,default value) ==== Bed roughness predictor ==== === Total bed roughness (uncoupled mode) === Different options are programmed in Sisyphe to predict the total bed roughness. It is recalled that the bed friction option of Sisyphe is not used in the case of internal coupling with Telemac . * For KFROT = 0 : the bed is assumed to be flat $k_s = k_s'= \alpha d_{50}$ * For KFROT = 1: the bed roughness is decomposed into different components (skin, bedforms, etc.) and each component is calculated by Sisyphe depending on the flow regime. For waves and combined waves and currents, bedform dimensions are calculated as a function of wave parameters following the method of Wiberg and Harris. The wave-induced bedform bed roughness $k_r$ is calculated as a function of the wave-induced bedform height $\eta_r$: \begin{equation*} k_r = \max(k_s', \eta_r). \end{equation*} For currents only, the van Rijn predictor has been implemented. The total bed roughness can be decomposed into a grain roughness $k_s'$, a small-scale ripple roughness $k_r$, a mega-ripple component $k_{mr}$, and a dune roughness $k_d$: \begin{equation*} k_s = k_s' + \sqrt{k_r^2 + k_{mr}^2 + k_d^2}. \end{equation*} Both small scale ripples and grain roughness have an influence on the sediment transport laws, while the mega-ripples and dune roughness only contribute to the hydrodynamic model (total friction). === Bedform correction factor === In order to account for the effect of ripples on the skin friction, the keyword ‘Skin friction correction’ should be set to ICR = 2. For currents only, the ripple bed roughness is a function of the mobility number, according to Van Rijn : \begin{equation*} \begin{align*} k_r &= d_{50}(85-65\tanh(0.015(\Psi-150))),\,\,\text{for\,\,} \Psi<250\\ k_r &= 20 d_{50},\,\,\text{for\,\,} \Psi \geq 250, \end{align*} \end{equation*} where $\Psi =U^2/(s-1)gd_{50}$ __Keywords__ * LAW OF BOTTOM FRICTION (KFROT =0 or 1, flat bed or rippled bed for total friction) * SKIN FRICTION CORRECTION (ICR=1 or 2 flat bed or rippled bed , for skin friction correction) \\ ===== Bed-load transport ===== \\ ==== Exner equation ==== === Equilibrium conditions === Bed-load occurs in a thin near-bed high-concentrated region. The bedload layer therefore adapts very rapidly to any changes in the flow conditions, such that equilibrium conditions can be considered to be valid. The bed-load transport rate can then be calculated by use of some equilibrium sediment transport formula, as a function of various flow and sediment parameters, assuming that the transport rate corresponds to saturation conditions. === Bed evolution === To calculate the bed evolution, Sisyphe solves the Exner equation: \begin{equation*} (1-n)\frac{\partial Z_f}{\partial t} + \grad\cdot Q_b =0, \end{equation*} where $n$ is the bed porosity ($n \approx 0.4$ for non cohesive sediment), $Z_f$ the bottom elevation, and $Q_b$ (m$^{2}$s$^{-1}$) the solid volume transport (bed-load) per unit width. The previous Equation states that the variation of sediment bed thickness can be derived from a simple mass balance. It is only valid for equilibrium conditions. The previous Equation is strictly valid for bed-load only. However, it can be extended to total load (including the suspended load) assuming quasi-steady and uniform flow conditions. Otherwise, a different treatment should be applied to both bed load and suspended load, in order to take non-equilibrium flow conditions into account properly. __Keywords__ * BEDLOAD = yes (default option) * SUSPENDED LOAD = No (default option) * COEFFICIENT FUNCTION OF THE POROSITY (XKV=1/(1- n)=1.6, default value) === Bedload transport formulae === For currents only (no wave effects), a large number of semi-empirical formulae can be found in the literature to calculate the bedload transport rate. Sisyphe offers the choice between 10 different bedload formulae including the Meyer-Peter and Müller, Engelund-Hansen and Einstein- Brown formulae. Most sediment transport formulae assume threshold conditions for the onset of erosion (e.g. Meyer-Peter andMüller, van Rijn and Hunziker). Other formulae are based on similar energy concept (e.g. Engelund- Hansen) or can be derived from a statistical approach (e.g. Einstein- Brown, Bijker, etc.). The non-dimensional current-induced sand transport rate, noted $\Phi_s$, can be expressed as a function of the non-dimensional skin friction or Shields parameter, noted $\theta'$ , as defined by: \begin{equation*} \theta'=\frac{\mu\tau_0}{(\rho_s-\rho)gd_{50}}, \end{equation*} $\mu$ , the correction factor for skin friction \begin{equation*} \Phi_s = \frac{Q_s}{\sqrt{g(s-1)d_{50}^3}} \end{equation*} __Keywords__ * BEDLOAD TRANSPORT FORMULA (default option is MP formula, ICF=1) ==== Sediment transport formulae ==== === Choice of formulae === The choice of a transport formula depends on the selected value of the model parameter ICF , as defined in the steering file by the keyword ‘bedload transport formula’ (ICF). The user can program a specific transport formula through the Qsform.f(user-subroutine), by setting ICF=0. Other sediment transport formulae, described in chapters 8 and 9, have been programmed in Sisyphe to account for the effects of waves (cf. Bijker, Bailard, Dibajnia and Watanabe, etc.) or sand grading (cf. Hunziker). They can nevertheless be used for currents only and uniform grain size. * ICF= 1: Meyer-Peter formula This classical bed-load formula has been validated for coarse sediments in the range ($0.4$ mm $< d_{50} < 29$ mm). It is based on the concept of initial entrainment: \begin{equation*} \begin{align*} \Phi_b &=0\quad\text{if}\,\theta'<\theta_c \\ \Phi_b &=8(\theta'-\theta_c)^{3/2}\quad\text{otherwise} \end{align*} \end{equation*} where $\theta_c$ is the critical Shields parameter ($\theta_c = 0.047$ by default). * ICF = 2: Einstein-Brown formula. This bed-load formula is recommended for gravel ($d_{50} > 2$ mm) and large bed shear stress $\theta > \theta_c$. The solid transport rate is expressed as: \begin{equation*} \Phi_b = F(D_*)f(\theta'), \end{equation*} with \begin{equation*} F(D_*) = \left(\frac{2}{3} +\frac{36}{D_*}\right)^{0.5} - \left( \frac{36}{D_*}\right)^{0.5}, \end{equation*} and \begin{equation*} \begin{align*} f(\theta') = 2.15\exp(-0.391/\theta') & \quad\text{if}\,\,\theta' \leq 0.2 \\ f(\theta') = 40\,\theta^{'3} & \quad\text{if}\,\,\theta' >0.2 \end{align*} \end{equation*} * ICF= 3 or 30: Engelund-Hansen formula. The Engelund-Hansen formula predicts the total load (bed load plus suspended load). It is recommended for fine sediments, in the range ($0.2$ mm $< d_{50} < 1$ mm but beware that the use of a total load formula is only suitable under equilibrium conditions (quasi steady and uniform flow). Two different forms of the same equation are programmed in Sisyphe: ICF = 30 corresponds to the original formula, where the transport rate is related to the skin friction without threshold: \begin{equation*} \Phi_s = 0.1\,\theta'^{5/2} \end{equation*} ICF = 3 corresponds to the version modified by Chollet and Cunge in order to account for the effects of sand dunes formation. The transport rate is here related to the total bed shear stress: \begin{equation*} \Phi_s = \frac{0.1}{C_d} \hat{\theta}^{5/2}, \end{equation*} the dimensionless bed shear stress $\theta$ is calculated as a function of the dimensionless skin friction $\theta'$: \begin{equation*} \begin{align*} \hat{\theta} = 0 &\quad \text{if}\,\theta' < 0.06\,\text{(flat bed regime - no transport)} \\ \hat{\theta} = \sqrt{2.5(\theta'-0.06)} &\quad \text{if}\, 0.06 < \theta' < 0.384\,\text{(dune regime)} \\ \hat{\theta} = 1.065\theta'^{0.176} &\quad \text{if}\, 0.384 < \theta' < 1.08\,\text{(transition regime)} \\ \hat{\theta} = \theta' &\quad \text{if}\,\theta' > 1.08\,\text{(flat bed regime - upper regime)} \end{align*} \end{equation*} {{:panneau_attention_40.png |}}**Notice**\\ In the case of coupling between bed-load and suspended load (SUSPENSION = YES), a total load formula (N$^\circ$3 and 30) should not be used, in order to avoid the suspended load to be calculated twice. * ICF = 7: Van Rijn formula. This formula was proposed by Van Rijn (1984a) in order to calculate the bed-load transport rate for particles of size $0.2$ mm $< d_{50} < 2$ mm: \begin{equation*} \Phi_b = 0.053 D_*^{-0.3} \left( \frac{\theta_p-\theta_{cr}}{\theta_{cr}} \right)^{2.1}. \end{equation*} {{:panneau_attention_40.png |}}**Validity range**\\ Most sediment transport formulae are based on experiments performed under steady-flow conditions. They demonstrate rapid variation of the bed-load transport rates as a function of the mean flow intensity: an increase of the current velocity by $10\%$ will result, depending on the formula being used, in an increase of the transport rate by over $30\%$ (Meyer-Peter), $60\%$ (Engelund-Hansen) or almost $80\%$ (Einstein-Brown). Therefore, any error made when calculating the hydrodynamics will be significantly amplified by the sediment transport rates estimates. On the other hand, under variable flow conditions (e.g. in tidal regime), the average transport will be highly influenced by the stronger currents and will not be directly related to the mean flow. The validity range of the different formulae can be summarized in the table below: ^ Formula ^ Meyer-Peter ^ Einstein-Brown ^ Engelund-Hansen ^ van Rijn | ^ IFC | 1 | 2 | 3 or 30 | 7 | ^ Mode of transport | Bed-load | Bed-load | Total-load | Bed-load | ^ Sand diameter | > 1 mm | 0.2-3 mm | 0.2-4 mm | 0.2-2 mm | ^ Rippled bottoms | Yes | Yes | Yes | XXX | ^ Low bed-load flows | Yes | Yes | No | XXX | ^ High bed-load flows | Yes | Yes | Yes | XXX | ==== Sloping bed effect ==== The effect of a sloping bottom is to increase the bedload transport rate in the downslope direction, and to reduce it in the upslope direction. In Sisyphe , we apply a correction factor to both the magnitude and direction of the solid transport rate, before solving the bed-evolution equation. This sloping bed effect is only activated if the keyword ‘SLOPE EFFECT’ is activated. Two different formulations for both effects are available depending on the choice of ‘FORMULA FOR SLOPE EFFECT’ (SLOPEFF), which chooses the magnitude correction and ‘FORMULA FOR DEVIATION’ (DEVIA), which chooses the direction correction. === Magnitude of bedload transport rate === * For SLOPEFF=1. The first correction method is based on the Koch and Flokstra’s formula (1981). The solid transport rate intensity Qb0 is multiplied by a factor: \begin{equation*} 1-\beta\frac{\partial Z_f}{\partial s}, \end{equation*} where $s$ is the coordinate in the current direction and $\beta $ is an empirical factor, which can be specified by keyword 'BETA' (default value = 1.3). This effect of bed slope is then similar to adding a diffusion term in the bed-evolution equation. It tends to smoothe the results and reduces unstabilities. * For SLOPEFF=2, the method of Soulsby (1997) is programmed: The threshold bed shear stress is modified as a function of bed slope and depends on the ‘FRICTION ANGLE OF THE SEDIMENT' $\phi_s$ \begin{equation*} \frac{\theta_c}{\theta_{co}} = \frac{\cos\psi \sin\beta + ( \cos^2\beta \tan^2\phi_s - \sin^2\psi \sin^2\beta)^{0.5}}{\tan \phi_s}, \end{equation*} where $\beta$ is the angle of the bed to the horizontal plane and $\psi$ the angle of the current to the upslope direction. === Direction of bedload transport rate === The change in the direction of solid transport is taken into account by the formula: \begin{equation*} \tan\alpha = \tan\delta -T\frac{\partial Z_f}{\partial n}, \end{equation*} where $\alpha$ is the direction of solid transport in relation to the flow direction, $\delta$ is the direction of bottom stress in relation to the flow direction, and $n$ is the coordinate along the axis perpendicular to the flow. : - DEVIA = 1: According to Koch and Flochstra, the coefficient T now depends on the Shields parameter: : \begin{equation*} T=\frac{4}{6\theta} \end{equation*} : - DEVIA = 2, the deviation is calculated based on Talmon et al. (1995) and depends on the empirical coefficient $\beta_2 = 0.85$. : \begin{equation*} T=\frac{1}{\beta_2\sqrt{\theta}} \end{equation*} This default value can be modified by changing the keyword: `PARAMETER FOR DEVIATION' = $\beta_2$. __Keywords__ * SLOPE EFFECT (=YES, default option) * FORMULA FOR SLOPE EFFECT (SLOPEFF=1, default option) * FORMULA FOR DEVIATION (DEVIA=1, default option) * FRICTION ANGLE OF THE SEDIMENT ($\phi_s = 40^\circ$, default value) * BETA (($\beta = 1.3$, default value) * PARAMETER FOR DEVIATION (($\beta_2 = 0.85$, default value) ==== Sediment slide ==== An iterative algorithm prevents the bed slope from becoming greater than the maximum friction angle ($\phi_s\approx 32^{\circ} - 40^{\circ}$). This option is activated by use of key word `SEDIMENT SLIDE'. Further developments are needed to introduce the bed stability, in order to represent mud slides. \begin{equation*} \left\| \overrightarrow{\grad(Z_f)} \right\| <\tan \phi_s \end{equation*} __Keywords__ * SEDIMENT SLIDE (= NO, default option) * FRICTION ANGLE OF SEDIMENT ($\phi_s=40$, default value) ==== Secondary currents ==== In meandering streams, centrifugal inertial forces are induced by the curvature $R^{-1}$ of the streamlines. The depth-averaged mean pressure gradient due to the cross-sectional variation of the free surface balances the depth -averaged centrifugal force, noted F: \begin{equation*} \begin{align*} F_y &= -\frac{1}{\rho}\frac{\partial P}{\partial y} \\ F_y &= \rho \frac{U^2}{R}-\frac{1}{\rho} \frac{\partial P}{\partial y} =-g\frac{\partial Z_s}{\partial y}, \end{align*} \end{equation*} where $y$ is the cross-sectional axis and $\partial Z_s/\partial y$ the cross-sectional variation in the free surface. The mass balance in the y-direction gives a zero-depth averaged cross-sectional velocity. \begin{equation*} \rho \frac{U^2}{R} =-g\frac{\partial Z_s}{\partial y} \end{equation*} However, in the vertical, the centrifugal force increases with distance from the bed whereas the pressure gradient is constant. P and F are therefore out of balance (P greater in the lower part of the flow, whereas F is greater in the upper part). This results in a cross-sectional circulation, with a cross-velocity directed near the inner bank, and flow velocity therefore follows a spiral motion which cannot be represented by a depth-averaged hydrodynamic model. The correction method, based on Engelund (? ), has been programmed in Sisyphe , in order to account for the effect of these radial stresses on the sand transport rates. This method is based on the following assumptions: * the bottom shear stress is constant in a cross-section * the bed roughness is also constant * the mean water depth h is constant The near bed angular deviation of the main flow $\delta$ can be shown to be proportional to $h/R$, as explained in Yalin and Ferreira da Silva (2001). We use here $\tan\delta = 7h/R$. The bend radius $R$ is generally unknown, but can be inferred from the cross sectional variation in the free surface: \begin{equation*} R=-\rho\alpha'\frac{U^2}{g\frac{\partial Z_s}{\partial y}}, \end{equation*} \begin{equation*} R=-\rho\alpha'\frac{U^2}{g\frac{\partial Z_s}{\partial y}}, \end{equation*} __Keywords__ * SECONDARY CURRENTS (=NO, default option) {{:panneau_attention_40.png |}}**Notice**\\ This option should only be activated when the flow is calculated by a depth-averaged model. ==== Treatment of rigid beds ==== === Different methods === Non-erodible beds are managed numerically by limiting bedload erosion and letting incoming sediment pass over. The problem of rigid beds is conceptually trivial but numerically complex. In finite elements the minimum water depth algorithm allows a natural treatment of rigid beds (cf. Hervouet et al. 2011). The sediment is managed as a layer with a depth that must remain positive, and the Exner equation is managed similarly to the shallow-water continuity equation. With finite volumes, two methods are implemented in Sisyphe , depending on the ‘OPTION FOR THE TREATMENT OF NON ERODABLE BEDS’ : * Option (0, default option), the bed is movable everywhere. The new algorithm for minimum water depth is used * Options (4) are mass-conservative and avoid erosion of the rigid beds. __Keywords__ * OPTION FOR THE TREATMENT OF NON ERODABLE BEDS (=0, default option) * MASS-LUMPING (=YES, default option) __Subroutine noerod.f__ The position of the rigid bed is set by default to a very large negative value and can be changed in subroutine noerod.f. ==== Effect of tidal flats ==== Tidal flats are areas of simulation domain where depth can become 0 during simulation can be managed by two approaches: * In the default option (OPTBAN = 1), equations are solved everywhere. The water depth is set to a minimum value HMIN (key-word MINIMAL VALUE OF THE WATER HEIGHT, default value 1.D-3) in order to avoid divisions by zero. However, in case of coupling with Telemac-2d or Telemac-3d , the original depth is kept for suspension, because a strict compatibility with Saint-Venant continuity equation is required which would be spoiled by clipping. This is the recommended option. * Another option (OPTBAN = 2), removes from the computation elements with points that have a depth less than a minimum HMIN. This option must be considered obsolete, it is actually not exactly mass-conservative and does not work properly in parallel mode. __Keyword__ * TIDAL FLATS (=YES , default option) * MINIMAL VALUE OF THE WATER HEIGHT (HMIN = 10−3m ,default value) * OPTION FOR THE TREATMENT OF TIDAL FLATS (OPTBAN = 1, default option) \\ ===== Suspended load ===== \\ ==== Suspended load transport equation ==== === Passive scalar hypothesis === For non-equilibrium flow conditions, the bedload and the suspended load are dealt with separately treatment. The interface between the bed-load and suspended load is located at $z=Z_{ref}$: * In the thin high-concentrated bed-load layer ($z Z_{ref}$), clear flow concepts still apply, and the sediment grains can be regarded as a passive scalar which follows the mean and turbulent flow velocity, with an additional settling velocity term. === Settling velocity === The settling velocity Ws is an important parameter for the suspension. It can be either specified, or calculated by the model as a function of grain diameter, by use of a semi-empirical formula. The Van Rijn formula which is valid for non-cohesive spherical particles and dilute suspensions, has been programmed in Sisyphe (cf. Van Rijn, 1993, cf. [31], p. 3.13): \begin{equation*} \begin{align*} W_{s} &= \frac{(s-1)g d_{50}^2}{18\nu}, \quad \text{if } d_{50} \leq 10^{-4} \\ W_{s} &= \frac{10\nu}{d_{50}} \left(\sqrt{1+0.01\frac{(s-1)gd_{50}^3}{18\nu^2}}-1\right), \quad \text{if } 10^{-4} \leq d_{50} \leq 10^{-3}\\ W_{s} &= 1.1 \sqrt{(s-1)gd_{50}}, \quad \text{if } 10^{-3} \leq d_{50}, \end{align*} \end{equation*} where $s=\rho_{s}/\rho_0$ is the relative density and $g$ is gravity. Outside the bed-load layer sediment particles are assumed to follow the mean flow velocity, noted $\overrightarrow U$, with an additional vertical settling velocity, noted $W_s$. === 3D transport equation === The 3D advection-diffusion equation for the suspended sediment concentration is: \begin{equation*} \frac{\partial C}{\partial t} + \grad\cdot \left(( \overrightarrow{U} +\overrightarrow{W_s}) C\right) = \grad\cdot (\epsilon_s \grad C), \end{equation*} where $\epsilon_s$ is the sediment diffusivity coefficient. The following boundary conditions must be satisfied both at the free surface ($Z_s$) and at the interface between the bed-load and the suspended load ($z=Z_{ref}$): \begin{equation*} \begin{align*} \text{at } z &= Z_s,\quad \epsilon_s \frac{\partial C}{\partial z} + W_s C = 0 \\ \text{at } z &= Z_{ref},\quad -\epsilon_s \frac{\partial C}{\partial z} - W_sC = (E-D)_{Z_{ref}} \end{align*} \end{equation*} where $Z_s$ is the free surface elevation $Z_{ref}$ the position of the interface between the bed-load layer and the suspension. $E-D$ is the net erosion minus deposition flux of sediment. The 3D transport equation and boundary conditions can be solved in the 3D sediment transport model of Telemac-3d. === 2D transport equation === The 2D transport equation for the depth-averaged mean concentration C is obtained by depth-integration : \begin{equation*} \frac{\partial hC}{\partial t} + \frac{\partial h\overline{uC}}{\partial x} + \frac{\partial h\overline{vC}}{\partial y} = \frac{\partial }{\partial x} \left(h\epsilon_s\frac{\partial C}{\partial x}\right) + \frac{\partial }{\partial y} \left(h\epsilon_s\frac{\partial C}{\partial y}\right) + E - D, \end{equation*} with: \begin{equation*} X = \overline{x} =\frac{1}{h} \int_{Z_{ref}}^{Z_s} x(z) \d z, \end{equation*} where $h=Z_s-Z_f \approx Z_s-Z_{zref}$ is the water depth, assuming the bed-load layer thickness to be small. After simplification of the advection terms and using the continuity equation, the following approximate depth-averaged transport equation can be solved in its non-conservative form: \begin{equation*} \frac{\partial C}{\partial t} + U_{conv} \frac{\partial C}{\partial x} + V_{conv} \frac{\partial C}{\partial y} = \frac{1}{h} \left(\frac{\partial }{\partial x} \left( h\epsilon_s \frac{\partial C}{\partial x} \right) + \frac{\partial }{\partial y} \left(h\epsilon_s \frac{\partial C}{\partial y} \right) \right) + \frac{(E-D)_{z=Zref}}{h}, \end{equation*} where $U_{conv}$ and $V_{conv}$ are the depth-averaged convective flow velocities in the x and y directions. * Convection: The convective velocity can be corrected from the depthaveraged mean velocity in order to account for the fact that most sediment is transported near the bed. * Diffusion The diffusion terms can be set to zero (‘DIFFUSION’ = NO) According to the choice of the parameter ‘OPDTRA’ (key word), the diffusion terms can be simplified, and equation becomes \begin{equation*}\ \frac{\partial C}{\partial t} + U_{conv} \frac{\partial C}{\partial x} + V_{conv} \frac{\partial C}{\partial y} = \left(\frac{\partial }{\partial x} \left( \epsilon_s \frac{\partial C}{\partial x} \right) + \frac{\partial }{\partial y} \left(\epsilon_s \frac{\partial C}{\partial y} \right) \right) + \frac{(E-D)_{z=Zref}}{h}, \end{equation*} The value of the dispersion coefficient $\epsilon_s$ depends on the choice of the parameter OPTDIF and different possibilities are offered (constant coefficient, by default). The Elder model (OPTDIF=2) allows to calculate both the longitudinal and transversal components as a function of the friction velocity. ($\epsilon_s = 6 u_*h$ and $\epsilon_s = 0.3 u_*h$). Different schemes are available for solving the non-linear advection terms depending on the choice of the parameter the classical characteristics schemes (to the diffusive schemes SUPG and PSI, it is recommended to use the conservative schemes __Keywords__ * SUSPENSION (= NO by default) * DIFFUSION (= YES, by default) * OPTION FOR THE DIFFUSION OF TRACER (OPDTRA=1,the default option corresponds to the simplified expression of the diffusion term) * OPTION FOR THE DISPERSION (=1, default option) * DISPERSION ALONG THE FLOW ($\epsilon_s = 10^{-2}$ m$^{2}$s$^{-1}$, default value) * DISPERSION ACROSS THE FLOW ($\epsilon_s = 10^{-2}$ m$^{2}$s$^{-1}$, default value) * TYPE OF ADVECTION (RESOL=1 , method of characteristics by default) {{:panneau_attention_40.png |}}**Tidal flats**\\ It is recommended in the presence of tidal flats to use the conservative scheme based on the calculation of flux per segments. Scheme 14 should be used if the convection velocity differs from the depth-averaged velocity and no longer satisfies the shallow-water continuity equation (‘correction on convection velocity’ =YES ). ==== Bed evolution ==== === Bed evolution due to the suspension === The bottom-evolution equation due to the suspension accounts for the vertical fluxes at the interface between the bed-load and the suspended load. The bed-evolution equation becomes: \begin{equation*} (1-n)\frac{\partial Z_f}{\partial t} + (E-D) = 0 \end{equation*} where $Z_f$ is the bottom elevation, $Z_{ref}$ the interface between the bed-load and suspended-load, and $n$ the bed porosity. Erosion and deposition fluxes are calculated at the bed-load/suspended load interface ($z = Z_{ref}$). {{:panneau_attention_40.png |}}**Units consistency**\\ The net sediment flux E-D is in m s$^{-1}$. For consistency with the bed-load units (m$^{2}$s$^{-1}$), concentrations and bed porosity are dimensionless. However, the user can choose mass concentration for graphic printouts, by use of key word `MASS CONCENTRATION'. The relation between volume concentration and mass concentration is : $C_m$ (kgm$^{-3}$) = $\rho_s C$. __Keywords__ * MASS CONCENTRATION (NO, default option) * BEDLOAD =NO (YES, default option) * SUSPENDED LOAD =YES (NO, default option) * ‘coefficient function of the porosity’ ==== Equilibrium concentrations ==== === Erosion/Deposition rate === For non-cohesive sediments, the net erosion minus deposition flux is determined based on the concept of equilibrium concentration (Celik and Rodi, 1988): \begin{equation*} \left(E-D \right)_{Z_{ref}} = W_s \left(C_{eq} - C_{Z_{ref}}\right) \end{equation*} where $C_{eq}$ is the equilibrium near-bed concentration and $C_{zref}$ is the near-bed concentration, calculated at the interface between the bed-load and the suspended load, $z=Z_{ref}$. Two different semi-empirical formula from the literature are programmed in Sisyphe : * ICQ= 1: Zyserman and Fredsoe formulation: \begin{equation*} C_{eq} =\frac{0.331(\theta'-\theta_c)^{1.75}}{1+0.72(\theta'-\theta_c)^{1.75}}, \end{equation*} where $\theta_c$ is the critical Shields parameter and $\theta'= \mu\theta$ ,the non-dimensional skin friction which is related to the Shields parameter, by use of the ripple correction factor. * 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: \begin{equation*} C_{eq} = \frac{Q_sB}{b Z_{ref} u_*} \end{equation*} with an empirical factor b = 6.34. {{:panneau_attention_40.png |}}**For ICQ= 2...**\\ For ICQ= 2, both bed-load and suspended load should be calculated at each time step. The ‘PERIOD OF COUPLING’ should be set to 1(coupled mode) and ‘BED-LOAD’ = YES. === Reference level === The reference elevation $Z_{ref}$ corresponds to the interface between bed-load and suspended-load. For flat beds, the bed-load layer thickness is proportional to the grain size, whereas when the bed is rippled, the bed-load layer thickness scales with the equilibrium bed roughness ($k_r$). The definition of the reference elevation needs also to be consistent with the choice of equilibrium bed-concentration formula. * ICQ= 1: According to Zyserman and Fredsoe, the reference elevation should be set to $Z_{ref} = 2d_{50}$. In \sisyphe we take $Z_{ref}=k'_s$, where $k'_s$ is the skin bed roughness for flat bed ($k_s' \approx d_{50}$), the default value of proportionality factor is KSPRATIO =3). * ICQ= 2: According to Bijker, $Z_{ref} = k_r$. === Vertical concentration profile === We assume here a Rouse profile for the vertical concentration distribution, which is theoretically valid in uniform steady flow conditions: \begin{equation*} C(z)=C_{Z_{ref}} \times \left(\frac{z-h}{z}\frac{a}{a-h}\right)^R, \end{equation*} where $R$ is the Rouse number defined by: \begin{equation*} R=\frac{W_s}{\kappa u_*}, \end{equation*} and $\kappa$ the von Karman constant ($\kappa = 0.4$), and $u_*$ the friction velocity corresponding to the total bed shear stress as commented before, and $a$ is the reference elevation above the bed elevation. Ratio between the reference and depth-averaged concentration By depth-integration of the Rouse profile, the following relation can be established between the mean (depth-averaged) concentration and the reference concentration: \begin{equation*} C_{Zref} = F C, \end{equation*} where: \begin{equation*} F^{-1} = \left(\frac{Z_{ref}}{h}\right)^R\int_{Zref/h}^1\left(\frac{1-u}{u}\right)^R \d u. \end{equation*} In Sisyphe, we use the following approximate expression for $F$: \begin{equation*} \begin{align*} F^{-1} &= \frac{1}{\left(1-Z\right)} B^R\left( 1-B^{(1-R)} \right), \quad \text{for } R \neq 1\\ F^{-1} &= -B \log B, \quad \text{for } R = 1, \end{align*} \end{equation*} with $B = Z_{ref}/h$. __Keywords__ * REFERENCE CONCENTRATION FORMULA (ICQ) ==== Convection velocity ==== A straightforward treatment of the advection terms would imply the definition of an advection velocity and replacement of the depth-averaged velocity U by: \begin{equation*} U_{conv} = \overline{UC}/C. \end{equation*} A correction factor is introduced in Sisyphe , defined by: \begin{equation*} F_{conv} =\frac{U_{conv}}{U} \end{equation*} The convection velocity should be smaller than the mean flow velocity ($F_{conv} \leq 1$) since sediment concentrations are mostly transported in the lower part of the water column where velocities are smaller. We further assume an exponential profile concentration profile which is a reasonable approximation of the Rouse profile, and a logarithmic velocity profile, in order to establish the following analytical expression for $F_{conv}$: \begin{equation*} F_{conv} =-\frac{I_2 - \ln\left(\frac{B}{30}\right) I_1}{I_1 \ln\left( \frac{eB}{30}\right)}, \end{equation*} with $B=k_s/h = Z_{ref}/h$ and \begin{equation*} I_1 =\int_B^1\left(\frac{(1-u)}{u}\right)^R \d u,\quad I_2 = \int_B^1 \ln u \left(\frac{(1-u)}{u} \right)^R \d u, \end{equation*} see Huybrechts et al., 2010. __Keywords__ * CORRECTION ON CONVECTION VELOCITY’ (= NO, default option) ==== Initial concentrations, boundary conditions ==== === Initial conditions === The initial concentration for the suspended load can be either imposed within condim_susp.f or specified in the steering file through the keyword ‘Initial suspension concentrations ’ initializes the value of the volume concentration for each class. === Boundary conditions === For the boundary conditions, the concentration of each class can be specified in the steering file through keyword: ‘concentration per class at boundaries’. It may be also convenient to use keyword ‘equilibrium inflow concentration =Yes’: the concentration at the entrance of the domain and at t = 0 is set to its equilibrium value, according to the choice of the ‘reference concentration formula’. Input concentrations can be also directly specified (user subroutines: conlit.f). === Equilibrium conditions === The concentrations at the entrance of the domain can be calculated by Sisyphe assuming equilibrium conditions in order to avoid unwanted bed-evolution at the entrance of the domain, and also at the first time step, it is possible to impose the concentration to its equilibrium value, by activating the keyword ‘EQUILIBRIUM INFLOW CONCENTRATION’. The equilibrium (depth-averaged) concentration is then calculated assuming equilibrium concentration at the bed and a Rouse profile correction for the F factor. __Keywords__ * INITIAL SUSPENSION CONCENTRATIONS (CS0 = 0, default value) * EQUILIBRIUM INFLOW CONCENTRATION (=NO, default option) * CONCENTRATION PER CLASS AT BOUNDARIES (=0, default value) === User subroutines === The subroutine condim_susp.f can be used to specify the initial conditions for the sediment concentration .The subroutine conlit.f can be used to specify the concentration at the entrance of the domain. \\ ===== Wave effects ===== \\ ==== Theoretical background === In the coastal zone, the effect of waves, superimposed to a mean (waveinduced or tidal) current modifies the bottom boundary layer structure. Due to the reduced thickness of the bed boundary layer, bottom shear stresses are largely enhanced and resulting sand transport rates are an order of magnitude greater than in the case of currents alone. In addition, the wave-generated ripples also have a strong effect on the bed roughness and sand transport mechanisms. These effects will only be accounted for when the keyword ‘effect of waves’ is activated. ==== Linear wave theory ==== The wave orbital velocity is calculated by Sisyphe , assuming linear theory: \begin{equation*} U_0=\frac{H_s \omega }{2 \sinh (kh)}, \end{equation*} where $\omega = 2\pi/T_p$ is the wave pulsation, $k = 2\pi/L$ , where $L$ is the wave length. The wave number is calculated from the wave dispersion relation: \begin{equation*} \omega^2 = gk\tanh (kh). \end{equation*} === Wave propagation === The wave parameters (wave height, wave period, wave direction) need to be specified. They can be simply imposed in the Fortran file (user subroutine: condim_sisyphe.f ) or calculated by a previous TOMAWAC simulation [10]. The ‘HYDRODYNAMIC FILE’ can contain all the wave parameters (which is recommended for the non-steady case). The waves can also be specified from an independent wave field. The ‘wave file’ gives the name of a previous wave (TOMAWAC) results file. Only the last record will be read. === Wave parameters === The wave-parameters, required as input parameters for the morphodynamic computation, are the significant wave height $H_s$ and the peak wave period $T_p$. The wave-angle $\theta_w$ needs also to be specified in some particular applications, e.g. when using for example the Bailard's formula. (If not, the default value is $90^\circ$, which means that the direction of propagation of the wave is parallel to the $x$-axis). If the wave field is calculated by TOMAWAC, the results hydrodynamic file should contain the following variables: * HM0: significant wave height (corresponds to $H_s$) * TPR5: peak period (corresponds to $T_p$) * DMOY: mean wave direction, relative to $y$-axis (corresponds to $\theta_w$) The (last time record of the) ‘wave file’ issued from previous TOMAWAC computation is read by Sisyphe . In addition, the mean current issued from previous Telemac computation is read from the ‘hydrodynamic file’. For unsteady flow conditions (time-varying velocity and wave field), the user should run Tomawac in an unsteady mode, and store all the wave and hydrodynamic relevant variables on the same input ‘hydrodynamic file’. __Keywords__ * WAVE EFFECTS (=NO, by default) * WAVE FILE User-subroutine condim_sisyphe.f ==== Wave-induced bottom friction ==== === Wave-induced bottom friction === The maximum stress due to waves is calculated at each time step as a function of the wave-orbital velocity $U_w$ by use of a quadratic friction coefficient $f_w$ due to waves: \begin{equation*} \tau_w = \frac{\rho}{2} f_w U_w^2. \end{equation*} The wave friction factor $f_w$ is calculated as a function of relative density: \begin{equation*} f_w = f_w \left( A_0/k_s \right), \end{equation*} where $A_0$ is the orbital amplitude on the bed ($A_0 = U_w/\omega$), and $k_s$ the bed roughness. In Sisyphe, we use the Swart (1974) formula: \begin{equation*} \begin{align*} f_w &= \exp \left(-6.0 + 5.2\left( \frac{A_0}{k_s} \right)^{-0.19}\right), \quad \text{if } \frac{A_0}{k_s} > 1.57\\ f_w &= 0.30, \quad \text{if } \frac{A_0}{k_s} \leq 1.57. \end{align*} \end{equation*} === Wave-current interaction === For combined waves and currents, the wave-induced (mean and maximum) bottom stresses are an order of magnitude larger than in the case of currents alone. Different models can be found in the literature to calculate the wave and current bottom stress, noted $\tau_{cw}$, as a function of the bottom shear stress due to currents only, $\tau_c$ , and the maximum shear stress due to waves only, $\tau_w$. In Bijker : \begin{equation*} \tau_{cw} = \tau_c + \frac{1}{2} \tau_w \end{equation*} The mean ($\tau_{mean}$) and maximum ($\tau_{\max}$) shear stress can also be calculated following Soulsby's method. Non-dimensional variables are defined: \begin{equation*} X=\frac{\tau_0}{\tau_0+\tau_w};\quad Y=\frac{\tau_{moy}}{\tau_0+\tau_w};\quad Z=\frac{\tau_{\max}}{\tau_0+\tau_w} \end{equation*} They can be parameterized based on different turbulence models: \begin{equation*} \begin{align*} Y &= X\left(1+bX^p(1-X)^q\right),\\ Z &= 1+aX^m(1-X)^n \end{align*} \end{equation*} The various model coefficients ($a, b, m, n, p, q$) are empirical functions of friction coefficients $f_w$, $C_D$ and also of the wave-current angle $\phi$: \begin{equation*} \begin{align*} a &= (a_1 +a_2 \left|\cos \phi \right|^I)+(a_3+a_4 \left|\cos \phi \right|^I)\log_{10}\left(\frac{2f_w}{C_D} \right),\\ b &= (b_1 +b_2 \left|\cos \phi \right|^J)+(b_3+b_4 \left|\cos \phi \right|^J)\log_{10} \left(\frac{2f_w}{C_D } \right), \end{align*} \end{equation*} and similar expressions for $m, n, p$, and $q$. The 26 fitted constants ($a_i$, $b_i$, $m_i$, $n_i$, $p_i$, $q_i$, $I$ and $J$) depend on the turbulence model selected. The coefficients which have been used in \sisyphe correspond to the model of Huynh-Thanh et Temperville (1991) [24b], $I = 0.82$, $J = 2.7$: ^ $i$ ^ $a_i$ ^ $b_i$ ^ $m_i$ ^ $n_i$ ^ $p_i$ ^ $q_i$ | ^1 | -0.07 | 0.27 | 0.72 | 0.78 | -0.75 | 0.89 | ^2 | 1.87 | 0.51 | -0.33 | -0.23 | 0.13 | 0.40 | ^3 | -0.34 | -0.10 | 0.08 | 0.12 | 0.12 | 0.50 | ^4 | -0.12 | -0.24 | 0.34 | -0.12 | 0.02 | -0.28 | === Friction coefficient under combined waves and current === The characteristic shear-stress to represent wave current interactions, noted $\left\langle \left| \overrightarrow{\tau_{cw}} (t)\right| \right\rangle$, is related to the time averaged mean-velocity: \begin{equation*} \left\langle \left| \overrightarrow{\tau_{cw}} (t)\right| \right\rangle =\rho f_{cw} \left\langle \left| \overrightarrow{U(t)} \right|^2 \right\rangle, \end{equation*} where \begin{equation*} \left\langle \left| \overrightarrow{U(t)} \right|^2 \right\rangle =U_c^2 +\frac{1}{2} U_w^2 \end{equation*} According to Camenen (2002, page 60), $\left\langle \left| \overrightarrow{\tau_{cw}} (t)\right| \right\rangle$ is taken as a weighted average between $\tau_{moy}$ and $\tau_{\max}$: \begin{equation*} \left\langle \left| \overrightarrow{\tau_{cw}} (t)\right| \right\rangle = X\tau_{moy} + (1-X)\tau_{\max}, \end{equation*} which is equivalent to: \begin{equation*} \left\langle \left| \overrightarrow{\tau_{cw}} (t)\right| \right\rangle = Y\tau_c + Z\tau_w. \end{equation*} The final expression for the wave-current interaction factor is: \begin{equation*} f_{cw} = \frac{Y\tau_c + Z\tau_w}{{U_c^2} + \frac{1}{2}{U_w^2}}. \end{equation*} This expression of the wave and current friction factor is used in Bailard and in Dibajnia and Watanabe sand transport formulae. ==== Wave-induced ripples ==== Equilibrium ripples are generally observed outside the surf zone. Their dimensions can be predicted as a function of waves (orbital velocity $U_0$ and wave period $T=2\pi/\omega$), for a given uniform sediment diameter $d_{50}$. We follow in \sisyphe the procedure of Wiberg and Harris ([34b], 1994). This formulation is only applicable for oscillatory flow conditions and does not account for the effect of a superimposed mean current. Ripples can be classified into three types depending on the value of the ratio of wave orbital diameter to mean grain diameter, $D_0/d_{50}$, where $D_0 = 2U_0/\omega$. The ripples dimensions, namely their wave length $\lambda$ and height $\eta$, can be calculated as follows: In the orbital ripples regime, under moderate wave conditions, ripples dimensions are simply proportional to $D_0$: \begin{equation*} \begin{align*} \lambda &= 0.62D_0 \\ \eta &= 0.17\lambda. \end{align*} \end{equation*} For larger waves, in the anorbital ripple regime, ripples dimensions scale with the sand grain diameter: \begin{equation*} \begin{align*} \lambda &= 535 d_{50} \\ \eta &= \lambda \exp \left(-0.095\left(\log\frac{D_0}{\eta} \right)^2 + 0.042\log\frac{D_0}{\eta}-2.28\right). \end{align*} \end{equation*} The effect of ripples is to increase the bed roughness. As in Bijker's formula, we assume in Sisyphe: \begin{equation*} k_s = \max \left(\eta ,3d_{50} \right). \end{equation*} The effect of the mean current superimposed to the waves is accounted for by introducing a correction factor, noted $\alpha$, to the orbital excursion, following Tanaka and Dang (1996): \begin{equation*} \alpha = 1 + 0.81\left(\tanh (0.3S_*^{2/3})\right)^{2.5} \left(\frac{U}{U_w}\right)^{1.9}, \end{equation*} with \begin{equation*} S_* = d_{50}\frac{\sqrt{(s-1)g d_{50}}}{4\nu}. \end{equation*} __Keywords__ * BOTTOM FRICTION LAW kfrot= 1 (3, default option) ==== Wave-induced sand transport ==== The main effect of waves is to increase the sand transport rates in comparison to the case of currents alone. This effect can be accounted for by using one of the wave sand transport formula which have been programmed in Sisyphe : the Bijker’s, the Soulsby – Van Rijn formula ,Bailard’s and the Dibajnia and Watanabe formulae. In Sisyphe , the solid flow is assumed to be oriented in the direction of the mean current. The possible deviation of the transport in the direction of waves, which can be important in the near shore zone due to non-linear effects, is not accounted for. * Bijker formula : The Bijker formula (1968) can be used for determining the total transport rate Qs (bed-load + suspension), with its two components, the bed-load Qsc and the suspended-load Qss being determined separately. The development is made on the scalar variables, since Bijker’s formula is an extension of a steady flow formula to account for the effect of the wave enhanced shear stress. : The suspended load transport is treated in a simplified manner: the concentration profile is assumed to be continuously in equilibrium. The inertia effects related to the water column loading and unloading will then in no way be modelled. In addition, no exchange takes place with the bed load layer, only the continuity of concentration is ensured in that case. : - Bed-load transport : Bijker extended the steady bed-load formula proposed by Frijlink (1952), adding the effect of the wave into it. In nondimensional variables, the bed-load transport rate is: : \begin{equation*} \Phi_b = b\theta_c^{0.5}\exp\left(-0.27\frac{1}{\mu \theta_{cw}}\right), \end{equation*} : where $\theta_c$ is the non-dimensional shear stress due to currents alone, $\theta_{cw}$ the non-dimensional shear stress due to wave-current interaction, and $\mu$ is a correction factor which accounts for the effect of ripples. : In the original formulation, we have b = 5. Recent studies showed that b = 2 provides better results in comparison to 43 field data and other engineering models. Then the formula is programmed in Sisyphe with b = 2. The ripple correction factor ì is calculated in the same way as in the Meyer Peter formula. : - Suspended load transport : The suspended load transport is obtained by vertically integrating, from a reference level corresponding to the thickness of the bed load layer, the product of the volume concentration by the mean horizontal velocity. After depth-integration and assuming a Rouse profile for the concentration and a logarithmic velocity profile for the mean velocity profile, the suspended load can be written as : \begin{equation*} Q_{ss} = Q_{sc} I, \end{equation*} : with : \begin{equation*} I=1.83\times 0.216\frac{B^{A-1}}{(1-B)^A} \int_B^1 \left(\frac{1-y}{y}\right)^A \ln\left(\frac{33y}{B}\right) \d y, \end{equation*} : and : \begin{equation*} \begin{align*} A &= \frac{W_s}{\kappa u_*},\quad u_*=\sqrt{\frac{\tau_{cw}}{\rho}}\\ B &= k_s/h. \end{align*} \end{equation*} : : __Keywords__ : * BED-LOAD TRANSPORT FORMULA: ICF= 4 : * WAVE EFFECTS’ = YES : : – The Soulsby - Van Rijn formula.The transport rate Qs (volume/volume) due to the combined action of waves and current is provided by the following equation: : : \begin{equation*} Q_{b,s} = A_{b,s} U\left[ \left( U^2+2\frac{0.018}{C_D} U_0^2\right)^{0.5}-U_{cr}\right]^{2.4}. \end{equation*} : : This formula can be applied to estimate both components of the total sand transport rate (bed-load + suspension), and it is suitable for rippled bed regimes. : : $A_{b,s}$ are the bed-load/suspended load coefficients, calculated by: : : \begin{equation*} \begin{align*} A_b &= \frac{0.005 h \left(d_{50}/h\right)^{1.2}}{\left((s-1)gd_{50}\right)^{1.2}} \\ A_s &= \frac{0.012 d_{50}D_*^{-0.6}}{\left((s-1)gd_{50}\right)^{1.2}}, \end{align*} \end{equation*} : : where $U$ is the depth-averaged current velocity, $U_0$ is the RMS orbital velocity of waves, and $C_D$ is the quadratic drag coefficient due to current alone. This formula has been validated assuming a rippled bed roughness: $k_s = 0.18$ m : : The critical entrainment velocity, $U_{cr}$ is given by: : : \begin{equation*} \begin{align*} U_{cr} &= 0.19 d_{50}^{0.1}\log_{10}\left(\frac{4h}{D_{90}}\right), \quad \text{if } 0.1 mm\leq d_{50}\leq 0.5 mm \\ U_{cr} &= 8.5 d_{50}^{0.6} \log_{10}\left(\frac{4h}{D_{90}} \right), \quad \text{if } 0.5 mm\leq d_{50}\leq 2.0 mm. \end{align*} \end{equation*} : : {{:panneau_attention_40.png |}}**Validity range** : $h = 1-20$ m, $U = 0.5-5$ms$^{-1}$, $d_{50}=0.1-2$mm. : : __Keywords__ : * BED-LOAD TRANSPORT FORMULA: ICF= 5 : * WAVE EFFECTS’ = YES : In the Van-Rijn-Soulsby formula, the diameter $D_{90}$, characteristic of the coarser grains, should also be specified by use of the keyword: `D90'. * Bailard’s formula Bailard’s formula [13] The Bailard’s formula is based on an energetic approach. The bedload and the suspended-load components of the sand transport rate are expressed respectively as the third- and fourth-order momentum of the near-bed time-varying velocity field $\overrightarrow{U(t)}$ \begin{equation*} \begin{align*} Q_b &= \frac{f_{cw}}{g(s-1)} \frac{\epsilon_c}{\tan \varphi} \left\langle \left| \overrightarrow{U} \right|^2 \overrightarrow{U} \right\rangle \\ Q_s &= \frac{f_{cw}}{g(s-1)} \frac{\epsilon_s}{W_s} \left\langle \left| \overrightarrow{U} \right|^3 \overrightarrow{U} \right\rangle \end{align*} \end{equation*} with $f_{cw}$ the friction coefficient which accounts for wave-current interactions, $\epsilon_s$, $\epsilon_c$ empirical factors ($\epsilon_s = 0.02$, $\epsilon_c = 0.1$), $\varphi$ sediment friction angle ($\tan \varphi$ = 0.63) and $<\cdot>$ time-averaged over a wave-period. Under combined wave and currents, the time-varying velocity vector $\overrightarrow{U(t)}$ is composed of a mean component $U_c$ and an oscillatory component of amplitude $U_0$, assuming linear waves; $\phi$ is the angle between the current and the wave direction, $\phi_c$ is the angle between the mean current direction and the x-axis, $\phi_w$ is the angle between the wave direction of propagation and the x-axis. The time-varying velocity field verifies: \begin{equation*} \overrightarrow{U(t)} =U_c \exp^{i\phi_c} + U_0 \cos \omega t \exp^{i\phi_w}. \end{equation*} For the third-order term, one can derive an analytical expression in the general case, whereas for the fourth-order momentum, we have to assume the waves and currents to be co-linear ($\phi = 0$), in order to find an analytical expression. (For the general case of a non-linear non-colinear wave, we would have to integrate numerically the fourth-order velocity term (cf. Camenen, 2002 [14a]). This method is however not efficient. * Third-order term: \begin{equation*} \begin{align*} \left\langle\left| \overrightarrow{U} \right|^2 \overrightarrow{U} \right\rangle &= \overrightarrow{U_c}(U_c^2+\frac{1}{2} U_0^2)+ \overrightarrow{U_0} (\overrightarrow{U_c}\cdot\overrightarrow{U_0})\\ \left\langle \left| \overrightarrow{U} \right|^2 \overrightarrow{U} \right\rangle_x &= \left[ U_c^3 + U_c U_0^2 (1+\frac{1}{2} \cos 2\phi )\right] \cos \phi_c -\frac{1}{2} U_c U_0^2\sin 2\phi \sin \phi_c\\ \left\langle \left| \overrightarrow{U} \right|^2 \overrightarrow{U} \right\rangle_y &= \left[ U_c^3 + U_c U_0^2 (1+\frac{1}{2} \cos 2\phi )\right] \sin \phi_c +\frac{1}{2} U_c U_0^2\sin 2\phi \cos \phi_c \end{align*} \end{equation*} * Fourth-order term: \begin{equation*} \begin{align*} \left\langle \left| \overrightarrow{U} \right|^3 \overrightarrow{U} \right\rangle_x = \frac{1}{8} (24 U_0^2 U_c^2 + 8 U_c^4 +3U_0^4)\cos \phi_c\\ \left\langle \left| \overrightarrow{U} \right|^3 \overrightarrow{U} \right\rangle_y = \frac{1}{8} (24 U_0^2 U_c^2 + 8 U_c^4 +3U_0^4)\sin \phi_c \end{align*} \end{equation*} __Keywords__ * BED-LOAD TRANSPORT FORMULA: ICF= 8 * WAVE EFFECTS’ = YES * Dibajnia and Watanabe formula The Dibajnia and Watanabe (1992) formula is an unsteady formula, which accounts for inertia effects. The sand transport rate is calculated by: \begin{equation*} \frac{\overrightarrow{Q_s} }{W_s d_{50}} =\alpha \frac{\overrightarrow{ \Gamma} }{\left| \Gamma \right|^{1-\beta}}, \end{equation*} with $\alpha=0.0001$ and $\beta=0.55$, empirical model coefficients. $\overrightarrow{\Gamma}$ represents here the difference between the amounts of sediments transported in the onshore and offshore directions. This formula is used to estimate the intensity of the solid transport rate, the direction is assumed to follow the mean current direction. Time-varying velocity field In the wave direction, the time-varying velocity is given by: \begin{equation*} U(t)=U_c \cos \phi + U_0 \sin \omega t, \end{equation*} $r$ is defined here as the asymmetry parameter: \begin{equation*} r=\frac{U_0}{U_c\cos \phi}. \end{equation*} The wave cycle is decomposed into two parts: 1. During the first part of the wave-cycle ($0 < t < T_1$), the velocities are directed onshore ($U(t) > 0$). 2. During the second part ($T_2 = T-T_1$), the velocities are negatives ($U(t) <0$). The sand transport rate in the wave direction is: \begin{equation*} \Gamma_{X'} = \frac{u_1T_1(\Omega_1^3+\Omega_2^{'3})-u_2T_2(\Omega_2^3+\Omega_1^{'3})}{(u_1+u_2)T} \end{equation*} where $\Omega_1$ and $\Omega_2$ represent the amount of sand transported respectively in the onshore and offshore directions which will be deposited before flow reversal; $\Omega'_1$ et $\Omega'_2$ represent the remaining part which stays in suspension after flow eversal. The inertia of sediment depends on the ratio $d_{50}/W_s$, and represented by parameter $\omega_i$ defined by: \begin{equation*} \omega_i = \frac{u_i^2}{2(s-1)gW_s T_i}, \end{equation*} with \begin{equation*} u_i^2 = \frac{2}{T_i} \int_{u>\text{or}<0} u^2 \d t \end{equation*} 1. all sediment is deposited before flow reversal $\omega_i\leq \omega_{crit}$ \begin{equation*} \Omega_i=\omega_i\frac{2W_s T_i}{d_{50}}, \quad \Omega'_i = 0 \end{equation*} 2. part of the sediment stays in suspension after flow reversal $\omega_i \geq \omega_{crit}$ \begin{equation*} \Omega_i = \omega_{crit} \frac{2W_s T_i}{d_{50}}, \quad \Omega'_i=(\omega_i-\omega_{crit})\frac{2W_sT_i}{d_{50}}. \end{equation*} The critical value of $\omega_{crit}$ is calculated as a function of the wave-current non-dimensional friction: \begin{equation*} \Theta_{cw} = \frac{<\tau_{cw} >}{\rho}\frac{1}{(s-1)gd_{50}} \end{equation*} According to Temperville et Guiza (2000): \begin{equation*} \begin{align*} \omega_{crit} &= 0.03, \quad\text{if }\Theta_{cw}\leq 0.2 \\ \omega_{crit} &= 1-\sqrt{1-((\Theta_{cw}-0.2)/0.58)^2}, \quad\text{if } 0.2\leq\Theta_{cw}\leq 0.4 \\ \omega_{crit} &= 0.8104\sqrt{\Theta_{cw}} - 0.4225, \quad\text{if } 0.4\leq\Theta_{cw}\leq 1.5\\ \omega_{crit} &= 0.7236\sqrt{\Theta_{cw}} - 0.3162, \quad\text{if } 1.5\leq\Theta_{cw}. \end{equation*} __Keywords__ * BED-LOAD TRANSPORT FORMULA: ICF= 9 * WAVE EFFECTS’ = YES \\ ===== Sand grading effects ===== \\ ==== Sediment bed composition ==== === Granulometry distribution === The number of classes of bed material is specified in the steering file and limited to $NSICLA < xxx$. For uniform sediment, the granulometry distribution is represented by one single value ($NSICLA =1$) for the whole domain. The mean grain diameter and corresponding settling velocity are defined in the steering file. For sediment mixtures, the granulometry distribution is discretized in a number of classes. Each class of sediment noted $j$ ($1 < j < NSICLA$) is defined by its mean diameter $d_{50}(j)$ and volume fraction in the mixture at every node, noted $AVAI(j)$. The characteristics of each class of sediment (mean diameter, Shields number $\theta_c (i)$, settling velocities $W_s(i)$ of each class can be specified in the steering file or calculated by the model (same Van Rijn formulae than for single class sediment). The percent of each class of material need to verify: \begin{equation*} \sum\limits_{j=1,NSICLA}AVAI(j)=1 \quad \text{and} \quad 0\leq AVAI(j) \leq 1. \end{equation*} The initial bed composition can be defined in the steering file for simple configuration (uniform bed). For a more complex bed configuration (spatial variation, vertical structure), the user subroutine {\texttt init\_compo.f} must be used to define the initial state. The mean diameter ($D_m$, noted $ACLADM$) is calculated based : \begin{equation*} D_m = \sum\limits_{j=1,NSICLA} AVAI(j) D(j). \end{equation*} === Bed structure === The bed is stratified into a number of layers ($NOMBLAY$) that enables vertical variations of the sediment bed composition. The percentage of each class $j$ of material (noted $AVAI(i,j,k)$ as well as the thickness of each layer $E_s$ are defined at each point $i$ and for each layer $k$. The composition of transported material is computed using the composition of the upper layer (see below the definition of an `active layer'). The initial composition of the bed (number of layers, thickness of each layer $E_s$, composition of each layer $AVAI$ are specified in user-subroutine init_compo.f. The subroutine init_avail.f checks that the structure of the initial bed composition compatible with the position of the rigid bed (as defined in user subroutine noerod.f): \begin{equation*} Z_f(i)-Z_r(i) = \sum_{k=1,NOMBLAY} E_s(k). \end{equation*} In general, the thickness of the bed is very large (by default, $100$ m), so the bottom layer thickness need to be increased. Assuming the porosity of each class to be identical (parameter $XKV =1/(n-1)$ is a constant), the total mass of sediments per unit area is: \begin{equation*} M_s(i) = \sum_{k=1,NOMBLAY}\rho_s E_s(k)(n-1). \end{equation*} In each layer, the AVAI (percent of each class) which are defined as the volume of each class of material per the total volume of material, are considered to be constant. In each layer (and at each point) the constraints on the AVAI need to be satisfied : \begin{equation*} \sum_{j=1,NSICLA} AVAI(i,k,j) = 1\quad \text{and}\quad 0 \leq AVAI(i,k,j) \leq 1. \end{equation*} __Keywords__\\ Initial sediment bed composition is defined by: * NUMBER OF SIZE CLASSES OF BED MATERIAL (nsicla = 1, by default) * SEDIMENT DIAMETERS (01;.01; . . . by default) * SETTLING VELOCITIES (if not specified calculated by SISYPHE) * SHIELDS PARAMETERS * INITIAL FRACTION FOR PARTICULAR SIZE CLASS (ava0 = 1; 0. ;0. . . . , by default) For more than one size-classes, the user should specify NSICLA values for the mean diameter, initial fraction. . . For example, for a mixture of two sands * NUMBER OF SIZE-CLASSES OF BED MATERIAL = 2 * INITIAL FRACTION FOR PARTICULAR SIZE CLASS =.5;.5 * MEAN DIAMETER OF THE SEDIMENT =.02;.01 * SHIELDS PARAMETERS =.045;.01 * SETTLING VELOCITIES =.10;.05 If not specified the settling velocities and Shield parameters are calculated as a function of grain size for each particles. ==== Transport of sediment mixtures ==== The method programmed in SISYPHE for the treatment of mixtures of sediment with variable grain sizes is classical and based on previous models from the literature (as for example the 1D model Sedicoup developed at Sogreah). There are two layer concepts implemented in Sisyphe. The active layer model based on Hirano’s concept ([XX], 1971) (default) and a newly developed continues vertical sorting model (C-VSM) from Merkel ([XX], 2012), ([XX], 2012). === Active layer/Active stratum === Since only a thin upper layer will be transported, the upper layer is therefore subdivided into a thin ‘active layer’ (k = 1) and an ‘active stratum’ (k = 2). The composition of the active layer is used to calculate the composition of transported material and the intensity of transport rates for each class of sediments: \begin{equation*} Q_s =\sum_{j=1,NSICLA} AVA0(j) Q_s(j), \end{equation*} where $AVA0(j)$ is the percentage of class $j$ in the active layer. The active layer thickness is not well defined and depend on the flow and sediment transport parameters (cf. Van Rijn, 1987). According to Yalin (1977), it is proportional to the sand diameter of the upper layer ($1.2d_{90}$). The active later thickness generally scales with the bed roughness, which is typically of the order of a few grain diameters for flat beds up to a few cm in the case of rippled beds. In the presence of dunes, the active layer thickness should be half the dune height (cf. Ribberink, 1987). In Sisyphe, the active layer thickness is an additional model parameter (noted $ELAY0$) which can be estimated by calibration and specified by the user in the steering file by keyword `ACTIVE LAYER THICKNESS'. With the option `CONSTANT ACTIVE LAYER THICKNESS' = NO, it is possible to use a time-and spatially varying active layer thickness during the simulation. In Sisyphe, we simply assume $ELAY0 = 3 D_m$ ($D_m$ is the mean diameter of the upper layer). The erosion rate is restricted by the total amount of sediments in the active layer, which therefore acts as a rigid bed. The same methods applied for rigid beds are applicable for active layer formulation. In order to avoid numerical problems such as negative sediment fractions, it is advised to use finite elements combined with the positive depth algorithm, as explained in \S 4.5 for the treatment of rigid bed, when dealing with graded sediment and a thin-sized active layer. The error message `TIME STEP SHOULD BE REDUCED' also appears in the listing file when the erosion is greater than half of the active layer thickness. It is stongly recommended to follow the message. === Continuous vertical sorting model === The C-VSM overcomes many limitations of the classic layer implementation of Hirano (HR-VSM). Even though it is a different way to manage the grain sorting, it is just another interpretation of Hiranos original idea with fewer simplifications. So still an active layer is defined. But the grain distribution in this layer is computed each time step from a depth dependent bookkeeping model for each grain size fraction with unlimited numerical resolution for each horizontal node. The new model doesn’t overcome the need to carefully calibrate the same input parameters as all other models, but the new interpretation has the following advantages: * It is possible to keep minor but prominent grain mixture variations even after a high number of time steps. Smearing effects and bookkeeping accuracy is defined by user defined thresholds or the computational resources, rather than through a fix value. * Various functions for the impact depth of the shear stress can be chosen to the projects demands. The result is a much more accurate vertical grain sorting, which results in better prognoses for bed roughness, bed forms and erosion stability. A couple of new keywords must be set in your sis.cas file. The full C-VSM output can be found in the new Selafin files VSPRES & VSPHYD in the tmp-folders. As the higher resolution of the C-VSM needs resources, you can reduce the printout period, or suppress the output at all. The common SISYPHE result files only show the averaged values for the active layer. Even more disk space can be saved, if only few points are printed out as .VSP.CSV files in the subfolder /VSP/. We recommend using between 200 and 1000 vertical sections. More will not improve the accuracy much, and less will lead to increasing data management, as the profile compression algorithms are called more often. __Keywords__\\ * 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 + 4x 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|1118|1750|2104|3316|1212|1280|2186|3187|1356|3027|1535|485 * 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 (Uses Keyword: ACTIVE LAYER THICKNESS) * 1 = Hunziker & Guenther \begin{equation*} ALT = 5 d_{MAX} \end{equation*} * 2 = Fredsoe & Deigaard (1992) \begin{equation*} ALT =\frac{2 \tau_B}{(1-n) g (\rho_S - \rho) \tan \Phi} \end{equation*} * 3 = van RIJN (1993) \begin{equation*} ALT =0.3 D_*^{0.7}\frac{\tau_B-\tau_C}{\tau_B}^{0.5} d_{50} \end{equation*} * 4 = Wong (2006) \begin{equation*} ALT =5 \frac{\tau_B}{(\rho_S - \rho) g d} - 0.0549)^{0.56} d_{50} \end{equation*} * 5 = Malcharek (2003) \begin{equation*} ALT =\frac{d_{90}}{1-n} \max(1,\frac{\tau_B}{\tau_C}) \end{equation*} * 6 = 3*d50 within last time steps ALT' \begin{equation*} ALT =3 d_{50} \end{equation*} === General formulation === Bed-load transport rates are computed separately for each class using classical sediment transport formulae, corrected for sand grading effects (see below). The Exner equation is solved for each class of material. The individual bed evolutions due to each class of bed material are then added to give the global evolution due to bed-load. Similarly, the suspended transport equation is solved for each class of material and the resulting bed evolutions for each class are added to give the global evolution due to the suspended load. The concentration of each class of sediments is computed, with the corresponding settling velocity, erosion and deposition flux. At each time step, the bed evolution (due to bed load and to suspension) is used to compute the new structure (composition and layer thickness) of the new sediment bed mixture. The algorithm which determines the new bed composition considers two possibilities (deposit or erosion) and ensures mass conservation for each individual class of material. The algorithm is explained in Gonzales (2002). === Hiding-exposure === In order to calculate the bed-load transport rate for each class of sediment, it is necessary to take into account the effect of hiding-exposure: in a mixture, the big particles will be more exposed to the flow than the smaller ones which will be protected. Classical sand transport formulae (for uniform sand) need to be corrected by use of a “hiding-exposure” correction factor (HIDFAC, as defined below in §6.4. __Keywords__\\ Sand grading effects are defined by: * HIDING FACTOR FORMULA (Hidfac=0, by default) * ACTIVE LAYER THICKNESS (Elay0= , by default) * CONSTANT ACTIVE LAYER THICKNESS (= YES, by default) * NUMBER OF BED LOAD MODEL LAYERS (NOMBAY= 2, by default) {{:panneau_attention_40.png |}}**Remarks** \\ \\ * The number of layers has to be less than 10. The user has to compare the bed layers thickness in comparison to the total bed thickness (given in noerod.f) to ensure that there will be no more than 10 layers during the computation. The composition of the sediment AVAIL is defined at each point and for each layer: AVAIL(NPOIN,10,NSICLA). The sum of AVAIL on the classes must be equal to one everywhere. * The number of layers is allowed to vary spatially and with time as the sediment bed undergoes erosion (empty sediment layers are suppressed) or deposition (sediment layers can be added). * A new procedure will be developed in release 6.1 with a constant number of layers which will be compatible with the consolidation algorithm. === Sediment transport formulae === The formula of Hunziker ([24], 1995) is an adaptation of the Meyer-Peter Mu"ller formula for fractional transport. The volumetric sediment transport per class is given by: \begin{equation*} \Phi_b = 5 p_i \left[\zeta_i\left( \mu \theta_i -\theta_{cm} \right) \right]^{3/2} \quad\text{if}\quad \mu \theta_i > \theta_{cr}, \end{equation*} with $p_i$ the fraction of class $i$ in the active layer, $\theta_i$ the Shields parameter,$\theta_{cm} = \theta_{cr}\left(D_{mo}/D_m\right)^{0.33}$ the corrected critical Shields parameter, $\theta_{cr}$ the critical Shields parameter ($\theta_{cr}=0.047$), $\xi_{i} =\left(D_i/D_m\right)^{-\alpha}$ the hiding / exposure factor, $D_i$ the grain size of class $i$, $D_m$ the mean grain size of the surface layer (m), and $D_{mo}$ the mean grain size of the under layer (m). Here the critical Shields parameter is calculated, as a function of the dimensionless mean grain size $D_*$. It should be noted that according to Hunziker, stability problems may occur outside the parameter range: \begin{equation*} \alpha \leq 2.3\quad\text{and}\quad \frac{D_i}{D_m} \geq 0.25. \end{equation*} __Keywords__\\ Sand grading effects are defined by: * BED-LOAD TRANSPORT FORMULA (ICF = 6) * HIDING FACTOR FORMULA (not used if ICF = 6) ==== Hiding/exposure factors ==== Sediment transport rates of each class can be calculated using classical sediment transport formulae. Those formulation are initially valid for uniform particles and need to be modified to account for sand grading effects: in a mixture, big particles are more exposed to the flow than if they are alone, and small particles more protected. The bed-load formulae can be modified by use of a hiding/exposure factor. Different formulations have been implemented in SISYPHE, depending on the choice of ‘HIDING FACTOR FORMULA’ (HIDFAC). Two formulae (Egiazaroff and Ashida & Michiue) have been written in association with the Meyer-Peter and Müller formula. The third formula (Karim et al)can be associated with any bed-load predictor. The first two formulae modify the critical Shields parameter that will be used in the Meyer-Peter formula, whereas the third formula modify directly the bed-load transport rate. * Formulation of Egiazaroff (HIDFAC=1) \begin{equation*} \begin{align*} \theta_{cr} &= \zeta_i 0.047, \\ \zeta_i &= \left[\frac{\log (19)}{\log (19 D_i /D_m)} \right]^2 \end{align*} \end{equation*} * Formulation of Ashida & Michiue (HIDFAC=2) * Formulation of Karim et al. (HIDFAC=4)