Rainfall--runoff modelling with TELEMAC-2D
In telemac: R Interface to the TELEMAC Model Suite

knitr::opts_chunk$set(collapse = TRUE)
Sys.setlocale("LC_ALL","en_GB")

This vignette outlines how to conduct rainfall--runoff simulations with TELEMAC-2D using the telemac R package as interface. This demonstration reproduces the pluie example that comes with the TELEMAC code. The case study along with the rainfall--runoff routine and how to use TELEMAC-2D with dynamic rainfall is described by Ligier (2016).

TELEMAC-2D initialisation

Make sure the following packages are installed and load them. In addition, we define the path to a directory where the TELEMAC-2D runs shall be conducted (default is a temporary directory).

library(tidyverse)
library(sf)
library(telemac)

path_proj <- tempdir()

The study area will be simple square of 100 x 100 m^2 and the model mesh consists of triangles with a maximum area of 4 m^2.

tin_obj <- tin(list(boundary = data.frame(x = c(0, 100, 100,    0),
                                          y = c(0,   0, 100,  100))),
               s = 2, a = 4, q = 30)
plot(tin_obj, pch = ".")

The runoff module in TELEMAC-2D is based on a simple Curve Number (CN) approach, where the CN values define the infiltrability of soil as described in more detail by Ligier (2016). The geometry file may not only contain mesh elevations but also further (so-called private) variables, such as the CN values. As given in the pluie example we define a constant bathymetry with zero elevation and different CN numbers.

geo_obj <- geo(tin_obj,
               dem = list(elevation = data.frame(x = tin_obj$points[,1],
                                                 y = tin_obj$points[,2],
                                                 z = rep(0, nrow(tin_obj$points))),
                          cn = list(unit = "-",
                                    values = data.frame(
                                      x = tin_obj$points[,1],
                                      y = tin_obj$points[,2]) %>%
                                      mutate(z = case_when(
                                        between(x, 0, 49.9) & between(y, 0, 49.9) ~ 80,
                                        between(x, 0, 49.9) & between(y, 50, 100) ~ 85,
                                        between(x, 50, 100) & between(y, 0, 49.9) ~ 90,
                                        between(x, 50, 100) & between(y, 50, 100) ~ 95,
                                      ))
                                    )),
               fname = "geo.slf", title = "Constant bathymetry (zero elevation) with CN values")
plot(geo_obj, s = 1)
plot(geo_obj, s = 1, v = "cn")

With the boundary conditions inferred from the geometry object we can conclude the basic setup. Here we define a closed boundary, i.e. runoff will accumulate and may not leave the study area.

cli_obj <- cli(geo_obj, "bnd.cli")

Example I: constant rainfall without runoff

The first example will show how rainfall can be added to a TELEMAC-2D simulation without infiltration. This can be done by setting the steering parameter RAIN OR EVAPORATION = YES and specifying RAIN OR EVAPORATION IN MM PER DAY (here 100 mm) and DURATION OF RAIN OR EVAPORATION IN HOURS (here 6 hours resulting in a total rainfall of 25 mm).

cas_pars <- list(
  # files
  "BOUNDARY CONDITIONS FILE" = "bnd.cli",
  "GEOMETRY FILE" = "geo.slf",
  "RESULTS FILE" = "results.slf",
  # General
  # H: water depth (m);
  # U, V: velocity along in x, y direction (m/s);
  # S, B: free surface, bottom elevation (m)
  # N, O, R, Z: additional 'private' variables defined in user code
  "VARIABLES FOR GRAPHIC PRINTOUTS" = "H",
  "DURATION" = 3600*7, # seconds
  "TIME STEP" = 60, # seconds
  "GRAPHIC PRINTOUT PERIOD " = 60, # number of time steps
  "LISTING PRINTOUT PERIOD" = 60, # number of time steps
  # rainfall
  "RAIN OR EVAPORATION" = "YES",
  "RAIN OR EVAPORATION IN MM PER DAY" = 100,
  "DURATION OF RAIN OR EVAPORATION IN HOURS" = 6,
  # Numerics
  # 1: conjugate gradient, recommended when using wave equation,
  # time step should be adapted until convergence is reached after 20 to 30 iterations
  "SOLVER" = 1,
  "MAXIMUM NUMBER OF ITERATIONS FOR SOLVER" = 200,
  "MASS-BALANCE" = "YES",
  "TIDAL FLATS" = "YES",
  # 1: surface gradient corrected (recommended);
  # 2: areas masked from computation;
  # 3: like 1 with porosity term added to half-dry elements to change water quantity
  "OPTION FOR THE TREATMENT OF TIDAL FLATS" = 1,
  # 0: no treatment; 1: smoothing negative depths;
  # 2: flux limitation by segment ensuring positive depths;
  # 3: flux limitation by element
  "TREATMENT OF NEGATIVE DEPTHS" = 2,
  # required if TREATMENT OF NEGATIVE DEPTHS = 2 (not documented, see forum)
  "MASS-LUMPING ON H" = 1,
  # required if TREATMENT OF NEGATIVE DEPTHS = 2
  "CONTINUITY CORRECTION" = "YES",
  # SUPG OPTION=...;0;... required if TREATMENT OF NEGATIVE DEPTHS = 2
  "SUPG OPTION" = "0;0",
  # 1: coupled; 2: wave equation (more stable)
  "TREATMENT OF THE LINEAR SYSTEM" = 2,
  # recommended for steep topography and tidal flats
  "FREE SURFACE GRADIENT COMPATIBILITY" = 0.9,
  # Initial conditions
  "INITIAL CONDITIONS" = 'ZERO DEPTH',
  # Friction
  # 0: no friction, 2: Chezy, 3: Strickler, 4: Manning
  "LAW OF BOTTOM FRICTION" = 3,
  # concrete: 100; straight stream: 30-40; natural stream with wood: <10
  "FRICTION COEFFICIENT" = 50
)
cas_obj <- cas(cas_pars, fname = "rainfall.cas")

