Steady 2d

Requirements

This tutorial is designed for advanced beginners and before diving into this tutorial make sure to complete the TELEMAC pre-processing tutorial.

The case featured in this tutorial was established with the following software:

Get Started

This section builds on the SELAFIN (*.slf) geometry and the Conlim (*.cli) boundary condition files that result from the TELEMAC pre-processing tutorial. Both files can also be downloaded from the supplemental materials repository of this eBook:

Consider saving both files in a new folder, such as /steady2d-tutorial/ that will contain all model files.

Download simulation files

All simulation files used in this tutorial are available at https://github.com/hydro-informatics/telemac/tree/main/steady2d-tutorial/.

Steering File (CAS)

The steering file has the file ending *.cas (presumably derived from the French word cas, which means case in English). The *.cas file is the main simulation file with information about references to the two always mandatory files (i.e., the SELAFIN *.slf geometry and the *.cli boundary files) and optional files, as well as definitions of simulation parameters. The steering file can be created or edited with a basic text editor or advanced GUI software such as Fudaa PrePro or BlueKenue. This tutorial uses a basic text editor (e.g., Notepad++ on Windows).

Fudaa PrePro

Fudaa PrePro comes with variable descriptions that facilitate the definition of boundaries, initial conditions, and numerical parameters in the steering file. However, Fudaa PrePro makes file directions according to the platform on which it is running (i.e., \ on Windows and / on Linux) and writes absolute file paths to the *.cas file, which often requires manual correction (e.g., if Fudaa PrePro is used for setting up a *.cas file on Windows for running a TELEMAC simulation on Linux). For working with Fudaa PrePro, follow the download instructions in the software chapter. To launch Fudaa Prepro, open Terminal (Linux) or Command Prompt (Windows) and tap:

  • cd to the installation (download) directory of Fudaa PrePro

  • Start the GUI (requires java):

    • Linux: sh supervisor.sh

    • Windows: supervisor.bat

For this tutorial, create a new text file in the same folder where qgismesh.slf and boundaries.cli live, and name it, for instance, steady2d.cas (e.g., /steady2d-tutorial/steady2d.cas). The next sections guide through parameter definitions that stem from the Telemac2d manual. The final steering file can be downloaded from the supplemental materials repository (download steady2d.cas).

Overview of the CAS File

The below box shows the provided steady2d.cas file that can be used for running this tutorial.

General Parameters

The general parameters define the computation environment starting with a simulation title and the most important links to the two mandatory input files:

  • BOUNDARY CONDITIONS FILE : boundaries.cli - with a MED file, use a BND boundary file

  • GEOMETRY FILE : qgismesh.slf

The model output can be defined with the following keywords:

  • RESULTS FILE : r2dsteady.slf - can be either a MED file or a SLF file

  • VARIABLES FOR GRAPHIC PRINTOUTS: U,V,H,S,Q,F - many more options can be found in section 1.317 (page 85) of the Telemac2d reference manual.

The velocities (U and V), the water depth (H), and the discharge (Q) are standard variables that should be used in every simulation. In particular, the discharge (Q) is required to check when (steady) s converge at the inflow and outflow boundaries. Moreover, the discharge (Q) enables to trace integrated fluxes along any user-defined line in the model. The procedure for verifying and identify discharges is described in the discharge verification section in the post-processing.

The time variables (TIME STEP and NUMBER OF TIME STEPS) define the simulation length and the printout periods (GRAPHIC PRINTOUT PERIOD and LISTING PRINTOUT PERIOD) define the result output frequency. The smaller the printout period, the longer will take the simulation because writing results is one of the most time-consuming processes in numerical modeling. The printout periods (frequencies) refer to a multiple of the TIME STEPS parameter and need to be a smaller number than the NUMBER OF TIME STEPS. Read more about timestep parameters in the Telemac2d manual in sections 5 and 12.4.2.

In addition, the MASS-BALANCE : YES setting will print out mass fluxes and errors in the computation region, which is an important parameter for verifying the plausibility of the model. Note that this keyword only enables mass balance printouts and does not enforce mass balance of the model, which must be achieved through a consistent model setup following this tutorial and the Telemac2d manual.

General Numerical Parameters

The following descriptions refer to section 7.1 in the Telemac2d manual.

Telemac2d comes with three solvers for approximating the depth-averaged Navier-Stokes equations (i.e., the Shallow water equations) [KC08] that can be chosen by adding the EQUATIONS keyword to the *.cas file:

  • EQUATIONS : SAINT-VENANT FE is the default that makes Telemac2d use a Saint-Venant finite element method,

  • EQUATIONS : SAINT-VENANT FV makes Telemac2d use a Saint-Venant finite volume method, and

  • EQUATIONS : BOUSSINESQ makes Telemac2d use the Boussinesq approximations (constant density except in the vertical momentum equation).

In addition, a type of discretization has to be specified with the DISCRETIZATIONS IN SPACE keyword, which is a list of five integer values. The five list elements define spatial discretization schemes for (1) velocity, (2) depth, (3) tracers, (4) \(k-\epsilon\) turbulence, and (5) \(\tilde{\nu}\) advection (Spalart-Allmaras), respectively. The minimum length of the keyword list is 2 (for velocity and depth) and all other elements are optional. The list elements may take the following values defining spatial discretization:

  • 11 (default) activates triangular discretization in space (i.e., 3-node triangles),

  • 12 activates quasi-bubble discretization with 4-node triangles, and

  • 13 activates quadratic discretization with 6-node triangles.

The Telemac2d manual recommend using the default value of DISCRETIZATIONS IN SPACE : 11;11 that assigns a linear discretization for velocity and water depth (computationally fastest). The option 12;11 may be used to reduce free surface instabilities or oscillations (e.g., along with steep bathymetry gradients). The option 13;11 increases the accuracy of results, the computation time, memory usage, and it is currently not available in Telemac2d.

In addition, the FREE SURFACE GRADIENT keyword can be defined for increasing the stability of a model. Its default value is 1.0, but it can be reduced close to zero to achieve stability. The developers propose a minimum value of 0., but this would lead to non-meaningful results, which is why this eBook recommends a value slightly larger than zero. For instance, the following keyword combination may reduce surface instabilities (also referred to as wiggles or oscillations):

DISCRETIZATIONS IN SPACE : 12;11
FREE SURFACE GRADIENT : 0.03

By default Advection is activated through the keyword ADVECTION : YES and it can be deactivated for particular terms only:

ADVECTION OF H : NO / deactivates depth advection
ADVECTION OF U AND V : NO / deactivates velocity advection
ADVECTION OF K AND EPSILON : NO / deactivates turbulent energy and dissipation (k-e model) or the Spalart-Allmaras advection
ADVECTION OF TRACERS : NO / deactivates tracer advection

The PROPAGATION keyword (default: YES) steers the simulation of propagation and related phenomena. For instance, disabling propagation (PROPAGATION : NO) will also disable Diffusion. The other way round, when propagation is enabled, Diffusion can be disabled separately. Read more about Diffusion in Telemac2d in the turbulence section.

Numerical Parameters for Finite Elements

The following descriptions refer to section 7.2.1 in the Telemac2d manual.

Telemac2d uses finite elements for iterative solutions to the Shallow water equations. The TREATMENT OF THE LINEAR SYSTEM keyword enables replacing the original set of equations (option 1) involved in TELEMAC’s finite element solver with a generalized wave equation (option 2). The replacement (i.e., the use of the generalized wave equation) is set to default since v8p2 and decreases computation time, but smoothens the results. This default (TREATMENT OF THE LINEAR SYSTEM : 2) automatically activates mass lumping for depth and velocity, and implies explicit velocity diffusion.

The Telemac2d manual state that the following scalar SCHEME FOR ADVECTION keywords apply instead of the soon deprecated TYPE OF ADVECTION list:

SCHEME FOR ADVECTION OF VELOCITIES : 1 / default
SCHEME FOR ADVECTION OF TRACERS : 1 / default
SCHEME FOR ADVECTION OF K-EPSILON : 1 / default

The three SCHEME FOR ADVECTION scalar keywords may take the following values:

  • 1 sets a not mass-conservative method of characteristics (default for all),

  • 2 sets a semi-implicit scheme and activates the Streamline Upwind Petrov Galerkin (SUPG) scheme (read more below),

  • 3, 4, 13, and 14 activate the so-called NERD scheme (these numbers activate different schemes in 3d only),

  • 5 sets a mass-conservative PSI distributive scheme (do not use with tidal flats), and

  • 15 sets the mass-conservative ERIA scheme that works with tidal flats.

Options 4 and 5 require that the CFL condition is smaller than 1.

Recommended SCHEME OF ADVECTION … keywords

The Telemac2d manual recommend specific combinations depending on the simulation scenario.

For models without any dry zones use:

SCHEME FOR ADVECTION OF VELOCITIES : 4 / alternatively keep 1
SCHEME FOR ADVECTION OF TRACERS : 5
SCHEME FOR ADVECTION OF K-EPSILON : 4

For models with tidal flats use (like in this tutorial):

SCHEME FOR ADVECTION OF VELOCITIES : 14 / alternatively keep 1
SCHEME FOR ADVECTION OF TRACERS : 5
SCHEME FOR ADVECTION OF K-EPSILON : 14

Without any SCHEME FOR ADVECTION … keyword, the SUPG OPTION (Streamline Upwind Petrov Galerkin) keyword defines if upwinding applies and what type of upwinding applies. The SUPG OPTION is a list of four integers, where every element may take one of the following values:

  • 0 disables upwinding,

  • 1 enables upwinding with a classical SUPG scheme (recommended when the CFL condition is unknown), and

  • 2 enables upwinding with a modified SUPG scheme, where upwinding equals the CFL condition (recommended when the CFL condition is small).

The default is SUPG OPTION : 2;2;2;2, where

  • the first list element refers to flow velocity (default 2),

  • the second to water depth (default 2 - set to 0 when MATRIX STORAGE : 3),

  • the third to tracers (default 2), and

  • the last (fourth) to the k-epsilon model (default 2).

Note that the SUPG OPTION keyword is not optional for many keyword combinations and this tutorial uses SUPG OPTION : 0;0;2;2.

Implicitation parameters (IMPLICITATION FOR DEPTH, IMPLICITATION FOR VELOCITIES, and IMPLICITATION FOR DIFFUSION OF VELOCITY) apply to the semi-implicit time discretization used in Telemac2d. To enable cross-version compatibility, implicitation parameters should be defined in the *.cas file. For DEPTH and VELOCITIES use values between 0.55 and 0.60 (default is 0.55 since v8p1); for IMPLICITATION FOR DIFFUSION OF VELOCITY set the v8p2 default of 1.0.

The default TREATMENT OF THE LINEAR SYSTEM : 2 involves so-called mass lumping, which leads to a smoothening of results. Specific mass lumping keywords and values are required for the flux control option of the TREATMENT OF NEGATIVE DEPTHS keyword and the default value for the treatment of tidal flats. To this end, the mass lumping keywords should be defined as:

MASS-LUMPING ON H : 1.
MASS-LUMPING ON VELOCITY : 1.
MASS-LUMPING ON TRACERS : 1.

In addition, MASS-LUMPING FOR WEAK CHARACTERISTICS : 1. may be defined, which will make Telemac2d using weak characteristics (see below). The default value of any MASS-LUMPING ... keyword is 0. and the maximum value is 1., which makes mass matrices diagonal.

The OPTION OF CHARACTERISTICS keyword defines the method of characteristics that can take a strong (default of 1) or a weak (2) form. A weak form decreases Diffusion, is more conservative, and increases computation time. Telemac2d automatically switches from the default strong (1) to the weak (2) form when

  • the TYPE OF ADVECTION is set to 1,

  • any SCHEME FOR ADVECTION ... is set to 1, or

  • any SCHEME OPTION FOR ADVECTION OF ... is set to 2.

None of these options should be used with tracers because they are not mass-conservative.

Numerical Parameters for Finite Volumes

The finite volume method is mentioned here for completeness and detailed descriptions are available in section 7.2.2 of the Telemac2d manual.