Now we can compile the full TELEMAC-2D project.

t2d_obj <- t2d("Example: static rainfall without runoff",
               wdir = paste(path_proj, "rainfall", sep = "/"),
               cas = cas_obj, geo = geo_obj, cli = cli_obj)

Finally we can export the TELEMAC-2D input files and run the model.

write_t2d(t2d_obj)
t2d_sim <- simulate_t2d(t2d_obj, exec = "telemac2d.py")

t2d_sim <- t2d_obj
t2d_sim$res <- results(system.file("telemac/rainfall_runoff/res_rainfall.slf", package = "telemac"),
                       times = 7*3600)

When plotting the water depth we see that at the end of the simulation period the area is covered by water at a depth of 25 mm.

res_df <- tin2grid(t2d_sim$res, s = 2, output = "data.frame")
ggplot(res_df, aes(x = x, y = y, fill = value)) +
  geom_raster() +
  coord_equal() +
  scale_y_continuous(expand = expansion(0,0)) +
  scale_x_continuous(expand = expansion(0,0)) +
  theme_bw()

Example II: constant rainfall with runoff and CN values read from a file

For this example we will use the runoff module of TELEMAC-2D. First we define CN values in an external file. We only need to provide the node points that are then used for interpolation to grid within TELEMAC's bief module. This procedure is for illustration, the CN values in the geometry file will still be ignored so far.

cn_dat <- rbind(
  matrix(c(   0, 49.9, 49.9,    0,    0,
              0,    0, 49.9, 49.9,    0,
            rep(80, 5)), ncol = 3),
  matrix(c(   0, 49.9, 49.9,    0,    0,
           50.1, 50.1,  100,  100, 50.1,
           rep(85, 5)), ncol = 3),
  matrix(c(50.1,  100,  100, 50.1, 50.1,
              0,    0, 49.9, 49.9,    0,
           rep(90, 5)), ncol = 3),
  matrix(c(50.1,  100,  100, 50.1, 50.1,
           50.1, 50.1,  100,  100, 50.1,
           rep(95, 5)), ncol = 3)
)
dir.create(paste(path_proj, "rainfall_runoff", sep = "/"), recursive = T)
write("#", paste(path_proj, "rainfall_runoff/cn.txt", sep = "/"), sep = "\n")
write("# x y cn", paste(path_proj, "rainfall_runoff/cn.txt", sep = "/"), append = T, sep = "\n")
write.table(cn_dat, paste(path_proj, "rainfall_runoff/cn.txt", sep = "/"),
            row.names = F, col.names = F, quote = F, sep = "\t", append = T)

The actual runoff module is included as user-defined Fortran code that has to be specified as steering parameter FORTRAN FILE. The code also defines some additional variables (accumulated rainfall and accumulated runoff), also referred to as private variables. The external input file providing the CN values has to be given as FORMATTED DATA FILE 2, also defined in the user code. The code was taken from the pluie example in TELEMAC-2D.

Note: Lines in the steering file must be shorter than 72 characters, whereas longer values can be split over multiple lines. However, when having a long file path it cannot be split over multiple lines as the slashes will then be interpreted as comments. Therefore you might have to copy the Fortran code from system.file("telemac/rainfall_runoff/code_rainstat_cnfile", package = "telemac") to a different location resulting in a shorter file path. You may also provide a path relative to the project directory.

# path might be too long, see note above
cas_obj[["FORTRAN FILE"]] <- system.file("telemac/rainfall_runoff/code_rainstat_cnfile",
                                         package = "telemac")