Use finite volumes only with v8p2 or later

Earlier versions of Telemac2d’s finite volume solver are buggy, but since major improvements were implemented wit v8p2, the newest versions run stable.

The finite volume method involves the definition of a scheme through the FINITE VOLUME SCHEME keyword that can take one of the following integer values:

  • 0 enables the Roe [Roe81] scheme,

  • 1 is the default and enables the kinetic scheme [ABP00],

  • 3 enables the Zokagoa and Soulaïmani [ZS10] scheme that is incompatible with tidal flats,

  • 4 enables the Tchamen and Kahawita [TK98] scheme for modeling wetting and drying of a complex bathymetry,

  • 5 enables the Harten Lax Leer-Contact (HLLC) scheme [Tor09], and

  • 6 enables the Weighted Average Flux (WAF) [Ata12] scheme for which parallelism is currently not implemented.

All finite volume schemes are explicit and potentially subjected to instability. For this reason, a desired CFL condition and a variable timestep are recommended to be defined:

DESIRED COURANT NUMBER : 0.9
VARIABLE TIME-STEP : YES / default is NO
DURATION : 10000

The DURATION keyword is required to terminate the simulation.

The variable timestep will cause irregular listing outputs, while the graphic output frequency is written as a function of the above-defined TIME STEP. Note that this tutorial uses VARIABLE TIME-STEP : NO, though variable timesteps do also work with FINITE ELEMENTS.

The FINITE VOLUME SCHEME TIME ORDER keyword defines the second-order time scheme, which is by default set to Euler explicit (1). Setting the time scheme order to 2 makes Telemac2d using the Newmark scheme where an integration coefficient may be used to change the integration parameter. Note that NEWMARK TIME INTEGRATION COEFFICIENT : 1 corresponds to Euler explicit. To implement these options in the steering file, use the following settings:

FINITE VOLUME SCHEME TIME ORDER : 2 / default is 1 - Euler explicit
NEWMARK TIME INTEGRATION COEFFICIENT : 0.5 / default is 0.5

Depending on the type of analysis, the solver-related parameters of SOLVER, SOLVER OPTIONS, MAXIMUM NUMBER OF ITERATION FOR SOLVER, and TIDAL FLATS may also be modified.

Numerical Solver Parameters

The following descriptions refer to section 7.3.1 in the Telemac2d manual.

The solver can be selected and specified with the SOLVER, SOLVER FOR DIFFUSION OF TRACERS, and SOLVER FOR K-EPSILON MODEL keywords where the following settings are recommended:

SOLVER : 1 / default is 3
SOLVER FOR DIFFUSION OF TRACERS : 1
SOLVER FOR K-EPSILON MODEL : 1

Setting the SOLVER to 1 instead of the default value of 3 is recommended with TREATMENT OF THE LINEAR SYSTEM : 2 (i.e., the default since v8p2) for consistent and backward-compatible steering files.

Every solver keyword can take an integer value between 1 and 8, where 1-6 use conjugate gradient methods:

  • 1 sets the conjugate gradient method for symmetric matrices,

  • 2 sets the conjugate residual method,

  • 3 sets the conjugate gradient on normal equation method,

  • 4 sets the minimum error method,

  • 5 sets the squared conjugate gradient method,

  • 6 sets the stabilized biconjugate gradient (BICGSTAB) method,

  • 7 sets the Generalised Minimum RESidual (GMRES) method, and

  • 8 sets the Yale University direct solver (YSMP) that is not compatible with parallelism.

The GMRES method may be enabled with the finite element scheme where the following solver options for the Krylov space:

SOLVER OPTION : 2 / hydrodynamic propagation
SOLVER OPTION FOR TRACERS DIFFUSION : 2 / tracer diffusion
OPTION FOR THE SOLVER FOR K-EPSILON MODEL : 2 /  k-e or Spalart-Allmaras

The solver options vary between values of 2 for a small mesh and 5 for a large mesh. Integers of 3 or 4 can be used for medium-sized meshes. The Telemac2d manual recommends running simulations multiple times for finding an optimum value, where higher values (close to 5) increase the time required for an iteration but lead to faster convergence.

Numerical Accuracy

The following descriptions refer to section 7.3.2 in the Telemac2d manual.

The accuracy keywords make Telemac2d stop an iteration when two consecutive solutions for the same element vary by less than an ACCURACY threshold. To this end, the following default accuracy thresholds may be varied (Telemac2d ignores non-relevant parameters):

SOLVER ACCURACY : 1.E-4 / propagation steps
ACCURACY FOR DIFFUSION OF TRACERS : 1.E-6 / tracer diffusion
ACCURACY OF K : 1.E-9 / diffusion and source terms of turbulent energy transport
ACCURACY OF EPSILON : 1.E-9 / diffusion and source terms of turbulent dissipation transport
ACCURACY OF SPALART-ALLMARAS : 1.E-9 / diffusion and source terms of the Spalart-Allmaras equation

In experience, the solver accuracy should not be larger than 1.E-3 (10\(^{-3}\)). In contrast, very small accuracies will lead to longer computation times. In addition or alternatively to the accuracy keywords, the following default numbers of maximum iterations can be modified to speed up calculations:

MAXIMUM NUMBER OF ITERATIONS FOR SOLVER : 100 / maximum number of iterations when solving the propagation step
MAXIMUM NUMBER OF ITERATIONS FOR DIFFUSION OF TRACERS : 60 / tracer diffusion
MAXIMUM NUMBER OF ITERATIONS FOR K AND EPSILON : 50 / diffusion and source terms of k-e or Spalart-Allmaras

Telemac2d will print out warning messages when convergence could not be reached with the defined combination of accuracy and maximum iteration number keywords. The warning message printouts can be deactivated with the INFORMATION ABOUT SOLVER keyword, though deactivating convergence warnings is not recommended.

Tidal Flats

The following descriptions refer to section 7.5 in the Telemac2d manual.

The TIDAL FLATS (default: YES) keyword applies to the finite elements scheme only (EQUATIONS keyword) and can be ignored with finite volumes. The term tidal may be slightly confusing because tidal flats can occur beyond coastal regions: Tidal flats can occur wherever there are flow transitions, such as when fast-flowing water enters a backwater zone. Flow transitions occur in almost all environments more complex than a square-like flume, and therefore, the activation of tidal flats in Telemac2d models is highly recommended. Though activating tidal flats leads to longer computation times, in most cases a calculation with tidal flats provides physically reasonable results.

The TIDAL FLATS keyword is linked with a couple of other Telemac2d keywords driving model stability and physical meaningfulness. The following keyword setups may be generally applied to (quasi) steady, real-world rivers and channels (as opposed to lab flumes with simplified geometries):

TIDAL FLATS : YES
CONTINUITY CORRECTION : YES / default is NO
OPTION FOR THE TREATMENT OF TIDAL FLATS : 1
TREATMENT OF NEGATIVE DEPTHS : 2 / value 2 or 3 is required with tidal flats

The OPTION FOR THE TREATMENT OF TIDAL FLATS accepts integer values between 1 and 3 to select one of the following options:

  • 1 detects tidal flats and corrects the free surface gradient.

  • 2 removes tidal flat elements by using a masking table that eliminates any contribution of concerned mesh elements. This option may affect the mass conservation of the model.

  • 3 resembles 1, but adds a porosity term to half-dry mesh elements. This affects the amount of water in the model, which equals here the depth integral multiplied by the porosity. A user Fortran file may be used to modify the porosity term in the USER_CORPOR subroutine.

The TREATMENT OF NEGATIVE DEPTHS (default: 1) keyword defines an approach for eliminating negative water depth values where the following integer numbers can be used:

  • 0 disables any treatment of negative water depths.

  • 1 conservatively smoothens negative water depths (default).

    • A float number keyword THRESHOLD FOR NEGATIVE DEPTHS (default: 0.) is available only for this option.

    • Setting the threshold to, for instance, -0.1 makes that negative water depths larger (e.g., -0.05 m) than -0.1 meters remain unchanged.

  • 2 imposes a flux limitation that strictly ensures positive water depths.

  • 3 acts similarly as 2 but for the ERIA Advection scheme (set SCHEME FOR ADVECTION OF TRACERS to 4 or 5). This option is appropriate for modeling conservative tracers.

TIDAL FLATS options require particular keyword combinations

The SCHEME FOR ADVECTION … keywords (see the finite element parameters section) must be set for TRACERS to LIPS (either 4 or 5), and for all others to either a NERD (13 or 14) or the ERIA (15) scheme.

When using LIPS (4 or 5) with NERD (13 or 14) use the following combination (used in this tutorial):

TIDAL FLATS : YES
OPTION FOR THE TREATMENT OF TIDAL FLATS : 1
TREATMENT OF NEGATIVE DEPTHS : 2

When using LIPS (4 or 5) with ERIA (15) use the following combination:

TIDAL FLATS : YES
OPTION FOR THE TREATMENT OF TIDAL FLATS : 1
TREATMENT OF NEGATIVE DEPTHS : 3

Read more about viable or trouble-making parameter combinations for tidal flats in section 16.5 in the Telemac2d manual.

Matrix Handling

The following descriptions refer to section 7.6 in the Telemac2d manual.

Telemac2d provides multiple options for matrix handling that need to be set up for particular solver schemes.

The MATRIX STORAGE keyword may be set to:

  • 1 for using classic element-by-element matrix storage.

  • 3 for using edge-based matrix storage (default). This default is required when any SCHEME FOR ADVECTION … keyword is set to 3, 4, 5, 13, 14, or 15, and when any direct SOLVER is set to 8.

The additional MATRIX-VECTOR PRODUCT keyword may be used to switch between multiplication methods for the finite element scheme. However, the default value of 1 (vector multiplication by a non-assembled matrix) should currently not be changed because the only alternative (2 for frontal assembled matrix multiplication) is not implemented for parallelism and quasi-bubble discretization.

Friction Boundary Conditions

The following descriptions of friction parameters refer to section 6.1 in the Telemac2d manual.

Parameters for Boundary Conditions enable the definition of roughness laws and properties of liquid boundaries.

The LAW OF BOTTOM FRICTION keyword defines a friction law for topographic boundaries, which can be set to:

  • 0 for no friction.

  • 1 for the Haaland [Haa83] equation, which is an implicit form of the Colebrook and White [CW37] equation that builds on the Darcy-Weisbach friction factor \(f_D\). This law involves a high degree of uncertainty that stems from the experimental dataset of the original author.

  • 2 for the Chézy [Che76] roughness that can be similarly used as 3 and 4.

  • 3 for Strickler [Str23] roughness \(k_{st}\) (read more, for example, in the 1d hydraulics exercise), which is the inverse of \(n_m\) (4).

  • 4 for Manning [Man91] roughness \(n_m\) (read more, for example, in the 1d hydraulics exercise), which is the inverse of \(k_{st}\) (3).

  • 5 for the Nikuradse [Nik33] roughness law, which should correspond to 3 \(\cdot D_{90}\) according to van Rijn [vR19].

  • 6 for the logarithmic law of the wall for turbulent flows. This option assumes that the average flow velocity is a logarithmic function of the distance from the wall beyond the viscous and buffer layers. The thickness of these layers is a function of the wall roughness length [VK30].

  • 7 for the Colebrook and White [CW37] equation that calculates the Darcy-Weisbach friction factor \(f_D\) for turbulent flows in smooth pipes.

With respect to the 2d applications in this eBook, the most relevant bottom friction laws are 3 [Str23], 4 [Man91], and 6 (log law). The Nikuradse [Nik33] roughness law (5) is recommended for 3d simulations (see the Telemac3d tutorial). Moreover, friction is more generally referred to with the general coefficient \(c_{f}\), which has a particular relevance for Bedload transport (cf. morphodynamic calculations with Gaia).