cas_obj[["FORMATTED DATA FILE 2"]] <- "cn.txt"
# 'acc. runoff' as defined in user code
cas_obj[["VARIABLES FOR GRAPHIC PRINTOUTS"]] <- "R"
# 1: on; 0: off
cas_obj[["RAINFALL-RUNOFF MODEL"]] <- 1
# 1: 0.2 (standard); 2: revised method, 0.05 with conversion of CN values
cas_obj[["OPTION FOR INITIAL ABSTRACTION RATIO"]] <- 1
# 1: dry, 2: normal, 3: wet
cas_obj[["ANTECEDENT MOISTURE CONDITIONS"]] <- 2

Now we can initialise a new t2d setup and run the model.

t2d_obj <- t2d("Example: static rainfall with runoff, CN from file",
               wdir = paste(path_proj, "rainfall_runoff", sep = "/"),
               cas = cas_obj, geo = geo_obj, cli = cli_obj)

write_t2d(t2d_obj)
t2d_sim <- simulate_t2d(t2d_obj, exec = "telemac2d.py")

t2d_sim <- t2d_obj
t2d_sim$res <- results(system.file("telemac/rainfall_runoff/res_rainfall_runoff.slf", package = "telemac"),
                       times = 3600 * c(0,1,4,7))

We will plot the user-defined variable acc. runoff, i.e. total accumulated runoff at each timestep. It can be seen how accumulated runoff depends on the CN values. All in all, the values are as reported by Ligier (2016).

res_df <- tin2grid(t2d_sim$res, s = 2, output = "data.frame")
ggplot(res_df, aes(x = x, y = y, fill = value)) +
  geom_raster() +
  coord_equal() +
  scale_y_continuous(expand = expansion(0,0)) +
  scale_x_continuous(expand = expansion(0,0)) +
  facet_wrap(~ timestep) +
  theme_bw()

Example III: dynamic rainfall with runoff and CN values from geometry file

In the last example it will be shown how to provide dynamic rainfall intensities and directly use the CN values we earlier stored in the geometry file. First we define a rainfall event with the same accumulated depth as in the last examples (25 mm) but with varying intensity.

dir.create(paste(path_proj, "rainfall_runoff_dyn", sep = "/"), recursive = T)
write("#", paste(path_proj, "rainfall_runoff_dyn/rain.txt", sep = "/"), sep = "\n")
write("# time (s) rainfall (mm)", paste(path_proj, "rainfall_runoff_dyn/rain.txt", sep = "/"), append = T, sep = "\n")
rain <- data.frame(t = c(0, 3600, 3*3600, 4*3600, 6*3600, 8*3600),
                   r = c(0,   10,      5,      7,      3,      0))
write.table(rain, paste(path_proj, "rainfall_runoff_dyn/rain.txt", sep = "/"),
            row.names = F, col.names = F, quote = F, sep = "\t", append = T)

The rainfall input is given as steering parameter FORMATTED DATA FILE 1 (defined in the user code). To use CN (or any other private variables) from the geo file the steering parameter NAMES OF PRIVATE VARIABLES has to be provided. Also the user code has to be adapted. In the example code parameter RAINDEF has to be set to a value of 3 (instead of 1) and lines 248 to 274 need to be adapted (see code in system.file("telemac/rainfall_runoff/code_raindyn_cngeo", package = "telemac") and compare with the user code of example II).

# path might be too long, see note above
cas_obj[["FORTRAN FILE"]] <- system.file("telemac/rainfall_runoff/code_raindyn_cngeo",
                                         package = "telemac")
cas_obj[["FORMATTED DATA FILE 1"]] <- "rain.txt"
cas_obj[["FORMATTED DATA FILE 2"]] <- NULL # not needed here
# as given in geometry file and used in user code
cas_obj[["NAMES OF PRIVATE VARIABLES"]] <- "CN"

Now we can initialise a new t2d setup, run the model, and plot again the accumulated runoff. Note how accumulated runoff over different time steps differs from example II, whereas at the end of the simulation total accumulated runoff is the same as total rainfall depths are equal.

t2d_obj <- t2d("Example: dynamic rainfall with runoff, CN from geo",
               wdir = paste(path_proj, "rainfall_runoff_dyn", sep = "/"),
               cas = cas_obj, geo = geo_obj, cli = cli_obj)

write_t2d(t2d_obj)
t2d_sim <- simulate_t2d(t2d_obj, exec = "telemac2d.py")

t2d_sim <- t2d_obj
t2d_sim$res <- results(system.file("telemac/rainfall_runoff/res_rainfall_runoff_dyn.slf", package = "telemac"),
                       times = 3600 * c(0,1,4,7))

res_df <- tin2grid(t2d_sim$res, s = 2, output = "data.frame")
ggplot(res_df, aes(x = x, y = y, fill = value)) +
  geom_raster() +
  coord_equal() +
  scale_y_continuous(expand = expansion(0,0)) +
  scale_x_continuous(expand = expansion(0,0)) +
  facet_wrap(~ timestep) +
  theme_bw()