The FRICTION COEFFICIENT FOR THE BOTTOM keyword sets the value for a characteristic roughness coefficient. For instance, when the friction law keyword is set to 3 [Str23], the friction corresponds to the Strickler roughness coefficient \(k_{st}\) (in fictive units of m\(^{1/3}\) s\(^{-1}\)). For rough channels (e.g., mountain rivers) \(k_{st} \approx 20\) m\(^{1/3}\) s\(^{-1}\) and for smooth concrete-lined channels \(k_{st} \approx 75\) m\(^{1/3}\) s\(^{-1}\). In fully turbulent flows, the Strickler roughness can be approximated with \(k_{st} \approx \frac{26}{D_{90}^{1/6}}\) [MPM48] where \(D_{90}\) is the grain diameter of which 90% of the surface grain mixture are finer. This tutorial features the application of a Manning roughness coefficient of \(n_m\)= 0.03, which is the inverse of \(k_{st}\) and implemented with:

LAW OF BOTTOM FRICTION : 4 / 4-Manning
FRICTION COEFFICIENT : 0.03 / Roughness coefficient

Liquid Boundary Conditions

The following descriptions of friction parameters refer to section 4.2 in the Telemac2d manual.

Liquid boundary keywords assign hydraulic properties to the spatially defined upstream and downstream liquid boundary lines in the Conlim (*.cli) file created with BlueKenue. This section features the assignment of steady liquid boundaries for one discharge of 35 m\(^3\)/s. To this end, the upstream boundary condition is set to a steady target inflow rate (Open boundary with prescribed Q) and the downstream boundary condition gets a Stage-discharge relation (Open boundary with prescribed Q and H) assigned (recall Fig. 74). Thus, for running this tutorial add the following keywords to the steering (*.cas) file:

  • The keyword PRESCRIBED FLOWRATES : 35.;35. assigns a flowrate of 35 m\(^3\)/s to the upstream and the downstream boundary edges, respectively.

  • The keyword PRESCRIBED ELEVATIONS : 0.;371.33 assigns a water surface elevation \(wse\) (or H in Telemac) of 371.33 m a.s.l. (above sea level) to the downstream boundary. The upstream \(wse\) prescription of 0.0 makes Telemac2d ignore this value corresponding to the assigned upstream boundary type (prescribed flowrate only).

The order of prescribed flowrates (Q) and \(wse\) (H) values depends on the order of the definition of the boundaries. Thus, the first list element defines values for the upstream and the second list element for the downstream open boundary.

How to find out the order of boundary conditions?

The order of open boundaries can be read from the *.cli file. The first open boundary that is listed in the *.cli file corresponds to the first list element in any PRESCRIBED … keyword. An open boundary node in the *.cli file is characterized by a line beginning with something like 4 5 5 or 5 5 5 (i.e., anything but 2 2 2, which corresponds to a closed wall boundary node) and BlueKenue also marks the names of open boundaries at the line ends (after the hashtag). Figure 78 illustrates the boundaries.cli file used in this chapter where the upstream open boundary is defined at line 7, before the definition of the downstream open boundary starting at line 313.

telemac 2d cli boundary conditions order cas steering file prescribed prescription

Fig. 78 The boundaries.cli file used in this tutorial starts with the upstream boundary defined in line 7. To find the downstream boundary scroll down to line 313.

Liquid boundary conditions may be assigned to any open boundary in the *.cli file.

A velocity profile type can be assigned to any prescribed Q (flowrate) or prescribed U (velocity) open boundary in the form of a list that has the same element order as the above-defined PRESCRIBED … keywords. For this purpose, upstream and downstream velocity profiles can be defined with the VELOCITY PROFILES keyword that accepts the following values:

  • 1 is the default option that defines the flow velocity direction at the boundary nodes normal to their edges. This option assigns a length of 1 to the vector and multiplies it with a numeric factor to yield a target flowrate.

  • 2 reads U and V velocity profiles from the boundary conditions (*.cli) file, which are multiplied with a constant to yield a target flowrate.

  • 3 imposes the velocity vector direction normal to the boundary and reads the value (UBOR) from the *.cli file, which is then multiplied by a constant to yield a target flowrate.

  • 4 imposes the velocity vector direction normal to the boundary and calculates the value’s norm proportional to the square root of the water depth. This option can only be used with a prescribed Q open boundary.

  • 5 imposes the velocity vector direction normal to the boundary and calculates the value’s norm proportional to the square root of a virtual water depth.

With the upstream boundary being a prescribed Q boundary, this tutorial uses VELOCITY PROFILES : 4;1 in the steering file. Read more about options for defining velocity profiles in section 4.2.8 of the Telemac2d manual.

Initial Conditions

The following descriptions refer to section 4.1 in the Telemac2d manual.

The initial conditions describe the state of the model at the beginning of a simulation. Telemac2d recognizes the following types of initial conditions, which can be defined in the steering file with the keyword INITIAL CONDITIONS : 'TYPE' where TYPE can be one of the following:

  • ZERO ELEVATION initializes the free surface elevation at 0 (default). Thus, the initial water depths correspond to the bottom elevation.

  • CONSTANT ELEVATION initializes the free surface elevation at a value defined with an INITIAL ELEVATION keyword that has a default value of 0.. Thus, the initial water depths correspond to the subtraction of the bottom elevation from the water surface elevation \(wse\). The initial water depth is set to zero at nodes where the bottom elevation is higher than defined by the INITIAL ELEVATION keyword.

  • ZERO DEPTH initializes the simulation with 0 (i.e., \(wse\) corresponds to bottom elevation). Thus, the model starts with dry conditions, similar as in the BASEMENT tutorial.

  • CONSTANT DEPTH initializes the water depths at a value defined by an INITIAL DEPTH keyword that has a default value of 0..

  • TPXO SATELLITE ALTIMETRY initializes the model using information provided by a user-defined database (e.g., the OSU TPXO model for ocean tides). Read more in section 4.2.12 of the Telemac2d manual on modeling marine systems.

This tutorial uses a constant water depth initial condition of 1 (integer to speed up calculations), which corresponds to a completely flooded initial model state (i.e., water volume surplus):

INITIAL CONDITIONS : 'CONSTANT DEPTH'
INITIAL DEPTH : 1

The simulation speed can be significantly increased when the model has already been running once at the same discharge. The result of an earlier simulation can be used for the initial condition with the COMPUTATION CONTINUED : YES (default is NO) and PREVIOUS COMPUTATION FILE : *.slf (provide the name of a *.slf file) keywords. This type of model initialization is also referred to as hotstart. Read more about hotstarts in the unsteady simulation and Gaia sections. Also section 4.1.3 in the Telemac2d manual provides descriptions for continuing (hotstart) calculations.

Turbulence

The following descriptions refer to section 6.2 in the Telemac2d manual.

Turbulence describes a seemingly random and chaotic state of fluid motion in the form of three-dimensional vortices (eddies). True turbulence is only present in 3d vorticity and when it occurs, it mostly dominates all other flow phenomena through increases in energy dissipation, drag, heat transfer, and mixing [KC08]. The phenomenon of turbulence has been a mystery to science for a long time, since turbulent flows have been observed, but could not be explained by the linear equations systems. Today, turbulence is considered a random phenomenon that can be accounted for in linear equations, for instance, by introducing statistical parameters. For instance, when turbulence applies to the depth-averaged Navier-Stokes equations a numerical solution for a quantity (e.g., flow velocity) corresponds to \(value = \overline{mean value} + value fluctuation'\). For this purpose, there are a variety of options for implementing turbulence in numerical models [NN93].

The horizontal and vertical dimensions of turbulent eddies can vary greatly, especially in rivers and transitions to backwater zones (tidal flats) where the wide horizontal flow dimension (river width \(w\)) is significantly larger than the vertical flow dimension (water depth \(h\)): \(w >> h\). Telemac2d provides multiple turbulence models that can be applied to the vertical and/or horizontal dimensions and defined with the TURBULENCE MODEL keyword being an integer number for one of the following options:

  • 1 to use a constant viscosity coefficient (default) for turbulent viscosity, molecular viscosity, and Diffusion. This closure option should not be used with Stage-discharge relation open boundaries (i.e., do not use with prescribed Q and H) [WBH02].

  • 2 to use the Elder formula for the Diffusion coefficient \(D\). The Elder turbulence closure also yields small errors for Stage-discharge relation open boundaries (i.e., do not use this option with prescribed Q and H) [WBH02].

  • 3 to use the \(k-\epsilon\) two-equation model solving the Navier-Stokes equations. The first equation represents a turbulence closure for the turbulent energy \(k\); the second equation is a turbulence closure for the turbulent dissipation \(\epsilon\). Both equations express that the sum of change of (I) \(k\) and \(\epsilon\) in time, and (II) Advection transport of \(k\) and \(\epsilon\) equal the sum of (1) Diffusion transport of \(k\) and \(\epsilon\), (2) the production rate of \(k\)/\(\epsilon\), and (3) the destruction rate of \(k\)/\(\epsilon\) [LS74]. The \(k-\epsilon\) model is a generalization of the mixing length model (see option 5) and assumes that the turbulent viscosity is isotropic (valid for many river applications, but not for circular-rotating flows or groundwater) [Bra87]. Thus, the \(k-\epsilon\) model introduces two additional equations and requires a finer mesh than the constant viscosity option 1, which leads to a longer computation time. Yet, the \(k-\epsilon\) model generally yields accurate results and small errors with Stage-discharge relation open boundaries [WBH02]. The following default keywords are associated with the \(k-\epsilon\) model:

    • VELOCITY DIFFUSIVITY : 1.E-6 corresponding to the kinematic viscosity \(\nu\) of water (10\(^{-6}\) m\(^2\)/s).

    • TURBULENCE REGIME FOR SOLID BOUNDARIES : 2 for rough walls of closed boundaries to apply the value chosen for the LAW OF BOTTOM FRICTION and ROUGHNESS COEFFICIENT OF BOUNDARIES keywords (recall section Friction Boundary Conditions). For smooth closed boundary walls set TURBULENCE REGIME FOR SOLID BOUNDARIES : 1.

    • INFORMATION ABOUT K-EPSILON MODEL : YES enables console output of information on the \(k-\epsilon\) closure solution.

  • 4 to use the Smagorinsky [Sma63] (also known as general circulation) model, which stems from climate modeling and is appropriate for modeling maritime systems with large eddies. The Smagorinsky [Sma63] model does not account for Diffusion.

  • 5 to use a mixing length model according to Prandtl’s theory that a fluid quantity conserves its properties for a characteristic length before it mixes with the bulk flow [Bra74].

  • 6 to use the Spalart and Allmaras [SA92] model that solves the Continuity equation for a viscosity-like, kinematic eddy turbulent viscosity. The Spalart and Allmaras [SA92] model was originally developed for aerodynamic flows with low Reynolds number and it has also shown good results for other applications.

This tutorial uses the \(k-\epsilon\) model (3) because of its low error rate and wide applicability (compared to other turbulence closures).

DIFFUSION OF VELOCITY : YES / enabled by default
TURBULENCE MODEL : 3

Run Telemac2d

With the steering (*.cas) file, the last necessary ingredient for running a steady hydrodynamic 2d simulation with Telemac2d is available. Make sure to put all required files in one simulation folder (e.g., ~/telemac/v8p2/mysimulations/steady2d-tutorial/). The required files can also be downloaded from this eBook’s steady2d tutorial repository and include:

With these files prepared, load the TELEMAC environment, and run Telemac2d following the explanations in the next sections.

Load environment and files

Go to the configuration folder of the TELEMAC installation (e.g., ~/telemac/v8p2/configs/) and load the environment (e.g., pysource.openmpi.sh - use the same as for compiling TELEMAC).

cd ~/telemac/v8p2/configs
source pysource.openmpi.sh

To start a simulation, change to the directory (cd) where the simulation files live and run the steering file (*.cas) with the telemac2d.py script:

cd ~/telemac/v8p2/mysimulations/steady2d-tutorial/
telemac2d.py steady2d.cas

As a result, a successful computation should end with the following lines (or similar) in Terminal:

[...]
                     *************************************
                     *    END OF MEMORY ORGANIZATION:    *
                     *************************************

 CORRECT END OF RUN

 ELAPSE TIME :
                             06  MINUTES
                             44  SECONDS
... merging separated result files

... handling result files
        moving: r2dsteady.slf
... deleting working dir

My work is done

Thus, Telemac2d produced the file r2dsteady.slf that can now be analyzed in the post-processing with QGIS or ParaView.

Post-processing

The post-processing of the steady 2d scenario uses QGIS and the PostTelemac plugin. Alternatively, TELEMAC results can also be visualized with ParaView or BlueKenue.

Load Results

Launch QGIS, create a new QGIS project, set the project CRS to UTM zone 33N, add a satellite imagery basemap, and save the project (e.g., as tm2d-postpro.qgis) in the same folder where the Telemac2d simulation results file (r2dsteady.slf is located), similar to the descriptions in the pre-processing tutorial. Then, open the PostTelemac plugin as indicated in Fig. 79.

qgis load open PostTelemac plugin

Fig. 79 Open the PostTelemac plugin in QGIS.

The PostTelemac plugin typically opens as a frame at the bottom-right of the QGIS window (maybe hard to find the first time). Detach the PostTelemac plugin from the main QGIS window by clicking on the resize window button in the top-right corner of the PostTelemac plugin frame (next to the close cross). In the detached window load the model results as follows (also indicated in Fig. 80):

  • Click on the File … button, navigate to the location where the simulation lives and select r2dsteady.slf.

  • Move the Time slider to the last timestep (e.g., 8000) and observe the main window, which will show by default the VELOCITY U parameter in this tutorial (depends on the variables defined with the VARIABLES FOR GRAPHIC PRINTOUTS keyword).

  • Familiarize with the PostTelemac plugin by modifying the display Parameter and the Color gradient.

qgis load simulation results slf PostTelemac plugin

Fig. 80 Load the Telemac2d simulation results file in the detached PostTelemac plugin window.

Find the r2dsteady layer in the QGIS Layers panel

Once imported, the r2dsteady layer is listed in the Layers panel of QGIS (typically in the bottom-left of the window). Double-clicking on the r2dsteady layer will re-open the PostTelemac plugin when it was closed (e.g., after restarting QGIS). Structurally, the r2dsteady layer is a mesh with a particular format.

SLF Geometries can be loaded in QGIS without the PostTelemac plugin

The latest QGIS versions enable direct loading of SLF geometry files as mesh layer. Try it: Drag and drop r2dsteady.slf from the Browser panel to the Layers panel. Moreover, the unsteady post-processing section features this direct mesh loading in QGIS.

Analyze Results

The PostTelemac plugin enables to export simulation results in the form of multiple formats including GeoTIFF rasters of simulation output variables or data along nodes. In addition, the evolution of a parameter over the simulation time can be exported to a video. To export or animate results, go to the Tools tab (light blue box in Fig. 80) and follow the descriptions in the following paragraphs.

Export GeoTIFF

This example features the export of a flow velocity raster at the simulation end time (8000). For this purpose, click on the RasterCreation entry of the Export menu in the Tools tab. To export a flow velocity GeoTIFF raster:

  • Set the time step to 8000 (use the field indicated in Fig. 81).

  • Select 6 : VITESSE for Parameter.

    • Vitesse is French for velocity and it is calculated as \(VITESSE = \sqrt{(VELOCITY\ U)^2 + (VELOCITY\ V)^2}\)

    • Note that VELOCITY U and VELOCITY V are the flow velocities in \(x\) and \(y\) directions, respectively.

  • In the Group frame set:

    • Cell size to 1, and

    • Extent to Full Extent.

  • Start the export by clicking on Create raster.

The processing frame can be found at the bottom of the window (scroll down by clicking on the dotted circle indicated in Fig. 81) and informs about the progress.

qgis export simulation results slf PostTelemac raster geotiff tif

Fig. 81 Export a flow velocity raster of simulation results with the PostTelemac plugin.

The successful raster creation results in a new layer called r2dsteady_raster_VITESSE, which is automatically saved as a GeoTIFF raster in the same folder where the QGIS project (*.qgz) and the r2dsteady.slf files are located.

Figure 82 shows the exported flow velocity raster in QGIS with a Magma color map (select in the layer symbology).

qgis flow velocity vitesse results slf PostTelemac raster geotiff tif

Fig. 82 The exported flow velocity (VITESSE) GeoTIFF raster in QGIS (background map: Google [Goond] satellite imagery).

Export to shapefile or mesh

The PostTelemac plugin also enables exporting to other geodata types such as vector shapefiles or meshes.

How reasonable are the results?

The flow velocity raster in Fig. 82 shows some non-zero pixels on the floodplains, beyond the riverbanks. However, because the simulated discharge of 35 m\(^3\)/s corresponds to low baseflow conditions, there should not be any water on the floodplains. These apparently wrongly modeled pixels are an artifact of the wet initial conditions that put a 1-m deep water layer all over the model. In local swales (i.e., hollows or small terrain depressions) beside the main channel, the water cannot run off and remains here until the end of the simulation. To avoid the unrealistic disconnected swales, use INITIAL DEPTH : 0 or INITIAL CONDITIONS : ZERO DEPTH, which corresponds to dry initial conditions (see the initial conditions section below). Note that initializing the model with zero depths (dry) will require open (liquid) boundary conditions of the type PRESCRIBED Q AND H. Thus, for starting this tutorial with dry conditions, go back to the pre-processing section on creating Conlim Boundary Conditions and assign prescribed Q and H to both upstream and downstream open boundaries (i.e., not only at the downstream open boundary). Alternatively, set the initial depth to a very small value such as INITIAL DEPTH : 0.01.

Verify Discharge Convergence

The discharge convergence can be observed during the simulation in the Terminal running Telemac2d with the PRINTING CUMULATED FLOWRATES : YES keyword, and also later to make sure fluxes are correctly entering and leaving the model. Furthermore, it can be determined at which timestep inflows and outflows converge to each other, which enables the simulation to be shortened to this moment of convergence in steady simulations.

To export flowrates along any line or at any node of the mesh, make sure that Q is in the list of the VARIABLES FOR GRAPHIC PRINTOUTS keyword. Then, go to the Tools tab of the PostTelemac plugin in QGIS (Fig. 83) and:

  • Click on the Flow ribbon.

  • In the Selection frame select Temporary polyline and move the mouse cursor on the map viewport where the cursor should turn into a black cross that enables drawing a (green) thick line anywhere in the mesh layer (r2dsteady). If the cursor does not enable drawing, go somewhere else in the PostTelemac plugin (e.g., to the Samplingtool ribbon), then go back to the Flow ribbon, click in the Selection frame, and re-try. To draw a line for exporting associated flows:

    • left-click with the mouse cursor somewhere on the r2dsteady mesh on the map (e.g., the left bank at the inflow open boundary indicated in Fig. 83), and

    • double left-click on another point on the r2dsteady mesh (e.g., the right bank at the inflow open boundary indicated in Fig. 83).

    • The PostTelemac plugin then automatically draws the shortest path between the two points along the mesh nodes.

  • The flowrate across the green line is now plotted in the graph of the PostTelemac plugin for the simulation time (e.g., timesteps 0 to 8000 in Fig. 83).

  • To save the values for comparison at another line, click on Copy to clipboard and paste the values into a spreadsheet (office software, such as Libre Office).

Repeat the procedure at the downstream open boundary and paste the values in another column of the spreadsheet used for the upstream open boundary.

qgis flow rate discharge control section Post Telemac convergence

Fig. 83 Draw polylines along mesh nodes and export associated flows (Copy to clipboard; background map: Google [Goond] satellite imagery).

The diagram in Fig. 84 plots the two columns of flows at the upstream and downstream open boundaries over time for the simulation setup in this tutorial. The diagram suggests that the model reaches stability after the 55th output listing (simulation time \(t \leq 5500\)). Thus, the simulation time could be limited to \(t = 6000\), but a simulation time of \(t = 5000\) would be too short.

telemac2d convergence steady model simulation discharge verification validation

Fig. 84 Convergence of inflow (upstream) and outflow (downstream) at the open model boundaries.

Control Sections are more consistent for exporting fluxes along lines

In practice, the use of CONTROL SECTIONS (unsteady tutorial) is more consistent to verify flows.

Note the difference between the convergence duration in this steady simulation with Telemac2d that starts with an initial condition of 1.0 m water depth (plot in Fig. 84) compared with the longer convergence duration in the BASEMENT tutorial (plot in Fig. 59) that starts with a dry model. This difference mainly stems from the type of initial conditions (initial depth versus dry channel) that also reflects in an outflow surplus of the Telemac2d simulation and a zero-outflow in the BASEMENT simulation at the beginning of the simulations. However, the faster convergence is at the cost of unrealistically wetted hollows in the Telemac2d simulation. Read more in the above comment: How reasonable are the results?

Re-Initialize Dry

For comparison, run the Telemac2d simulation with initial dry conditions. To this end, change the downstream boundary type to 5 4 4 (prescribed H only) in the boundaries.cli file. For making this modification, it is sufficient to open boundaries.cli in any text editor and use its find-and-replace function (e.g., CTRL + H keys in NotepadPlusPlus (Text Editor), or CTRL + F keys in Atom):

  • In the Find field type 5 5 5.

  • In the Replace with field type 5 4 4.

  • Click on Replace until all downstream boundary node types are changed.

In addition, the upstream boundary requires a prescription of flowrate and water depth to avoid supercritical fluxes at the beginning of the simulation. To this end, find-and-replace the upstream boundary settings in boundaries.cli:

  • In the Find field type 4 5 5.

  • In the Replace with field type 5 5 5.

  • Click on Replace until all upstream boundary node types are modified.

  • Save and close boundaries.cli.

In the Telemac2d steering (*.cas) file comment out the INITIAL DEPTH keyword, change the INITIAL CONDITIONS keyword to ZERO DEPTH, and change the PRESCRIBED ELEVATIONS keyword to 374.80565;371.33. In addition, the PRESCRIBED FLOWRATES needs to be changed to prescribe the discharge only at the inflow boundary with 35.;0..

5-5-5 boundaries everywhere

5 5 5 (prescribed Q and H) boundaries can alternatively applied in the here presented case, also to the downstream boundary, which would require to keep PRESCRIBED FLOWRATES : 35.;35. in the steering file. However, this setting may lead to an overdetermination of boundary conditions which may cause non-sense results (e.g., standing waves) in some model setups (in particular, simple geometries).

Moreover, to use the results of the dry initialization for morphodynamic calculations, the bottom elevation (B) must be added to the list of graphical printout variables.

In summary, the steering file should involve the following changes for a dry model initialization:

/ ... header
VARIABLES FOR GRAPHIC PRINTOUTS : U,V,H,S,Q,F,B
/ ...
PRESCRIBED FLOWRATES  : 35.;0.
PRESCRIBED ELEVATIONS : 374.80565;371.33
/ ...
INITIAL CONDITIONS : 'ZERO DEPTH'
/ INITIAL DEPTH : 1
/ ... footer

Alternatively, download the modified files:

Re-run Telemac2d and open the resulting r2d...slf file in QGIS. The following video features the flow velocity vector evolution during the dry-initialized model run:

Sebastian Schwindt @ Hydro-Morphodynamics channel on YouTube.

Compare the results of the dry and wet model initializations with regard to the following questions:

  • How does the mass balance evolve during the simulation?

  • How long did the simulations take to converge and did you need to modify the NUMBER OF TIME STEPS?

  • How reasonable are the results (double-check with the above comment)?

  • Which of the two initial conditions would you use in practice to show that the simulation is correct?

2d Calibration Parameters

The following parameters may be used for calibrating a 2d model to measurements (e.g., water surface elevation, water depth, or flow velocity data):

Avoid accuracy-reducing keyword settings

Keyword settings such as MASS-LUMPING ... : ... lead to increased smoothing (i.e., reduced accuracy) of results to increase computation speed. However, in most cases, it is worth accepting longer computation times and yielding higher accuracy, which will reduce efforts for model calibration, and thus, saves more time in the end.

What next?

A steady discharge almost never occurs in reality and can be used at maximum to calibrate the model based on (field) measurements. Once the model is well-calibrated for 2-3 steady discharges, the steady model results may be used for initializing an unsteady simulation, possibly with sediment transport.