#  Basement v3: a modular freeware for river process modelling over multiple
computational backends

Davide Vanzo  Samuel Peter  Lukas Vonwiller  Matthias Buergler  Manuel
Weberndorfer  Annunziato Siviglia  Daniel Conde  David F. Vetsch  Laboratory
of Hydraulics, Hydrology and Glaciology. ETH, Swiss Federal Institute of
Technology, Zürich, Switzerland.  Dept. of Surface Waters - Research and
Management. Eawag, Swiss Federal Institute of Aquatic Science and Technology,
Kastanienbaum, Switzerland.  Dept. Civil, Environmental and Mechanical
Engineering. University of Trento, Trento, Italy.  ID Scientific IT Services.
ETH, Swiss Federal Institute of Technology, Zürich, Switzerland.  Axpo Group,
Baden, Switzerland.  TK CONSULT AG , Zürich, Switzerland.

###### Abstract

Modelling river physical processes is of critical importance for flood
protection, river management and restoration of riverine environments.
Developments in algorithms and computational power have led to a wider spread
of river simulation tools. However, the use of two-dimensional models can
still be hindered by complexity in the setup and the high computational costs.
Here we present the freeware Basement version 3, a flexible tool for two-
dimensional river simulations that bundles solvers for hydrodynamic,
morphodynamic and scalar advection-diffusion processes. Basement leverages
different computational platforms (multi-core CPUs and graphics processing
units GPUs) to enable the simulation of large domains and long-term river
processes. The adoption of a fully costless worflow and a light GUI facilitate
its broad utilization. We test its robustness and efficiency in a selection of
benchmarks. Results confirm that Basement could be an efficient and versatile
tool for research, engineering practice and education in river modelling.

###### keywords:

GPU-CUDA , river modelling , unstructured grid , shallow water , sediment
transport , pollutant transport

††journal: Environmental Modelling & Software

## Software availability

Name of software:

Basement, version 3 (v3)

Website:

www.basement.ethz.ch

E-mail:

<EMAIL_ADDRESS>
Developer:

Numerical Modelling division at the Laboratory of Hydraulics, Hydrology and
Glaziology (VAW), ETH Zürich

Language:

C++, C, CUDA

Interface:

graphical user interface (GUI), command-line interface (CLI)

Hardware:

CPUs, CUDA-enabled GPUs (optional)

OS:

Windows, Linux (Ubuntu)

Availability:

Freeware

##  1 Introduction

In the last decades the usage of numerical tools widely spread in several
subjects of the environmental sciences. River science [sensu 1] is no
exception in this trend, with a number of tools been developed to address
variegate research questions [e.g. 2, 3]. Modelled river physical processes
span from flood simulation, hydraulic and sediment dynamics, pollutant and
temperature transport, to vegetation and flow interactions, just to mention a
few [e.g. 4, 5, 6, 7, 8]. Such river processes occur at different spatial and
temporal scales, hence influencing the development and choice of suitable
modelling tools.

Advances in computational power, numerical algorithms and optimization
routines that occurred in the last decades allowed for the spread of more and
more sophisticated numerical tools. In the context of river science, two-
dimensional depth-averaged (hereinafter 2D) models are nowadays of common use
in research and engineering practice. This is particularly true for some
applications such as flood modelling and river morphodynamics [e.g. 3, 9]. The
increasing usage of 2D river models is also closely bonded with the growing
availability of high-resolution river datasets. In particular, advances in
LiDAR, UAV-Photography and others remote-sensed survey technologies enable
river topographic scans at an unprecedented level of detail [e.g. 10, 11].

The increased computational capabilities and refined datasets open the gates
for near-census [sensu 12] numerical modelling of several river processes.
Indeed, 2D river models have the capability to simulate fine spatial
(centimeters to meters) and temporal scales (seconds to days). At such scales,
relevant hydro-morphodynamic processes (e.g. bar formation) and also
ecohydraulic processes (e.g. habitat dynamics) can hence be modelled [e.g. 13,
14, 15, 16]. Nevertheless, 2D river models can still be computationally
demanding, with simulations lasting several days. This is particularly true
when complex physical processes, such as morphodynamics, are accounted for
[17]. Moreover, large-scale or near-census applications (i.e. with millions of
computational cells) and/or long-term simulations (i.e. years) all concur to
increase overall computational costs. Such drawbacks particularly apply for
the investigation of highly unsteady river processes such as artificial or
natural flood waves, where explicit numerical schemes are preferable. In such
cases the overall computational time scales exponentially with the number of
computational cells due to stability constraints [e.g. 18].

Increasing the efficiency and the computational performance of river models
represents yet a challenge. Pitfalls arise with the number of computational
cells, but also with the inherent complexity of 2D models. For example,
challenges are to be found in the setup of the computational domain [19] but
also in the definition of particular boundary conditions [20, 21]. Increasing
the computational performance is also sought by developing alternative
numerical solution strategies for the underlying physical governing equations.
Examples are variegate, spanning from the adoption of a local timestep for the
numerical integration [e.g. 22, 23], the automatic adaptation of the
computational mesh [e.g. 24], the use of acceleration factors for the hydro-
morphodynamic problem [e.g 25, 26], to the reformulation of the governing
mathematical equations for water quality simulations [e.g. 27], to mention a
few.

Parallel computing solutions are the most popular strategies to reduce
computational time. They historically benefit from the continuous improvements
of computational performances of both single CPUs and clusters. The general
aim of parallelization techniques is to split the total computational load
into tasks that can be executed simultaneously by different computational
units [e.g. 28]. The use of the Graphics Processing Unit (GPU) as a general-
purpose computational resource developed rapidly in the last decade [29]. For
many parallelizable workloads, offloading work to GPUs is a relatively cheap
and efficient high-performance computing strategy that is also easily
upgradable in standard desktop workstations.

By means of GPU parallelization, numerical models can potentially be
accelerated by an order of magnitude and more. The efficiency of such a
parallelization depends on the data exchange between the main memory and the
processors, with a complex memory hierarchy and bandwidth bottlenecks [30].
These low-level constraints can limit computational speedup and depend on the
model data/memory handling, the underlying model complexity and the nature of
the governing equations to be solved. In applications such as 2D river models
the type of computational mesh, i.e. structured or unstructured, has a
significant influence on the final computational speedup.

In the last decade, river simulation models have benefited from GPU
parallelization. Specific and ad-hoc implementations of GPU-based models for
2D hydrodynamic [e.g. 31, 32, 33, 34], and occasionally morphodynamic
simulations [e.g. 35] have become available. The vast majority of these models
are based on structured grids which allow for an easier implementation and for
relatively higher computational speedups. This is due to the fact that, for
structured meshes, the data structure is inherently simpler, which reduces the
need for mappings and indirections. To the best of the Authors’ knowledge, few
hydrodynamic models implement GPU-acceleration on unstructured meshes [36, 37,
38, 39], with very limited ad-hoc implementations for transient flows
morphodynamics [e.g. 40].

Bundled river modelling software that supports GPU acceleration is available
for commercial use (e.g. RiverFlow2D (hydronia.com/riverflow2d), TUFLOW
(tuflow.com), but costless ones are still few [see 41]. An increase in
availability of freeware GPU-based river models would be beneficial for
environmental modelers in academic research and education, but also in
consultancy and engineering offices.

In this paper we introduce the Basement software (version 3), a freeware
application developed at the Laboratory of Hydraulics, Hydrology and
Glaciology of ETH Zürich. The software can simulate two-dimensional
hydrodynamic, morphodynamic, and scalar advection-diffusion processes of
scientific and practical interest. It can seamlessly run on GPU-enabled
workstations, as well as on more standard multi-core CPUs. This flexibility in
the choice of the computational backend is achieved by integrating the OP2
framework [30, 42, 43]. This framework provides an additional abstract layer
for the acceleration of numerical models on unstructured computational meshes,
and has been successfully implemented in similar modelling context [44]. The
obtained parallelization performance alleviates the computational limitations
when simulating high resolution (or large) computational domains and/or long
term processes [e.g. 45]. This is particularly relevant when aiming at the
calibration [46] or at the uncertainty evaluation [47, 11] of deterministic
models. As proof of concept, a flood wave uncertainty propagation analysis
with Basement has been proposed in [48].

In the current version Basement is available for both Windows and Linux-based
(Ubuntu) environments. It is provided with a Command Line Interface (CLI) to
easily perform batch simulations, but also with a light Graphical User
Interface (GUI). The Basement software aims to enable a broad range of
potential users to skilfully simulate river processes in the domain of river
engineering and research on state-of-the-art computational hardware. Moreover,
with accompanying scholar programs and extended documentation material and
tutorials, the software is designed to be a valuable didactic tool for
engineering and river science students.

The paper is structured as follows: §2 provides the software application
context that justifies the adopted mathematical and numerical strategies. §3
to §5 report the mathematical basis, the numerical strategies and main
features of the basic modules of Basement. The software design, the modelling
workflow and the parallelization solutions are presented in §6. A selection of
benchmarks are reported in §7, whilst conclusions and outlooks are drawn in
§8.

##  2 Application context

One of the main goals of the novel software design of version 3 is the
capability to tackle river processes at different spatial and temporal scales.
For example, Basement can be used to simulate large scale (i.e. basin scale)
flood propagation, but also reach scale morphodynamic processes such as
formation and evolution of fluvial bars. Moreover, it can be applied together
with high-resolution topographies (in the order of centimeters) to simulate
ecohydraulic processes at different ecological scales (e.g. habitat
modelling). This range of application possibilities is enabled by specific
characteristics of the software. In particular:

unsteady and transitional flows:

Basement can deal with strongly unsteady flows and different flow regimes
(sub- and super-critical). For this reason, Basement is particularly suitable
for simulations of river flows in Alpine contexts, the propagation of natural
flood waves as well as hydropeaking events. This is ensured by the adoption of
a robust and accurate shock-capturing explicit solver for the hydrodynamic
problem (§4.2);

accurate front propagation:

it is possible to simulate extreme events such as dam-break induced floods,
but also ecologically-relevant processes such as the wetting-drying of
riparian areas and in-channel morphologies due to artificial flow alterations.
This is achieved by an implemented shock-wave capturing numerical scheme
complemented with a robust treatment of wet-dry interfaces (§4.2);

complex river topographies:

the use of unstructured grid for the computational domain discretization
enables for an accurate description of complex river morphologies and
riverrine structures (§4.1). The adoption of an unstructured mesh also reduces
the strong anisotropy of structured meshes, which can be crucial for
particular applications.

large problems:

the software adopts a parallelization strategy tailored to the acceleration of
problems on unstructured meshes (§6.4). Moreover, Basement simulations can
efficiently be executed on different computational backends. Those backends
inlude GPU cards, therefore allowing for the simulation of large domains
(millions of computational cells) on standard workstations, having a limited
cost.

multiple river processes:

the software is designed in a modular way, so different river processes such
as hydrodynamics, sediment or advection-diffusion of a scalar (e.g. a non-
reactive pollutant) can be simulated by activating specific modules at setup
time (§6). Different types of boundary conditions (§5) and closure
relationships are available to simulated, for example, simple hydraulic
structures (e.g. weirs) or flow inputs/outputs (e.g. water intakes). The
modular design (§6) allows to retain good parallelization performances in the
simulation of different river processes, as shown in (§7.7).

The basic modules available in Basement are i) hydrodynamics, ii)
morphodynamics and iii) advection-diffusion of scalar quantities. Each module
is composed by different sets of hyperbolic equations describing the
conservation and evolution of the water flow (hydrodynamics), the fluvial
sediment (morphodynamics) and the concentration of passive solutes (scalar
advection-diffusion). The governing equations represent a so-called Initial-
Boundary Value Problem [18], where process-specific initial and boundary
conditions are required to be set. The following Sections presents the main
governing equations and closure relationships (§3), the numerical strategies
(§4) and finally the initial and boundary conditions (§5) for the three basic
modules. The main module features are also listed in Table 1 of the
Supplementary Material.

##  3 Mathematical formulation

###  3.1 Hydrodynamics

The hydrodynamic module solves the so-called shallow water equations
(hereinafter SWE) [e.g. 18]. The two-dimensional SWE are of practical interest
with regard to water flows with a free surface under the influence of gravity.

Considering a Cartesian reference system $(x,y,z)$ where the $z$ axis is
vertical and the $(x,y)$ plane is horizontal (Fig. 1a), the system of
governing equations can be expressed as:

$\begin{cases}\partial_{t}H+\partial_{x}q_{x}+\partial_{y}q_{y}=S_{h}\\\
\partial_{t}q_{x}+\partial_{x}\left(\dfrac{q_{x}^{2}}{h}+\dfrac{1}{2}gH^{2}-gHz_{B}\right)+\partial_{y}\left(\dfrac{q_{x}q_{y}}{h}\right)=-gH\partial_{x}z_{B}-ghS_{fx}\\\
\partial_{t}q_{y}+\partial_{x}\left(\dfrac{q_{x}q_{y}}{h}\right)+\partial_{y}\left(\dfrac{q_{y}^{2}}{h}+\dfrac{1}{2}gH^{2}-gHz_{B}\right)=-gH\partial_{y}z_{B}-ghS_{fy}\\\
\end{cases}$ (1)

where the system unknowns are the water surface elevation $H$ [m], and the two
directional components of ${\mathbf{q}}=(q_{x},q_{y})$
[$\mathrm{m}^{2}\mathrm{/}\mathrm{s}$], representing the flow discharge per
unit width. With $z_{B}$ [m] we indicate the bottom elevation, whilst
$h=(H-z_{B})$ [m] is the water depth, and $g$
[$\mathrm{m}\mathrm{/}\mathrm{s}^{2}$] the acceleration due to gravity. Note
that the depth-averaged velocity vector can be consequently expressed as
${\mathbf{u}}=(u,v)=(q_{x}/h,q_{y}/h)$ [$\mathrm{m}\mathrm{/}\mathrm{s}$].
Finally $S_{fx}$ and $S_{fy}$ [-] represent the dimensionless friction terms
in $x$ and $y$ direction, whilst $S_{h}$ [$\mathrm{m}\mathrm{/}\mathrm{s}$]
represents potential external contribution/subtraction of flow discharge to
the mass conservation equation.

$z$
---
$x$
---
$y$
---
$H$
---
$h$
---
$z_{B}$
---
$\varphi_{c}$
---
$\varphi_{b}$
---
$\mathbf{q}$
---
$\mathbf{q}_{B}$
---
$s$
---
$n_{q}$
---
$R_{c}$
---
(a)
---
(b)
---
(c)
---
$z_{rel}$
---
Figure 1: Notation of scalar and vectorial quantities. (a) Reference system
$(x,y,z)$ with water surface elevation $H$, water depth $h$, bed elevation
$z_{B}$ and non-erodible fixed bed depth $z_{rel}$. (b) Bed load transport
($\mathbf{q_{B}}$) deviation angle $\varphi_{b}$ from the flow direction
$\mathbf{q}$ due to gravitational effects caused by the local lateral slope
$s$. (c) Bed load transport ($\mathbf{q_{B}}$) deviation angle $\varphi_{c}$
from the flow direction $\mathbf{q}$ due to the spiral flow motion caused by
the curvature of radius $Rc$.

####  3.1.1 Hydrodynamic closure relationships

To solve the system (1), closure relationships for the friction terms
$S_{fx},S_{fy}$ and the contribution of external inflow/outflow discharge
$S_{h}$ must be provided.

##### Friction terms

Under the hypothesis of turbulent flow, hence under the assumption that the
energy line slope is proportional to the square of the flow velocity, the
friction terms $S_{fx},S_{fy}$ can be written as:

$S_{fx}=\frac{u\lVert{\mathbf{u}}\rVert}{ghc_{f}^{2}}\;;\;S_{fy}=\frac{v\lVert{\mathbf{u}}\rVert}{ghc_{f}^{2}},$
(2)

where $c_{f}$ is the dimensionless friction coefficient. Several formulae are
available for $c_{f}$. Basement implements four well know formulations of
power or logarithmic type, given in Table 1.

Table 1: Friction closure relationships for the hydrodynamic problem.
Formulations for the dimensionless friction coefficient $c_{f}$; for both
Chézy and Bezzola entries, $d_{90}$ is the 90th percentile of the sediment
grain size distribution.

Closure Expression Parameters Range Ref Strickler
$c_{f}=k_{str}h^{1/6}/\sqrt{g}$ $k_{str}$ [$\mathrm{m}^{1/3}\mathrm{s}^{-1}$]
7-40 [49] Manning $c_{f}=h^{1/6}/(n\sqrt{g})$ $n$
[$\mathrm{m}^{-1/3}\mathrm{s}$] 0.025-0.143 [49] Chézy
$\left.\begin{array}[]{lccl}c_{f}=6.25+\log{\left(h/K_{s}\right)}&\quad\text{for}&&h>K_{s}\\\
c_{f}=6.25&\quad\text{for}&&h\leq K_{s}\end{array}\right.$ $K_{s}$ [m]
$K_{s}=n_{k}d_{90}$, [50] with $n_{k}=2\div 3$ Bezzola
$\left.\begin{array}[]{lll}c_{f}=2.5\sqrt{1-y_{R}/h}\;\ln{\left(10.9h/y_{R}\right)},&\quad\text{for}&\quad
h/y_{R}>2\\\
c_{f}=1.25\sqrt{h/y_{R}}\;\ln{\left(10.9h/y_{R}\right)},&\quad\text{for}&\quad
0.5\leq h/y_{R}\leq 2\\\ c_{f}=1.5,&\quad\text{for}&\quad
h/y_{R}<0.5\end{array}\right.$ $y_{R}$ [m] $y_{R}=nd_{90}$, $n\approx 1$ [51]

##### External inflow/outflow discharge

The term $S_{h}$ [$\mathrm{m}\mathrm{/}\mathrm{s}$] is used to represent
additional sources of water like rainfall and springs, or water abstraction
(sink), and can be defined over subsets of the computational domain. The
external water source can be provided by the user as total discharge
[$\mathrm{m}^{3}\mathrm{/}\mathrm{s}$] or as intensity
[$\mathrm{m}\mathrm{m}\mathrm{/}\mathrm{h}$], per squared meter. Different
behaviour can be imposed for each external source/sink area:

  * 1.

Exact (source/sink): the exact given water volume is added (source) or
extracted (sink) from the surface. This is the only option for water addition.
In case of water abstraction, the simulation might end abruptly if the
available water volume is smaller than the volume prescribed for subtraction.
This option allows to have the full control on the water entering/leaving the
computational domain, and is useful to simulate e.g. managed hydraulic
structures, such as regulated water intakes.

  * 2.

Available (sink): the given water volume to extract is limited by the
available water volume in the single element (i.e. computational cell). With
this abstraction option, the simulation proceeds with no interruptions because
the water conservation is ensured. This option is useful to simulate
particular unmanaged hydraulic structures, such as diversion spillways.

  * 3.

Infinity (sink): all available water will be abstracted from the computational
domain.

###  3.2 Morphodynamics

The basic morphodynamic module solves the so-called Exner equation [52]. It
describes the bed evolution due to erosion or deposition, which results in the
elevation change of the actual bed level $z_{B}$. Assuming the same coordinate
reference of Fig. 1a, it reads

$(1-p)\partial_{t}z_{B}+\partial_{x}q_{B_{x}}+\partial_{y}q_{B_{y}}=S_{b},$
(3)

where $p$ [-] is the bed sediment porosity, assumed constant in space and
time, $S_{b}$ [$\mathrm{m}\mathrm{/}\mathrm{s}$] is an external source term
specifying local inputs or outputs of sediment material (e.g. slope collapse
or excavation) and ${\mathbf{q}_{B}}=(q_{B_{x}},q_{B_{y}})$
[$\mathrm{m}^{2}\mathrm{/}\mathrm{s}$] is the specific sediment transport
flux.

The morphodynamic module in its basic form accounts only for sediment
transport occurring in the form of bed-load or total-load [49, 53]. The
simulation of sediment transport as suspended load is delegated to a specific
module, planned in future versions of the software.

####  3.2.1 Morphodynamic closure relationships

Two closure relationships are needed to numerically solve the problem of (3):
a sediment transport formula and the external source/sink of sediments.

##### Sediment transport formula

Basement implements three different types of sediment transport formulae as
given in Table 2. The first two, Meyer-Peter and Müller like (MPM-like), and
Grass like (Grass-like) expressions are adequate to simulate bed-load
dominated sediment transport condition. The Engelund and Hansen formulation
allows for the estimation of total sediment transport (i.e. suspended and bed-
load), whilst Smart and Jäggi account for bedload transport in steep channels.

The expressions and the typical parameter values to calculate the specific
sediment transport magnitude $\lVert{\mathbf{q}_{B}}\rVert$ are given in Table
2. In the MPM-like formulation, $\theta$ is the dimensionless bed shear stress
(i.e. Shields parameter [49]), $\theta_{cr}$ is the critical dimensionless bed
shear stress, $d_{m}$ is the representative grain diameter, $s=\rho_{s}/\rho$
is the relative density of the sediment with respect to water. The
coefficients $\alpha$, $m$ and the critical threshold $\theta_{cr}$ can be
assigned by the user or adopted from literature (see Table 2). The Grass-like
model proposes a simple bedload transport formula, where
$\lVert{\mathbf{q}_{B}}\rVert$ is a function of the flow velocity magnitude,
with $u_{cr}$ as critical threshold velocity. The coefficients $\alpha$, $m$
and the critical threshold $\mathbf{u}_{cr}$ can be assigned by the user or
adopted from literature (Table 2). Engelund and Hansen [54] proposed a total
transport formula for uniform bed material without a threshold condition for
incipient motion.

Table 2: Sediment transport closure relationships for the morphodynamic
problem. Expressions provide an estimation of the specific sediment transport
magnitude $\lVert{\mathbf{q}_{B}}\rVert$
[$\mathrm{m}^{2}\mathrm{/}\mathrm{s}$].

Type Expression Parameters Ref MPM-like
$\alpha(\theta-\theta_{cr})^{m}\sqrt{(s-1)gd_{m}^{3}}$
$\alpha=8,m=1.5,\theta_{cr}=0.047$ [55] $\alpha=4.93,m=1.6,\theta_{cr}=0.047$
[56] Grass-like $\alpha(\lVert{\mathbf{u}}\rVert-u_{cr})^{m}$
$\alpha\approx\mathcal{O}(-2;-3),m=3,u_{cr}=0.0$ [57] Engelund and Hansen
$0.05c_{f}^{2}\theta^{5/2}\sqrt{(s-1)gd_{m}^{3}}$ $-$ [54] Smart and Jäggi
$\alpha\left(\frac{d_{90}}{d_{30}}\right)^{0.2}J^{0.6}\lVert{\mathbf{u}}\rVert(\theta-\theta_{cr})d_{m}$
$\alpha=8,\theta_{cr}=0.05$ [58]

##### Local corrections of the sediment transport

The morphodynamic module implements three corrections to the basic Exner
equation (3) to account for the influence of local characteristics of the flow
and the bottom on the sediment transport. Namely, i) the influence of local
slope on incipient motion, ii) the effect of lateral bed slope and iii) of the
flow curvature on the sediment transport direction.

The threshold condition for incipient motion of grains, by Shields [59], is
valid for an almost horizontal bed. In case of a sloped bed in flow direction
or transverse to it, the stability of grains is either increased or reduced
due to the gravity. The critical shear stress value can be adapted
consequently to account for the influence of local longitudinal and
transversal slopes. A common approach is to scale the critical shear stress
for almost horizontal bed $\theta_{cr}$ with a correction factor $k$:

$\theta^{*}_{cr}=k\theta_{cr}.$ (4)

Basement implements the correction factor $k$ as proposed in [60] and [61].
Implementation details are given in the official documentation.

The bedload direction can be corrected to account for two relevant
morphodynamic processes linked to the slope of the bed and the curvature of
the flow. The deviation of the bedload direction from the flow direction can
thus be modelled as a deviation angle $\varphi=\varphi_{b}+\varphi_{c}$, sum
of the correction angle for bed slope ($\varphi_{b}$) and curvature
($\varphi_{c}$), as depicted in Fig. 1b and c. The bedload vector is then
rotated with the rotation matrix ${\bf T}(\varphi)$, being

${\bf T}=\begin{bmatrix}cos\varphi&-sin\varphi\\\ sin\varphi&cos\varphi\\\
\end{bmatrix},$ (5)

where the angle is positive counterclockwise.

The angle $\varphi_{b}$ is estimated with the approach proposed in [62] and
[63] for the effect of the local transversal bed slope. In particular, the
bedload direction deviates from the flow direction in presence of a local
transversal bed slope, due to the gravity acting on the bedload sediment
particles (Fig. 1b). The bed load deviation $\varphi_{b}$ with respect to the
flow is therefore evaluated as

$\tan\varphi_{b}=-N_{l}\sqrt{\dfrac{\theta_{cr}}{\theta}}\cdot\mathbf{s}\cdot\mathbf{n}_{q},\;\text{for}\;\mathbf{s}\cdot\mathbf{n}_{q}<0,$
(6)

where $N_{l}$ is an experimental lateral transport factor ($0.75\leq N_{l}\leq
2.63$), $\mathbf{s}=\left(\partial_{x}z_{B},\partial_{y}z_{B}\right)$ is the
local bed slope and $\mathbf{n}_{q}$ is the unit vector perpendicular to
$\mathbf{q}$ in downhill direction (Fig. 1b).

The angle $\varphi_{c}$ accounts for the effect of a marked flow curvature.
Due to three dimensional spiral flow motion that establishes in curved flows,
the bed load direction tends to point towards the inner side of the curve,
while the flow direction points towards the outer side (Fig. 1c). This
curvature effect is taken into account according to an approach proposed in
[64], where the deviation angle $\varphi_{c}$ is determined as

$\tan{\varphi_{c}}=-N_{*}\dfrac{h}{R},$ (7)

where $h$ is the water depth, $N_{*}$ is a curvature factor, and $R$ denotes
the radius of the river bend, positive for curvature in counterclockwise
direction. The curvature factor $N_{*}$ mainly depends on bed roughness and
assumes values $N_{*}\approx 7$ for natural streams [64], and values up
$N_{*}\approx 11$ for laboratory channels [65].

##### External sediment input/output

The source term $S_{b}$ represents additional sediment mass input or output
(sink) that can be defined on subsets of the computational domain. The source
can be specified as total volume flux including porosity
[$\mathrm{m}^{3}\mathrm{/}\mathrm{s}$]. Similarly to the hydrodynamic case
(§3.1.1), different approaches are adopted for the sediment sink, namely
exact, available and infinity.

####  3.2.2 Fixed bed concept

Morphodynamic simulations generate deposition and erosion patterns of the
riverbed. Erosion processes, if not limited, can proceed indefinitely in the
vertical direction. To account for the presence of non-erodible river bottom,
as in case of bedrock or concrete cover, a non-erodible fixed bed depth
$z_{rel}$ (Fig. 1a) can be set. This threshold also determines the volume of
sediment available for transport. The fixed bed elevation is defined relative
to the initial bottom elevation $z_{B}$ with $z_{rel}\leq 0$.

###  3.3 Scalar advection-diffusion

A number of environmental processes, such as pollutant, temperature or
nutrient transport, can be modelled assuming the passive advection and
diffusion of a scalar quantity, in the form of dissolved or particulated
particles [e.g. 27]. The scalar advection-diffusion module allows for the
simultaneous simulation of multiple passive species, up to a maximum of 5. The
transport of a generic species $c$ can be described by the following
advection-diffusion equation:

$\partial_{t}q_{c}+\partial_{x}\left[\dfrac{q_{x}q_{c}}{h}-h\left(K_{xx}\partial_{x}\phi_{c}+K_{xy}\partial_{y}\phi_{c}\right)\right]+\partial_{y}\left[\dfrac{q_{y}q_{c}}{h}-h\left(K_{yx}\partial_{x}\phi_{c}+K_{yy}\partial_{y}\phi_{c}\right)\right]=S_{\phi_{c}},\text{with}\,c=[1,5],$
(8)

where the unknown is $q_{c}$, the specific mass of the species $c$. It can be
expressed as $q_{c}=h\phi_{c}$, with $\phi_{c}$ the volumetric concentration
and $h$ the water depth, as usual. The term $S_{\phi_{c}}$ is a net source of
$c$ and $K_{ij}$ [$\mathrm{m}^{2}\mathrm{/}\mathrm{s}$] are the components of
the 2D diffusion tensor.

####  3.3.1 Scalar advection-diffusion closure relationships

For the scalar advection-diffusion module, the closure relationships are used
to model the contribution of external scalar input and output. In particular,
the term $S_{\phi_{c}}$ represents an additional scalar mass flux that is
added within a set of elements defined by regions. The source can be specified
either as an imposed concentration value or a total volumetric flux
[$\mathrm{m}^{3}\mathrm{/}\mathrm{s}$]. The behavior is analogous to the case
of hydro- and morphodynamics, §3.1.1 and §3.2.1.

The terms $K_{ij}$ of the diffusion tensor vary considerably with respect to
the physical nature of the transported species. Diffusive transport is
modelled in terms of both molecular diffusion $K^{m}$ and turbulent dispersion
$K_{ij}^{t}$, such that $K_{ij}=K^{m}I_{ij}+K_{ij}^{t}$, with $I_{ij}$ the
identity matrix. The molecular diffusion is assumed as an isotropic Fickian
process with constant coefficient $K^{m}$. Turbulent dispersion is non-
isotropic ($K^{t}_{ij}$) and scale with the friction velocity
$u_{*}=\lVert{\mathbf{u}}\rVert/c_{f}$ and water depth via a longitudinal
$\alpha_{L}$ and transversal $\alpha_{T}$ non-dimensional coefficient.
Suitable values for open channel flows in natural environments are
$\alpha_{L}$=13 and $\alpha_{T}$=1.2 [27].

##  4 Numerical solution

The numerical solution of the governing equations (1), (3) and (8) is sought
in a finite volume framework, with a spatial discretization based on
unstructured meshes (§4.1). For the temporal integration, an explicit first
order Euler scheme is used. In its basic configuration, the temporal
integration proceeds in a synchronous-decoupled way for all the modules,
meaning that the modules are independently integrated in time with the same
timestep (§4.5). The following sections detail the domain discretization
strategy and the adopted numerical solver for the fluxes calculation of the
three basic modules. The interested reader should refer to the provided
references for specific implementation details.

###  4.1 Domain discretization

The problem is discretised adopting a finite volume approach over unstructured
triangular meshes. A conforming triangulation $T_{\Omega}$ of the
computational domain $\Omega\subset\mathbb{R}^{2}$ by elements $\Omega_{i}$
such that $T_{\Omega}=\bigcup\Omega_{i}$, is assumed. Given a finite volume
element $\Omega_{i}$ (Fig. 2), $j=1,2,3$ is the set of indexes such that
$\Omega_{j}$ is a neighbour of $\Omega_{i}$. $\Gamma_{ij}$ is the common edge
of two neighbour cells $\Omega_{i}$ and $\Omega_{j}$, and $l_{ij}$ its length.
${\mathbf{n}}_{ij}=(n_{ij,x},n_{ij,y})$ is the unit vector which is normal to
the edge $\Gamma_{ij}$ and points toward the cell $\Omega_{j}$.

Figure 2: Sketch of the triangular discretization. Main notations adopted for
the generic computational cell $i$ and its $j$-th neighbour (with $j$=1,2,3).

###  4.2 Hydrodynamics

The system of governing equations (1) can be cast in vectorial form as

$\partial_{t}\boldsymbol{U}+\partial_{x}\boldsymbol{F}_{x}+\partial_{y}\boldsymbol{F}_{y}=\boldsymbol{S},$
(9)

where left-handside terms of (9) are

$\begin{array}[]{*{3}ccc}\boldsymbol{U}=\begin{bmatrix}H\\\ q_{x}\\\
q_{y}\end{bmatrix}&,&\boldsymbol{F}_{x}=\begin{bmatrix}q_{x}\\\
\dfrac{q_{x}^{2}}{h}+\dfrac{1}{2}gH^{2}-gHz_{B}\\\
\dfrac{q_{x}q_{y}}{h}\end{bmatrix}&,&\boldsymbol{F}_{y}=\begin{bmatrix}q_{y}\\\
\dfrac{q_{x}q_{y}}{h}\\\
\dfrac{q_{y}^{2}}{h}+\dfrac{1}{2}gH^{2}-gHz_{B}\end{bmatrix}\end{array}.$ (10)

The vector of source terms can be written as
$\mathbf{S}(\boldsymbol{U})=\mathbf{S}_{h}+\mathbf{S}_{fr}(\boldsymbol{U})+\mathbf{S}_{bed}(\boldsymbol{U})$,
where

$\begin{array}[]{*{3}ccc}\boldsymbol{S}_{h}=\begin{bmatrix}S_{h}\\\ 0\\\
0\end{bmatrix}&,&\boldsymbol{S}_{fr}=\begin{bmatrix}0\\\ -ghS_{fx}\\\
-ghS_{fy}\end{bmatrix}&,&\boldsymbol{S}_{bed}=\begin{bmatrix}0\\\
-gH\partial_{x}z_{B}\\\ -gH\partial_{y}z_{B}\end{bmatrix}\end{array}.$ (11)

By integrating the governing system of equations (9) in the control volume
$V=[\Omega_{i}]\times[t^{n},t^{n+1}]$, we obtain the general update formula
for the triangular element $i$:

$\boldsymbol{U}_{i}^{n+1}=\boldsymbol{U}_{i}^{n}-\frac{\Delta
t}{\left|\Omega_{i}\right|}\sum_{j=1}^{3}l_{ij}\left[\boldsymbol{F}_{ij}\right]+\Delta
t\boldsymbol{S}_{i}.$ (12)

Problem unknowns at cell $i$ and discrete time $n$ are represented by cell
averages $\boldsymbol{U}_{i}^{n}$; the numerical solution sought at time
$t^{n+1}=t^{n}+\Delta t$ is denoted by $\boldsymbol{U}_{i}^{n+1}$. In (12),
$\boldsymbol{F}_{ij}$ are the hydrodynamic fluxes estimated at the cell
interface $ij$, as in Fig. 2.

To compute the fluxes $\boldsymbol{F}_{ij}$ for the hydrodynamic system (9),
several well-established solvers are available. Here we adopt the well-known
HLLC approximate Riemann solver [66] which is a modification of the basic HLL
scheme to account for the influence of intermediate contact waves. Further
details on the HLLC approach are available in Chapter 10 of [18]. The solver
is proved to be robust and efficient in simulating unsteady flows and the
advection of passive tracers [27].

The numerical discretization of the three terms of the source vector
$\mathbf{S}$ (11) is conducted separately, according to the nature of each
term. The external inflow/outflow contribution $\mathbf{S}_{h}$ is added
explicitly to the continuity equation, as it is not a function of the problem
unknowns. The stiff friction source terms $\mathbf{S}_{fr}(\boldsymbol{U})$
are integrated with Runge-Kutta 2 [e.g. 67] in a semi-implicit fashion after
adopting a splitting technique. The implementation is analogous to the ones
proposed in [14, 68, 27]. The topographical terms
$\mathbf{S}_{bed}(\boldsymbol{U})$ are discretized using the modified-state
approach proposed by [69]. This results in an easy and robust treatment of
complex topographies and wetting and drying problems [27].

###  4.3 Morphodynamics

The Exner equation is solved in a synchronous-decoupled way with respect to
the shallow water problem (§4.2), meaning that the numerical integration of
the Exner equation (3) adopts the same integration timestep $\Delta t$ of the
hydrodynamic problem. The general update formula for the Exner problem reads:

$z_{Bi}^{n+1}=z_{Bi}^{n}-\frac{1}{1-p}\left[\frac{\Delta
t}{\left|\Omega_{i}\right|}\sum_{j=1}^{3}l_{ij}\left[q_{Bij}\right]+\Delta
tS_{bi}\right],$ (13)

with the same symbols introduced for (3). The term $q_{Bij}$ represents the
normal sediment flux at the cell interface $ij$ (Fig. 2).

For the numerical estimation of the term $q_{Bij}$, a number of approaches are
available in literature. In the current version, Basement implements an
Approximate Riemann Solver of HLL-type [sensu 18], as in [70]. The sediment
flux is thus calculated as

$q_{Bij}=\frac{\lambda_{s}^{+}q_{Bi}-\lambda_{s}^{-}q_{Bj}+\lambda_{s}^{+}\lambda_{s}^{-}(z_{Bj}-z_{Bi})}{\lambda_{s}^{+}-\lambda_{s}^{-}},$
(14)

where pedix $i$ ($j$) refers to quantities evaluated at the corresponding cell
(Fig. 2), and $\lambda_{s}^{+}$, $\lambda_{s}^{-}$ are speed estimation of the
morphological problem. We adopt the following closure [70]:

$\lambda_{s}^{-}=min(\lambda_{1i},\lambda_{1j})\;,\lambda_{s}^{+}=max(\lambda_{2i},\lambda_{2j}).$
(15)

The expression for the terms $\lambda_{1}$ and $\lambda_{2}$, calculated for
both cell $i$ or $j$ reads:

$\lambda_{1,2}=\frac{1}{2}\left(u_{n}-c\pm\sqrt{(u_{n}-c)^{2}+4\frac{\partial
q_{B,n}}{\partial q_{n}}c^{2}}\right),$ (16)

where $u_{n}$ is the normal velocity at the cell interface $ij$ (Fig. 2) and
$c=\sqrt{(}gh)$ is the so-called wave celerity.

###  4.4 Scalar advection-diffusion

The scalar advection-diffusion problem is solved in a synchronous-decoupled
way with respect to the shallow water problem (§4.2). We reformulate the
governing equation (8) via a Cattaneo-type relaxation technique, as proposed
by [27]. Two additional scalar conservation equations are then added to (8):

$\partial_{t}\psi_{x}^{c}-\partial_{x}\frac{\phi_{c}}{\zeta}=-\frac{\psi_{x}^{c}}{\zeta}\,,\;\;\;\;\partial_{t}\psi_{y}^{c}-\partial_{y}\frac{\phi_{c}}{\zeta}=-\frac{\psi_{y}^{c}}{\zeta}$
(17)

where the new symbols are $\zeta$, a positive and small relaxation time, and
$\psi_{x}^{c}$ and $\psi_{y}^{c}$ are two auxiliary variables that recover
$\partial_{x}\phi_{c}$ and $\partial_{y}\phi_{c}$, respectively, for a
sufficiently small $\zeta$ [27]. After a trivial substitution of
$\psi_{x}^{c}\approx\partial_{x}\phi_{c}$ and
$\psi_{y}^{c}\approx\partial_{y}\phi_{c}$ into (8), the system composed by (8)
and (17) can be rewritten in vectorial form as

$\partial_{t}\boldsymbol{U}+\partial_{x}\boldsymbol{A}_{x}+\partial_{y}\boldsymbol{A}_{y}+\partial_{x}\boldsymbol{D}_{x}+\partial_{y}\boldsymbol{D}_{y}=\boldsymbol{S}_{c}+\boldsymbol{S}_{rel}$
(18)

where the vectors $\boldsymbol{U}$, $\boldsymbol{A}_{x}$,
$\boldsymbol{D}_{x}$, $\boldsymbol{S}_{c}$ and $\boldsymbol{S}_{rel}$ read

$\begin{array}[]{*{4}cccc}\boldsymbol{U}=\begin{bmatrix}q_{c}\\\
\psi_{x}^{c}\\\
\psi_{y}^{c}\end{bmatrix},&\boldsymbol{A}_{x}=\begin{bmatrix}\dfrac{q_{c}q_{x}}{h}\\\
0\\\
0\end{bmatrix},&\boldsymbol{D}_{x}=\begin{bmatrix}-h(K_{xx}\psi_{x}^{c}+K_{xy}\psi_{y}^{c})\\\
-\frac{q_{c}}{\zeta h}\\\
0\end{bmatrix},&\boldsymbol{S}_{c}=\begin{bmatrix}S_{\phi_{c}}\\\ 0\\\
0\end{bmatrix},&\boldsymbol{S}_{rel}=\begin{bmatrix}0\\\
-\frac{\psi_{x}^{c}}{\zeta}\\\
-\frac{\psi_{y}^{c}}{\zeta}\end{bmatrix},\end{array}$ (19)

with $\boldsymbol{Q}$ representing the conserved scalar quantities,
$\boldsymbol{A}_{x}$ is the advective fluxes vector, $\boldsymbol{D}_{x}$ is
the diffusive-relaxed fluxes vector, both in x direction. The scalar source
terms are $\boldsymbol{S}_{c}$, whilst the source terms arising form the
relaxation are $\boldsymbol{S}_{rel}$. For brevity, we omit the formulation
for the y direction ($\boldsymbol{A}_{y}$ and $\boldsymbol{D}_{y}$), which is
analogous. The interested reader can refer to [27] for a step-by-step
derivation.

The scalar fluxes in (19) are solved through the SVT solver introduced by
[27]. The scheme presents a flux-splitting approach combining the advective
and diffusive-relaxed fluxes and employs different solvers for each. The HLLC
solver, applied also for the hydrodynamic fluxes (§4.2) provides the advective
component of the scalar fluxes at the cell interface $\boldsymbol{A}_{ij}$.
For the diffusive-relaxed component, the SVT technique derives the fluxes at
the interface $\boldsymbol{D}_{ij}$ directly from the Riemann invariants of a
two non-linear wave Riemann problem.

Similarly to the hydro- and morphodynamic problems, the control volume
$V=[\Omega_{i}]\times[t^{n},t^{n+1}]$ is used to integrate the governing
system (18), in order to obtain the following scalar update formula at the
element $i$

$\boldsymbol{U}_{i}^{n+1}=\boldsymbol{U}_{i}^{n}-\frac{\Delta
t}{\left|\Omega_{i}\right|}\sum_{j=1}^{3}l_{ij}\left[\boldsymbol{A}_{ij}+\boldsymbol{D}_{ij}\right]+\Delta
t(\boldsymbol{S}_{c}+\boldsymbol{S}_{rel}),$ (20)

where the fluxes $\boldsymbol{A}_{ij}$ and $\boldsymbol{D}_{ij}$ are computed
at each cell interface $ij$ (Fig. 2).

The numerical integration of the two source terms vectors is conducted
separately, according to the nature of the terms. The simple scalar sources
$\boldsymbol{S}_{c}$ are computed with a first-order Euler scheme, while the
stiff relaxation source terms $\boldsymbol{S}_{rel}$ are integrated by means
of a locally implicit Euler method.

###  4.5 Stability condition

Numerical integration proceeds with a dynamic timestep $\Delta t$, evaluated
at each time loop (Fig. 6b) that fulfills the well-known Courant-Friedrichs-
Lewy stability conditions [18]. In the current implementation, the condition
is expressed as:

$\Delta t=CFL\min_{1\leq i\leq N}\left(\min_{1\leq j\leq
3}\left(\frac{\rho_{ij}}{{\lambda_{ij}}}\right)\right)\;,$ (21)

where $\rho_{ij}$ is twice the distance between the edge $j$ and the centroid
of the cell $i$ (Fig. 2), and $N$ is the total number of domain elements. The
term $\lambda_{ij}$ is an estimation of the largest eigenvalue of the
hydrodynamic problem (1), namely $\lambda_{ij}=|u_{n}|+\sqrt{gh}$ with the
simbology already introduced. The CFL coefficient ranges between 0 and 1: by
default it is set to 0.9, if not specified otherwise.

##  5 Initial and boundary conditions

###  5.1 Initial conditions

All modules require the user to define the initial conditions of the
simulation. Two types of initial conditions are similarly available for all
the modules:

  * 1.

region defined: user explicitly defines the initial values of the problem
unknowns (e.g. water depth and specific discharge for hydrodynamics).
Different values can be assigned to different region of the computational
domain;

  * 2.

continue: values are taken from the result file of previous simulations.

In addition, the hydrodynamic module allows also to set dry conditions (no
water in the domain) as initial conditions. In this case, the domain will
progressively fill with water, in relation to the assigned inflow boundary
conditions or internal sources (§3.1.1).

###  5.2 Boundary conditions

The boundary conditions (hereinafter BCs) have different specifications for
each core module (see following Sections), but they all classify in three
common types: external standard, external linked and internal BCs.

Fig. 3 exemplifies the main concepts adopted for the BCs. The computational
domain $\Omega$ is defined by the domain boundaries, as $\Gamma_{1,2,3}$. An
external standard BC is dependent only on the local flow conditions and on
some user-defined rules. This represents the most common case, for example to
define impermeable walls or river inflows and outflows. By default, all the
external boundaries are all set as wall. The wall BC consists of a fixed,
frictionless, reflective impermeable wall. In external linked BCs instead, the
local BCs are defined also with information from a linked boundary. Typical
example is a weir, where the flow discharge at the downstream side of the weir
depends on the water stage on the upstream side.

The third type of BCs, internal, are defined within the computational domain
$\Omega$, and not at the edges (Fig. 3). This BC type comes in handy in case
of very large domain application, because it allows to test different
configurations of hydraulic structures (e.g. different locations of a weir or
training wall), without the need of regenerate the entire computational mesh
for every configuration.

A summary of the main feature of the BCs for the three core modules follows
here. The interested reader can refer to the official documentation for
further details.

Figure 3: Example of modelling domain with different types of boundary
conditions. The computational domain $\Omega$ can include the river channels
but also the surrounding floodplains. The domain is delimited by the external
BCs: impermeable walls (default type, $\Gamma_{1,2,3}$) are depicted in solid
black, while standard inflow and outflow are in dashed blue. A weir (dotted
brown) is modelled with an external linked BC. A training wall (thick solid
green) is modelled with an internal BC.

####  5.2.1 Hydrodynamic BCs

The hydrodynamic module implements different types of BCs, with a different
level of customization. Depending on the BC type, user-assigned data is
requested, as single constant value in time (e.g. lake level, constant
discharge), as time series (e.g. hydrograph), or as set of parameters
describing a dynamic behaviour (e.g. weir activation rule). In particular:

  * 1.

Standard BCs: in addition to wall BC, inflows (upstream BCs) and outflows
(downstream BCs) can be assigned. As standard inflows, three options are
provided (namely uniform, explicit, zhydrograph). For all cases, given a total
volume discharge $Q$ [$\mathrm{m}^{3}\mathrm{/}\mathrm{s}$] or water surface
elevation [m] as input, the inflows condition for the mathematical unknowns of
system (1) is set. For the standard outflows a value for the water depth $h$
must be specified. Possible options are: uniform conditions, hydraulic weir,
rating curve, hydrograph and zerogradient (i.e. Neumann BC). It is worth
remarking that the specific type of upstream and downstream BCs should be
selected depending on the local flow conditions (i.e. sub- or super-critical).

  * 2.

Linked BCs: this type of boundaries establish a link between two certain
region of the domain where the governing equations are not solved. This type
of BCs are particularly designed to simulate the behaviour of hydraulic
structures within the river channel, such as weirs, gates, bridges, spillways.

  * 3.

Internal BCs: they are fictitious boundaries defined as segments at the
interfaces of some computational cells. On these segments, three different
conditions can be enforced, instead of the solution of the SWE (1). Options
are: static walls, dynamic walls and rating curve. With the static wall, the
standard wall condition is applied on both sides of the internal boundary.
With the dynamic wall, the wall conditions are applied until reaching a given
threshold value (time or water depth) after which the wall is removed, and the
SWE are solved. With the rating curve option (or h-Q relation) a give flow
relation is applied on one side of the internal boundary, while on the other
side, wall conditions apply (unidirectional flow). The internal BCs are
particularly useful when simulating, for example, the collapse of hydraulic
structures: after the collapse, the simulation can proceed by calculating the
actual free unsteady flow (i.e. governing equations) over the (former)
boundary.

####  5.2.2 Morphodynamic BCs

The sediment flow is defined as a specific bedload flux, which is averaged and
evenly distributed at the domain boundary conditions over the boundary length.
In analogy with the hydrodynamic module, the morphodynamic boundaries are of
type external standard and linked.

  * 1.

Standard BCs: for the upstream BCs, Basement implements three versions that
allow to simulate: i) a given input of sediment as time series (i.e.
sedimentograph), ii) a sediment input derived from the hydrodynamic conditions
under transport capacity conditions or iii) bed equilibrium condition, where
the upstream bed elevation is kept constant. Two downstream BCs are available,
allowing the simulation of i) equilibrium condition and ii) check-dam. In this
second option an equilibrium boundary condition is activated only if the bed
level reaches a given threshold value, otherwise a wall type boundary is
assumed.

  * 2.

Linked BCs: one BC is available. It allows for the simulation of sediment
transport through given hydrodynamic linked conditions, hence to ensure
sediment continuity in the simulated channel.

####  5.2.3 Scalar advection-diffusion BCs

Scalar BCs are defined in terms of concentration of total volumetric rate
[$\mathrm{m}^{3}\mathrm{/}\mathrm{s}$], evenly distributed throughout the
length of the relevant domain boundary. The implemented types are:

  * 1.

Standard BCs: three types are available. i) scalar inflow as a constant value;
(ii) scalar inflow as a time-series and (iii) zerogradient (i.e. Neumann BC)
outflow.

##  6 Software design

###  6.1 Modelling workflow

The standard modelling procedure involves three phases: the pre-processing
phase, the numerical simulation phase and the post-processing phase (Fig. 4).
Basement is designed to integrate into this workflow. Moreover, the entire
workflow relies on open-source or freeware tools. In the following we list the
phases and provide a short description of the different configuration and
results file formats as used by Basement.

  1. 1.

Pre-processing: in this phase the user is required to define the model domain
and the input data. The mesh file (customized 2dm format, "MyMesh.2dm" in Fig.
4) contains the description of the triangular unstructured computational mesh.
The file can be generated with BASEmesh, a Python script as well as a QGIS
plugin (see software website for details), or via grid generator software that
supports the 2dm format. In addition, further input data such as time series
of water and sediment discharge (or other quantities) to be used as BCs can be
provided (ASCII format, "MyData.txt" in Fig. 4).

  2. 2.

Numerical simulation: the actual simulation can be run via either CLI or GUI.
In Basement version 3, the numerical simulation is split into three steps
(description follows in §6.2). The final simulation results are stored in a
general purpose binary container (Hierarchical Data Format HDF5,
www.hdfgroup.org). Basement generates also an XDMF file (eXtensible Data Model
and Format, http://www.xdmf.org) which contains a machine-readable description
of the data stored in the HDF5 file.

  3. 3.

Post-processing: the XDMF file ("output.xdmf" in Fig. 4) can be opened with
the Crayfish plugin for QGIS or with Paraview for final results visualization
and further post-processing. In addition, ad-hoc Python scripts can be used to
manipulate results directly from the binary container (some scripts are
provided at software website).

Figure 4: Modelling workflow. Pre-processing: generation of the computational
mesh from topographical data and definition of input time series. Numerical
simulation: Basement. Post-processing: elaboration and visualization of the
results.

###  6.2 Simulation steps

The numerical simulation phase consists of three steps: the pre-simulation,
the simulation, and the post-simulation (Fig. 5). Each step can be completed
by running a corresponding Basement executable via the Graphical User
Interface or Command Line Interface. This modular design allows a
customization of the simulation workflow by the user and an efficient batch
processing of Basement steps. For instance, the programs can be run from a
scripting language like Python.

The different executables are configured using a dedicated command file, as
detailed in the following list. These command files use the standardized JSON
file format (JavaScript Object Notation) (Fig. 5). The Basement GUI is
designed to support the user with creating the command files, running and
monitoring the three simulation steps. In particular, the GUI validates the
configuration parameters and automatically adds required parameters where
default values are available.

  1. 1.

The pre-simulation step focuses on the model definition. In particular, the
"model.json" command file contains: i) physical properties, ii) initial
conditions and iii) boundary conditions of the physical problem, and further
iv) numerical parameters. The setup executable first reads the computational
mesh "MyMesh.2dm", the external required data "MyData.txt" and the command
file "model.json". Then it validates and stores the model inside the binary
container "setup.h5".

  2. 2.

The simulation is carried out on a selected computational backend (§6.4). It
is driven by the command file "simulation.json" that contains the simulation
parameters (e.g. execution time, output time and desired output quantities).
The program reads and executes the model "setup.h5" generated in the previous
step. The results of the simulation are stored in a second binary container:
"results.h5".

  3. 3.

The post-simulation step is configured using a command file "results.json"
that contains the selected output format (currently only XDMF is supported).
The output is then available for the post-processing phase. (§6.1).

Figure 5: Software components and simulation steps. The Basement software is
composed of a set of executables (red rectangles) driven by JSON configuration
files (grey labels). Data is stored in HDF5 containers (green cylinders). In
dashed arrows: the special actions of simulation re-run and restart.

When it is necessary to run a new simulation starting from the results of a
previous one, two options are available: Restart and Re-run. When performing a
Restart the pre-simulation step is executed again, i.e. a new model is
generated from scratch with (potentially) a new set of parameters. The user
indicated an existing "results.h5" file that is to be used to fetch the
initial conditions for the new model. The Re-run option does not generate a
new model, but uses the existing model ("setup.h5") with initial conditions
taken from the current results file. In this scenario the user can only modify
the simulation and results parameters (i.e. duration, output timestep and
output type), but not the the model parameters. This second option is
particularly useful in the case of large models (i.e. millions of
computational cells): in this situation the pre-simulation step can be
computationally demanding (i.e. tens of minutes). If the user only needs to
extend the simulation duration, for example, then the Re-run option allows to
skip the pre-simulation step.

###  6.3 Modularity and Sequencing

Basement aims to simulate different river processes with a high level of
flexibility and efficiency: with this in mind, we designed the software
adopting a modular approach. The two core concepts of this design are Modules
and Kernels, described here below.

Modules take care of the simulation of specific river processes (e.g. module
"Hydrodynamics" in Fig. 6a). They can be nested to simulate processes with an
increasing level of detail/complexity (e.g. module "HYD External Source", Fig.
6a). Modules are activated by the user in the pre-simulation step. An
activated module triggers the execution of a number of Kernels throughout the
simulation.

A Kernel is a set of operations to be executed on each entity (e.g. a cell $i$
or an edge $j$, Fig. 2) of the computational domain or a subset of it.
Depending on the specific task, Kernels can be scheduled for a single
execution (i.e. initialization Kernels) or for repeated execution in each
iteration of the integration time loop (Fig. 6b). While the time loop is
executed with a timestep $\Delta t$ that satisfies certain stability
conditions (§4.5), the Kernels can be scheduled for execution at different
(larger) time intervals to reflect the nature of the simulated processes.

The architecture based on Modules and Kernels has two main advantages. First,
it is flexible in that it allows users (and software developers) to easily add
or remove specific modules without interfering with other existing modules. In
particular, this permits an integration of further modules as development
continues (§8). Second, it is efficient, because only the necessary Kernels
are scheduled for execution at setup time (pre-simulation step).

Figure 6: Examples of modules and associated kernels. (a) At model setup,
depending on the simulated physical problem, the user triggers the activation
of a set of modules. (b) The active modules trigger a unique kernel sequence
to be executed at simulation time to correctly simulate the requested
processes. Each module corresponds to a different set of kernels.

###  6.4 Parallelization Strategy and Computational Backends

The parallelization strategy of the Basement numerical core addresses two main
aspects: i) the use of different technologies (i.e. computational backends)
generated from the same, unique software source code. This allows for an
easier source code maintenance and integration of future/different backends.
ii) An efficient and heavy parallelization of the numerical core following the
concept of data parallelism. To this end the numerical core of Basement
integrates OP2 [30, 43], which is an open-source framework for the development
of unstructured grid applications. Using source-to-source translation, OP2
generates the appropriate code for different target platforms by introducing
an additional level of abstraction between the numerical algorithm and its
execution. It supports multi-core CPUs, GPUs, and even clusters via MPI
(Message Passing Interface, http://www.mpi-forum.org).

Basement currently supports multi-core CPUs and GPUs. When starting the
simulation, the user can select to compute on the CPU, the GPU, or a
combination of both. All the currently supported backends (Table 3) are
available for both Windows and Linux (Ubuntu) operating systems. It is
important to note that the choice of graphics processing units is currently
limited to Nvidia (CUDA) cards. The precise requirements are provided in the
official documentation. All the backends can execute the numerical simulations
in double (default) or single precision, with different performance
characteristics (§7.7).

Table 3: Description of available backend types for Basement v3. Type | Description
---|---
seq | sequential execution on the CPU
omp | multi-threading using OpenMP technology
cuda | GPU
cudaC | GPU with some kernels running sequentially on the CPU
cudaO | GPU with some kernels running in parallel (OpenMP) on the CPU

##  7 Results

A set of selected test cases (T1-T6) are proposed here to test the robustness,
accuracy and efficiency of the three basic modules. Table 4 summarizes the key
features of each test case. The interested reader can refer to the official
software documentation for further examples. Finally, §7.7 focuses on the
software performance and scalability. All the test cases are freely available
at (link provided after paper acceptance).

Table 4: Main features of the adopted benchmarks. ID | Module | Comparison | Key features
---|---|---|---
T1 | hydrodynamic | field data |  i) highly unsteady shock wave generation and propagation, ii) wet-and-dry processes, iii) performance of the hydrodynamic module
T2 | morphodynamic | lab data |  i) sediment transport with a transcritical 1D flow, ii) upstream BCs: uniform inflow, iii) downstream BCs: imposed water level
T3 | morphodynamic | lab data |  i) 2D dam-break in a complex domain, ii) sediment transport with an advancing wet-and-dry front, iii) downstream BCs: free outflow, iv) performance of the morphodynamic module
T4 | morphodynamic | lab data |  i) erosion and deposition in a channel bend, ii) sediment transport direction correction, iii) upstream BCs: unsteady hydrograph
T5 | scalar advection-diffusion | numerical sol |  i) 1D strong rarefaction waves, ii) conservation of a steady discontinuity of the scalar quantity
T6 | scalar advection-diffusion | numerical sol |  i) 2D complex domain, ii) advection and diffusion of two different scalar quantities, iii) conservation and mixing of the scalars, iv) performance of the scalar advection-diffusion module

###  7.1 T1: Malpasset dam collapse

The scope of this test is to assess the robustness and accuracy of the
hydrodynamic solver when simulating a shock-type hydrodyamic wave travelling
on a highly irregular and dry domain. The collapse of the Malpasset dam, in
the Reyran River Valley (Fréjus, France), represents a well-established
hydrodynamic benchmark for numerical models [e.g. 71, 72, 73]. In 1959, the
66.5 m high dam collapsed almost instantaneously, generating an up to 40 m
high flood wave that propagated down the Reyran valley, destroyed the two
villages Malpasset and Bozon and reached the Mediterranean Gulf 21 minutes
later [72]. The propagation of the flood wave was reconstructed via the
maximum water level and the flood arrival time, recorded at multiple
locations. In particular, the maximum water level is available from a police
survey for 17 survey points, marked as P1 to P17 in Fig. 7 and the flood
arrival time is known from three electric transformer stations which have been
destroyed by the flood wave. The locations of the transformer stations are
indicated as A, B and C in Fig. 7. Coordinates and recorded arrival times are
listed in Table 5. We make use of such field data to test the performance of
the hydrodynamic module.

The computational domain is discretized with 499,059 triangular elements. The
domain boundaries are set to walls (§5.2.1), with exception of the downstream
boundary located in the Mediterranean Gulf, where a fixed water level was set
to 0 m. The initial conditions are a fixed water surface elevation of 100 m in
the reservoir, and dry conditions in the rest of inland domain. The initial
velocity was set to $0.0\text{\,}\mathrm{m}\mathrm{/}\mathrm{s}$ in the entire
domain. In accordance with [71], the Manning’s friction coefficient was set to
$0.033\text{\,}\mathrm{m}^{-1/3}\mathrm{s}$ for the whole domain. The CFL
number was set to 0.9.

Figure 7: T1: planar view of the Malpasset test case. Visualization of the
computational domain illustrating the bottom elevation (gray scale) and the
initial water depth (blue scale). Letters indicate the location of the
transformer stations (A-C), and the maximum water level survey points
(P1-P17). The left box shows a magnification of the computational mesh, as
reference.

The simulated maximum water levels are compared to the 17 field observations
in Fig. 8, with overall good agreement. The average relative error is 7.15%,
with the largest observed at P13, with an overestimation of 30.6%. To
highlight the effects of the topographical approximations of the Digital
Elevation Model (and hence of the computational mesh), we compared recorder
and simulated water level values as follows. For each punctual maximum water
level recorded in field observations (blue triangles in Fig. 8), we compared
the maximum simulated values of the spatial mean, maximum and minimum among
the computational cell containing the observation point and its three
neighbours (black and red series in Fig. 8). We expect lower discrepancies
between recorded and simulated values where the numerical values, hence the
topographical elevations, are spatially homogeneous. As a matter of fact,
points P1, P7 and P13 (Fig. 8) have the large discrepancy between measured and
simulated values, but also the largest spatial variability of the numerical
values (red shaded area). This suggests that such discrepancies relates more
to the local topographical approximations of the DTM rather than to the
numerical model.

Figure 8: T1: Malpasset dam-break wave maximum water level. Numerical values
are compared to the field data at the survey points P1 to P17: simulated
values are given as maximum, minimum (red lines) and mean value (black line)
of the computational cell containing the survey point and its three
neighbours.

Observed and simulated times of flood arrival are given in Table 5. Simulated
values are in good agreement with measured ones for all the electrical
transformer stations (ET). Simulated arrival times have a maximum relative
error of 3.8% for ET B, corresponding to an absolute delay of 47 s. It is
worth mentioning that the friction value influences the simulated arrival
times.

Table 5: T1: Malpasset dam-break wave arrival times. Observed and simulation time of flood arrival (TFA) and relative error (Err) for the three electrical transformer stations (ET) destroyed by the flood wave. ET | x | y | TFAobs | TFAsim | Err
---|---|---|---|---|---
0 [-] | [m] | [m] | [s] | [s] | [%]
A | 5550 | 4400 | 100 | 103 | 3
B | 11900 | 3250 | 1240 | 1287 | 3.8
C | 13000 | 2700 | 1420 | 1435 | 1

###  7.2 T2: Propagation of a sediment bore

Scope of the test is to assess the robustness of the de-coupled hydro-
morphodynamic solver approach, particularly when simulating the sediment
transport over a transcritical flow. This represents a critical test,
especially when adopting de-coupled approaches [e.g. 74]. Moreover, the
simulation tests the morphological solver capability in well reproducing the
dynamics of an advancing sediment bore.

In this test case, the flume experiment proposed in [75, run 2] is reproduced
numerically. The computational domain is a composed by a straight 6.9 by 0.5 m
channel, representing the lower part of the original experimental flume, and
is discretized with 24,612 triangular elements. The sediment has a
characteristic diameter of 1.65 mm. The water and sediment discharge at the
upstream boundary are set to
$0.012\text{\,}{\mathrm{m}}^{3}\text{\,}{\mathrm{s}}^{-1}$ and
$0.196\text{\,}{\mathrm{m}}^{3}\text{\,}{\mathrm{s}}^{-1}$ (with porosity)
respectively. The flume is at initial uniform flow conditions, characterized
by a supercritical flow. At $t$=0 s, a fixed water level of 0.2093 m is
imposed at the downstream boundary and sediment transport out of the domain is
stopped. This results in the formation of an hydraulic jump moving upstream in
the flume, and a subsequent downstream propagation of a sediment bore. The CFL
number is set to 0.9, the bed porosity is assumed constant and equal to 0.42,
and the simulation duration is 500 s.

Figure 9: T2: propagation of a sediment bore. Initial (dashed lines) and final
(solid lines) longitudinal profiles of bed elevation (black) and water
elevation (blue) for the propagation of a sediment bore test.

Fig. 9 shows the initial and final profiles of the simulated bed and water
elevations. The solver reproduces well the sharp transition between super- and
sub-critical flow conditions. The position of the sediment front in time is
shown in Fig. 10, with good agreement between simulated and experimental
values.

Figure 10: T2: evolution in time of the sediment front position. Blue
triangles are the experimental values [75], while the black line is the
numerical solution.

###  7.3 T3: Dam-break over a mobile bed with a sudden enlargement

Scope of the test is to assess code robustness in simulating sediment
transport at wet-dry interface, and the accuracy in reproducing
scour/deposition patterns. The experiment illustrated in [76] represents a
well-know morphodynamic test for numerical models [e.g. 14, 70, 77]. The
domain consists of a flat flume with a non-symmetrical sudden enlargement
(Fig. 11). The bed is composed of a coarse uniform sand with a median diameter
of $d_{m}$=1.82 mm. The initial conditions are defined by an horizontal layer
of fully saturated sand of thickness 0.1 m over the whole domain and an
initial water storage of depth 0.25 m upstream of the dam, located at section
$x$=3.0 m. At time $t$=0 s, the dam is suddenly removed, resulting in the
propagation of a dam break wave with consequent sediment transport.

Figure 11: T3: planar view of the dam-break over mobile bed setup (deformed
axis). Experimental and numerical results are compared at survey points P1 to
P6 during the simulation, and at cross-sections CS1 and CS2 at the end of
simulation.

The computational domain is discretized by unstructured triangular cells at
different resolutions (follows in Tab. 6). Inviscid wall boundary conditions
are set at the upstream and lateral domain boundaries, while a free-outflow
condition is used at the downstream outlet. The Manning coefficient is set to
$0.0167\text{\,}\mathrm{m}^{-1/3}\mathrm{s}$, the sediment density and
porosity are set to
$2680\text{\,}\mathrm{k}\mathrm{g}\mathrm{/}\mathrm{m}^{3}$ and 0.47,
respectively. The sediment transport is evaluated with the MPM-like formula
(Table 2), setting $\theta_{cr}$= 0.0495 for the critical Shields stress,
$\alpha$ = 3.97 and $m$ = 1.5 for the remaining parameters [56]. The CFL
number is set to 0.9. The numerical simulations last 12 seconds.

The evolution of the water elevation during the simulation is shown in Fig. 12
for the six survey points. The simulated series show a fairly good agreement
with the experimental values: the dam break wave arrival time is well captured
and the maximum elevation values are comparable with the measured ones.
Moreover, the simulated series show minor discrepancies with the experimental
ones after the arrival of the first wave. As already pointed out by previous
works [14, 78], discrepancies are due to the extremely complex flow pattern
generate by multiple wave reflections while simulation proceeds in time, which
potentially generate tri-dimensional flow structures. Nevertheless, obtained
series are coherent with the ones of [14], where a second-order accuracy model
was employed.

Figure 12: T3: experimental and numerical water elevation at six survey
points. Sub-panels represent the survey points P1 to P6 as in Fig. 11:
experimental points are given as full circles, whilst numerical results as
lines, with four different mesh sizes (27k, 54k, 108k, 216k computational
cells).

Numerical bed elevations after 12 s are compared to the experimental results
in Fig. 13. The simulated scour and deposition patterns are well reproduced.
The magnitude of the scour at cross section CS1 ($y\approx 0.25$ m) matches
well, with an underestimation of the deposition pattern ($y\approx 0.35$ m).
At cross section CS2, the simulated deposition magnitude matches well with the
experimental one, but with a small shift toward the lateral boundary.

Figure 13: T3: experimental and numerical bed elevation at different cross-
sections. Sub-panels represent cross-section CS1 and CS2 as in Fig. 11:
experimental points are given as full circles, whilst numerical results as
lines, with four different mesh sizes (27k, 54k, 108k, 216k computational
cells).

###  7.4 T4: Scour and deposition on a channel bend

The scope of the simulation is to test the correct reproduction, both in term
of positioning and magnitude, of a river point bar generated by a channel
bend. In this test we numerically reproduced one experiment from [79], already
adopted as morphodynamic benchmark test [e.g. 80]. The flume is U-shaped, with
a bend of $None$ having a costant radius along the center line of $R_{c}$=4 m.
The cross section is rectangular with width $W$=1 m and slope $S$=0.2%. The
two straight reaches before and after the bend are 11.5 m long. The median
diameter of the bed material was $d_{m}$=1 mm. In the experimental run the
flume was fed with a simplified (triangular) flood hydrograph, having a base
flow of $0.02\text{\,}{\mathrm{m}}^{3}\text{\,}{\mathrm{s}}^{-1}$, and peak
flow of $0.053\text{\,}{\mathrm{m}}^{3}\text{\,}{\mathrm{s}}^{-1}$. The rising
and falling limbs last 100 and 200 minutes, respectively. Afterwards a
constant baseflow was kept for another 100 minutes. During the experimental
run, a steady point bar in the inner side of the bed develops and grows, with
a corresponding erosion on the outer side.

In the numerical setup, the domain is discretized with 24,523 computational
cells. We set the porosity to 0.4, and used the MPM-like formula with
parameters from [56], as in Table 2. The lateral slope factor $N_{l}$ is equal
to 1.4 and the curvature factor $N*$ is set to 11 (6 and 7). At the numerical
domain boundaries uniform flow and equilibrium sediment transport conditions
are imposed. The simulation, as the experimental run, lasts 400 minutes.

A planar view comparison between numerical and experimental run is depicted in
Fig. 14. The final bed variation with respect to the initial flat
configuration is scaled with the approaching (i.e. upstrem reach) flow depth
$h_{0}$. The magnitude of scours and depositions for the numerical run ranges
between -0.75 and 0.75, matching fairly well with the experimental values.
Also the positioning of the point bar, with the maximum deposition
anticipating the middle of the bend ($None$) is well reproduced numerically.

Figure 14: T4: scour and deposition on a channel bend. Planar view of relative
bed change $\Delta$ $z_{b}$/ $h_{0}$, with $h_{0}$ approaching flow depth.
Numerical results on the left, laboratory results from [79] on the right.

Fig. 15 shows the cross-sectional profile of the relative bed change for the
numerical simulation and the experimental run in the middle of the flume bend
(at $None$). The numerical profile reproduces well the experimental trend.
This test demonstrates the software capability in simulating an unsteady
morphological building process, i.e. a point bar development during a flood.
Such process can be well reproduced only by implementing suitable corrections
of the sediment transport direction due to gravity and curvature, as presented
in §3.2.1.

Figure 15: T4: scour and deposition at one cross-section. Cross-sectional view
(at $None$) of relative bed change $\Delta$ $z_{b}$/ $h_{0}$, with $h_{0}$
approaching flow depth. In the $x$-axis the radial coordinate $R$ is scaled
with the center-line curvature radius $R_{c}$; black line is the numerical
solution, blue triangles are from the laboratory experiment of [79].

###  7.5 T5: Steady scalar discontinuity with two diverging hydrodynamic waves

The test assesses the correct advection of scalar concentration. This is
assessed with a challenging test: a steady discontinuity of a scalar
concentration subjected to strongly variable flow conditions. The chosen test
is an idealized one-dimensional problem, but nevertheless it is particularly
challenging for a pletora of numerical scheme [18]. The domain is a simple,
flat-bed, channel 100 m long and 0.1 m wide. The domain is deliberately chosen
very narrow, to mimic a 1D setup, given that the exact solution of the problem
is available in one-dimension case. As initial conditions, the water depth $h$
is set even in all the domain, whilst the initial longitudinal specific
discharge $q_{x}$ and the concentration of a generic scalar $\phi$ present a
discontinuity:

$\left\\{\begin{aligned} q_{x}&=-3.0\leavevmode\nobreak\
$\mathrm{m}^{2}\mathrm{/}\mathrm{s}$&\text{if }x&<50\leavevmode\nobreak\
$\mathrm{m}$,&q_{x}&=3.0\leavevmode\nobreak\
$\mathrm{m}^{2}\mathrm{/}\mathrm{s}$&\text{\>otherwise},\\\ \phi&=1.0&\text{if
}x&<50\leavevmode\nobreak\ $\mathrm{m}$,&\phi&=0.0&\text{\>otherwise},\\\
h&=1.0\leavevmode\nobreak\ $\mathrm{m}$&\forall x&.\end{aligned}\right.$ (22)

The domain is discretized with 1362 computational cells, lateral walls are
reflective and inviscid, whilst transparent boundary conditions are set at
beginning and end of the channel. The CFL is set to 0.95 and the simulation
timeout is $t$=2.5 s. As the simulation starts, two strong rarefaction waves
start diverging from the center of the domain towards the two extremities,
suddenly forming a water depression in the center. Despite the strong
unsteadiness of the hydrodynamic quantities during the simulation, a steady
contact wave persists in the domain, avoiding the scalar quantity to mix in
the domain.

The numerical solution at simulation timeout is compared with the exact
solution of the problem in Fig. 16. The hydrodynamic exact solution is
obtained by resolving the two-rarefaction Riemann Problem [18], whilst the
exact solution for the scalar advection is identical to the given initial
conditions. Fig. 16 underlines how the numerical solution correctly
approximates the exact solution in all the domain sections. The scalar
discontinuity is perfectly maintained throughout the simulation, confirming
the accurate resolution of the steady contact wave.

(a)  (b)

Figure 16: T5: comparison between numerical and exact solution of a steady
contact wave. Solution is given at timeout $t$=2 s. Exact solutions are
depicted with lines, whereas numerical values are symbols, decimated for the
sake of visualization. Panel (a) shows the water depth (blue diamonds) and the
scalar concentration (red circles), panel (b) shows the specific discharge on
the x direction (green circles).

###  7.6 T6: Scalar advection and diffusion in a dam-break over a complex
domain

With this test we assess the solver capability in correctly preserving the
liquid and scalar species mass when simulating the advection and diffusion of
species during a dam-break phenomena. The test is particularly harsh, due to
the presence of fast wetting-drying fronts and multiple discontinuous flow
regions.

The test is an ad-hoc setup inspired to a common benchmark for hydrodynamic
codes [e.g. 27, 81]. The domain is composed by a rectangle
$[0;75]\times[-15;15]$ m. The bottom $\eta(x,y)$ is fixed during the
simulation, and defined as

$\displaystyle\eta(x,y)$
$\displaystyle=\mathrm{max}\left(0,\eta_{1}(x,y),\eta_{2}(x,y),\eta_{3}(x,y)\right)\text{,
with}$ (23) $\displaystyle\eta_{1}(x,y)$
$\displaystyle=1-\frac{1}{8}\sqrt{(x-30)^{2}+(y+9)^{2}},$
$\displaystyle\eta_{2}(x,y)$
$\displaystyle=1-\frac{1}{8}\sqrt{(x-30)^{2}+(y-9)^{2}},$
$\displaystyle\eta_{3}(x,y)$
$\displaystyle=3-\frac{3}{10}\sqrt{(x-47.5)^{2}+y^{2}}.$

The initial conditions are given by:

$\left\\{\begin{aligned} h&=1.0\leavevmode\nobreak\ $\mathrm{m}$&&\text{if
}x<16\leavevmode\nobreak\ $\mathrm{m}$,\>\>&h&=0.125\leavevmode\nobreak\
$\mathrm{m}$&\text{\>otherwise},\\\ \phi_{1}&=0.5&&\text{if
}x<16\leavevmode\nobreak\ $\mathrm{m}$,&\phi_{1}&=0&\text{\>otherwise},\\\
\phi_{2}&=0&&\text{if }x<16\leavevmode\nobreak\
$\mathrm{m}$,&\phi_{1}&=0.1&\text{\>otherwise},\\\
q_{x}&=q_{y}=0\leavevmode\nobreak\
$\mathrm{m}^{2}\mathrm{/}\mathrm{s}$&&\text{everywhere},\end{aligned}\right.$
(24)

presenting a virtual dam at $x$=16 m separating two discontinuous volumes of
water and scalar mass. Here the domain is discretized with 492,277 triangular
cells with a maximum of characteristic length of 0.1 m. The hydrodynamics
setup features reflective wall boundaries, a CFL of 0.95 and frictional
sources compatible with a Manning coefficient of $n$ =
$0.01\text{\,}\mathrm{m}^{-1/3}\mathrm{s}$. The scalar setup features two
initially unmixed species, both with a constant and isotropic diffusion
coefficient $K_{c}$=$0.25\text{\,}\mathrm{m}^{2}\mathrm{/}\mathrm{s}$. Fig.
17a illustrates the hydrodynamic (left) and scalar solutions (right) at the
initial condition ($t$=0 s).

At simulation start, the virtual dam collapses instantaneously, with an
advancing wave that overtops the two small lateral humps, fully circumvents
the larger hump and reaches the opposite wall in about $t$=20 s. At this time
the interface between the two species, in what would otherwise be a contact
discontinuity on flat topography, is still lagging by approximately 15 m (Fig.
17b). At this point, the reflected bores propagate upstream and further mix
both species, symmetrically around the $x$ axis. By $t$=50 s these bores
overcome the two smaller obstacles and propagate upstream on flat ground (Fig.
17c).

After a continuous sloshing and interaction of reflected waves, topography and
lateral walls, the friction sources gain relevance and dissipate most of the
kinetic energy in the flow, with a near-static solution being obtained at
approximately $t$=20 min. The scalars continue to mix, now due mostly to
molecular diffusion, in what is a much slower process, that only vanishes at
around 3.5 hours as both scalars become fully homogeneous across the domain
(Fig. 18).

The model is fully conservative, with the total liquid and scalar mass
preserved during the entire simulation. As the simulation approaches the lake-
at-rest conditions, the observed quantities correctly converge to their
resting values of $h$=0.364 m, $\phi_{1}$=0.386 and $\phi_{2}$=0.023 (Fig.
18).

$0$$15$$30$$45$$60$$75$$-10$$0$$10$$x\ (m)$$y\ (m)$

a) $t$=0 s
$0$$15$$30$$45$$60$$75$$-10$$0$$10$$0$$15$$30$$45$$60$$75$$-10$$0$$10$$x\
(m)$$y\ (m)$
b) $t$=20 s
$0$$15$$30$$45$$60$$75$$-10$$0$$10$$0$$15$$30$$45$$60$$75$$-10$$0$$10$$0$$15$$30$$45$$60$$75$$-10$$0$$10$$x\
(m)$$y\ (m)$
c) $t$=50 s
$0$$0.2$$0.4$$0.6$$0.8$$1$ $H$ (m) 0.50.250.00.050.1 $\phi_{1}$ (-) $\phi_{2}$
(-)

Figure 17: T6: dam break over complex topography. 3D visualization of the
water surface elevation $H$ (left panels) and planar view of concentration
distribution ($\phi_{1}$ and $\phi_{2}$) for both scalar species 1 and 2
(right panes). Subpanels (a,b,c) report different simulation timeout.

0.00.511.522.53$0.0$$0.1$$0.2$$0.3$$0.4$$0.5$$t$ (h)

0.00.511.522.53$0.0$$0.1$$0.2$$0.3$$0.4$$0.5$$t$ (h)

$h$ ($m$)$q_{x}$ ($m^{2}/s$)$\phi_{1}$ (-)$\phi_{2}$ (-)

Figure 18: T6: time series at two locations. Evolution of hydrodynamic
quantities ($h$ and $q_{x}$) and scalar concentrations $\phi_{1}$ (red) and
$\phi_{2}$ (blue) in time at position $(x,y)=(7.5,0.0)$ m (left) and
$(x,y)=(67.5,0.0)$ m (right). Dashed lines represent theoretical scalar
concentration values, at rest.

###  7.7 Performance and scalability

The performance and scalability of the software depends not only on the
implemented parallelization strategies but also on the physical model to be
reproduced. In general, "simpler" models, i.e. only few simulated physical
processes, are likely to show higher computational performances. To test
Basement’s computational performance, we selected the benchmarks T1
(hydrodynamic), T3 (morphodynamic) and T6 (scalar advection-diffusion).

Each of the selected numerical experiments (T1, T3, T6) has been conducted
with four different computational meshes. The sizes of these meshes are given
in Table 6 and have been chosen to cover a broad range of spatial resolutions,
ranging from thousands to hundred of thousands computational cells. The
simulations have been run with a set of computational backends. In particular,
CPU-based simulations have been performed on an Intel Xeon Gold 6154 (3.00GHz)
workstation equipped with 36 cores (two sockets with 18 physical cores each),
whilst GPU-based simulations have been run on three GPUs (GeForce GTX 1050 Ti,
GeForce GTX 1080 Ti, and Tesla P100; see Table 7 for the main
characteristics). Moreover, the simulations have been benchmarked in both
single and double precision mode. The GPUs were integrated in a workstation
with a 32-core Intel Xeon Gold 5218 (2.30GHz) processor (two sockets with 16
physical cores each).

Table 6: Number of computational cells for the performance and scalability benchmarks. Mesh ID | T1 | T3 | T6
---|---|---|---
1 | $24\,945\text{\,}$ | $27\,444\text{\,}$ | $24\,388\text{\,}$
2 | $52\,102\text{\,}$ | $47\,187\text{\,}$ | $49\,155\text{\,}$
3 | $101\,417\text{\,}$ | $109\,344\text{\,}$ | $98\,163\text{\,}$
4 | $499\,060\text{\,}$ | $218\,912\text{\,}$ | $196\,829\text{\,}$
Table 7: Characteristics of the benchmarked Nvidia GPU cards. Further specifications available at www.nvidia.com/en-gb/geforce/10-series/ and www.nvidia.com/en-gb/data-center/tesla-p100/. Type | Architecture | CUDA Cores
---|---|---
GeForce GTX 1050 Ti | Pascal | 768
GeForce GTX 1080 Ti | Pascal | 3584
Tesla P100 | Pascal | 3584

For a given a mesh size, the speedup achieved by a parallelized backend $p$ is
computed using the formula $\textrm{speedup}=T_{s}/T_{p}$, where $T_{s}$
($T_{p}$) is the total computational time used by the serial (parallel)
backend. The results for all the investigated benchmarks are depicted in Fig.
19.

(a) T1

(b) T3

(c) T6

Figure 19: Speedup of backends for varying mesh size, test cases T1, T3 and
T6. The benchmarks were executed with different degrees of parallelism on the
CPU (using OpenMP) and the GPU. The final "S" and "D" in GPU series denote
single and double precision, respectively.

As anticipated, the speedup depends on the simulated processes. Comparing the
speedup values among different benchmarks in Fig. 19, cases T1 (hydrodynamics)
and T3 (morphodynamics) show higher values on average than T6 (advection-
diffusion). Such results are expected, given the increased complexity (number
of equations and operations to be solved) of T6.

The performance benefits of Basement’s parallelization can be evaluated in
more detail by comparing the speedup values along the vertical axis of the
plots. In the following we focus on benchmark T3 (morphodynamics) which shows
an "intermediate" scalability among the three benchmarks (Fig. 19b). Looking
at the mesh with 47k elements as an example, the CPU-based family (i.e. OpenMP
on multiple cores) shows a speedup efficiency (i.e.
$\textrm{speedup}/\textrm{ncores}\cdot 100$) of 87% with 2 CPU cores
(Speedup=1.7), and of 58% with 32 cores (Speedup=18.5), with an average
efficiency of 74%. Basement performs even better on some GPU cards. The least
performing card (GTX 1050 Ti with double precision) has a speedup of about 8.
However, note that speedup jumps to 20 when using the single precision
version. Overall, the speedup provided by the tested GPUs ranges between 7 and
60.

The benchmarks in Fig. 19 also show how the maximum speedup changes with mesh
size, computational backend and simulated processes. For all three cases, the
CPU-based parallelizations show a mild speedup increase with an increasing
number of computational cells. The dependency on the mesh size is slightly
more pronounced when the number of computational cores is increased. This
reflects the fact that CPU-based solutions have shared memory and minimal
overhead (for multi-threading handling), thus the domain size (i.e. the data
size) does not represent a potential performance bottleneck. Focussing on T3
(Fig. 19b), the parallelization efficiency for 2-4 CPU cores is almost
constant for all mesh sizes and above 80%. On the other hand, the efficiency
for 16-32 CPU cores is larger than 70% only for mesh 4 (218k).

The speedup of the GPU-accelerated solutions shows not only a more marked
dependency on the problem size, but also on the simulated processes. In
benchmark T1 (Fig. 19a), the speedup clearly increases with problem size. This
can be explained with the overhead of GPU parallelization, which becomes more
and more negligible with increasing domain size. Conversely, benchmark T6
(Fig. 19c) shows little impact of the domain size on the speedup. In this case
the scalability is limited by data transfers (i.e. data bandwidth) given the
larger amount of data (compared to test case T1) needed for this simulation.

It is worth remarking that the GPU-accelerated backends show an average
speedup difference greater than 10 between the single and double precision
versions. Of course, the adequate choice depends on the requirements of the
specific application.

The analysis above shows Basement’s performance on different computational
backends and underlines the differences when simulating different processes.
The results summarized in Fig. 19 can also serve as a guideline for the
interested reader/user when choosing an appropriate computational
configuration for a given application. Finally it is worth highlighting that
all the tested hardware configurations can be easily installed in standard
office workstations.

##  8 Conclusions

In this paper we introduced the main features of Basement version 3, a
freeware tool for river simulation. Basement allows the simulation of a wide
variety of hydro-, morphodynamic, and scalar advection-diffusion scenarios. As
illustrated with the test cases, the software is able to efficiently capture
large scale hydrodynamic processes modelled with several hundreds of thousand
elements in good agreement with the measurements. On the opposite end of the
spectrum, the morphological solver is able to handle demanding sediment
transport scenarios well, albeit with known limitations. With the scalar
advection-diffusion module a further set of physical processes such as the
fate of river pollutants can be accurately modelled.

The impact of this flexibility on the software performance is minimized by
activating feature sets on request in Basement’s pre-simulation step. The
advantage of this approach is that only the required tasks are scheduled for
execution. This, together with OP2’s ability to generate executable code for
both multi-core CPU’s and GPUs, permit Basement to scale with both available
features and available computational power. Such advantages are reflected in
the presented benchmarks. Given a large enough domain, the software shows a
good parallel efficiency on the CPU and an even higher speedup when using
GPUs.

The Basement project is in continuous advancement to optimize and include
further features in the existing basic modules. As an example, the modelling
of the sediment transport in presence of non uniform sediment size and the
simulation of water temperature dynamics are in implementation phase. On the
other hand, efforts are dedicated also to develop novel modelling solutions
for river processes such as the bio-morphodynamic feedbacks between vegetation
and sediment transport. Table 2 of the Supplementary Material provides an
overview of under development features. The modularity of the development
framework allows also for further refinement of single specific numerical
solvers and the implementation of high-order schemes when needed.

Overall, the combination of different river processes that can be modelled,
the computational efficiency, the flexibility in the backend choice, but also
the availability of a light Graphical User Interface, make Basement a valuable
tool for a broad family of river modelers in both Academia and Practice.

## Acknowledgments

The Authors greatly thank the many former collaborators and developers within
the Basement project. Particular thanks to Aurélie Koch, for her valuable
contribution in testing and documenting the software.

The design of Basement was conducted by DFV, DV, SP, LV. The software
prototyping, development and implementation was done by DV, SP, LV, MB, MW,
with the coordination and supervision of DFV and AS. Implementation and
testing of the reported features was done by MB, DV, MW and DC. The manuscript
was conceptualized by DV, AS and DFV, and drafted by DV, with support of MB
and MW. All Authors contributed to the manuscript review.

## Conflicts of interest

The Authors declare that there are no conflicts of interest.

## Funding

The development of the software Basement is financially supported by the Swiss
Federal Office for the Environment (BAFU).

## References

## References

  * Gilvear et al. [2016] D. J. Gilvear, M. T. Greenwood, M. C. Thoms, P. J. Wood (Eds.), River Science, John Wiley & Sons, Ltd, doi:10.1002/9781118643525, 2016.
  * Brewer et al. [2018] S. K. Brewer, T. A. Worthington, R. Mollenhauer, D. R. Stewart, R. A. McManamay, L. Guertault, D. Moore, Synthesizing models useful for ecohydrology and ecohydraulic approaches: An emphasis on integrating models to address complex research questions, Ecohydrology 11 (7) (2018) e1966, ISSN 19360584, doi:10.1002/eco.1966.
  * Shimizu et al. [2019] Y. Shimizu, J. Nelson, K. A. Ferrel, K. Asahi, S. Giri, T. Inoue, T. Iwasaki, C.-L. Jang, T. Kang, I. Kimura, T. Kyuka, J. Mishra, M. Nabi, S. Patsinghasanee, S. Yamaguchi, Advances in computational morphodynamics using the International River Interface Cooperative (iRIC) software, Earth Surface Processes and Landforms 45 (1) (2019) 11–37, doi:10.1002/esp.4653.
  * Crosato and Saleh [2011] A. Crosato, M. S. Saleh, Numerical study on the effects of floodplain vegetation on river planform style, Earth Surface Processes and Landforms 36 (6) (2011) 711–720.
  * Sharma and Kansal [2012] D. Sharma, A. Kansal, Assessment of river quality models: a review, Reviews in Environmental Science and Bio/Technology 12 (3) (2012) 285–311, doi:10.1007/s11157-012-9285-8.
  * Williams et al. [2016] R. D. Williams, J. Brasington, D. M. Hicks, Numerical Modelling of Braided River Morphodynamics: Review and Future Challenges, Geography Compass 10 (3) (2016) 102–127, doi:10.1111/gec3.12260.
  * Dugdale et al. [2017] S. J. Dugdale, D. M. Hannah, I. A. Malcolm, River temperature modelling: A review of process-based approaches and future directions, Earth-Science Reviews 175 (2017) 97–113, doi:10.1016/j.earscirev.2017.10.009.
  * Teng et al. [2017] J. Teng, A. J. Jakeman, J. Vaze, B. F. W. Croke, D. Dutta, S. Kim, Flood inundation modelling: A review of methods, recent advances and uncertainty analysis, Environmental Modelling and Software 90 (2017) 201–216, ISSN 13648152, doi:10.1016/j.envsoft.2017.01.006.
  * Zischg et al. [2018] A. P. Zischg, M. Mosimann, D. B. Bernet, V. Röthlisberger, Validation of 2D flood models with insurance claims, Journal of hydrology 557 (2018) 350–361.
  * Marcus and Fonstad [2010] W. A. Marcus, M. A. Fonstad, Remote sensing of rivers: the emergence of a subdiscipline in the river sciences, Earth Surface Processes and Landforms 35 (15) (2010) 1867–1872, doi:10.1002/esp.2094.
  * Thomas Steven Savage et al. [2016] J. Thomas Steven Savage, F. Pianosi, P. Bates, J. Freer, T. Wagener, Quantifying the importance of spatial resolution and other factors through global sensitivity analysis of a flood inundation model, Water Resources Research 52 (11) (2016) 9146–9163.
  * Pasternack [2011] G. Pasternack, 2D modeling and ecohydraulic analysis, University of California at Davis, California, ISBN 9781466320093, 2011\.
  * Maddock [1999] I. Maddock, The importance of physical habitat assessment for evaluating river health, Freshwater Biology 41 (2) (1999) 373–391, doi:10.1046/j.1365-2427.1999.00437.x.
  * Siviglia et al. [2013] A. Siviglia, G. Stecca, D. Vanzo, G. Zolezzi, E. F. Toro, M. Tubino, Numerical modelling of two-dimensional morphodynamics with applications to river bars and bifurcations, Advances in Water Resources 52 (2013) 243–260, ISSN 03091708, doi:10.1016/j.advwatres.2012.11.010.
  * Wyrick et al. [2014] J. Wyrick, A. Senter, G. Pasternack, Revealing the natural complexity of fluvial morphology through 2D hydrodynamic delineation of river landforms, Geomorphology 210 (2014) 14–22, doi:10.1016/j.geomorph.2013.12.013.
  * Guan et al. [2016] M. Guan, N. Wright, P. Sleigh, S. Ahilan, R. Lamb, Physical complexity to model morphological changes at a natural channel bend, Water resources research 52 (8) (2016) 6348–6364, doi:10.1002/2015WR017917.
  * Siviglia and Crosato [2016] A. Siviglia, A. Crosato, Numerical modelling of river morphodynamics: Latest developments and remaining challenges, Advances in Water Resources 93 (2016) 1–3, doi:10.1016/j.advwatres.2016.01.005.
  * Toro [2001] E. Toro, Shock-capturing methods for free-surface shallow flows, Wiley and Sons Ltd, ISBN 0471987662, doi:10.1080/00221680309499935, 2001\.
  * Nahorniak et al. [2018] M. Nahorniak, J. Wheaton, C. Volk, P. Bailey, M. Reimer, E. Wall, K. Whitehead, C. Jordan, How do we efficiently generate high-resolution hydraulic models at large numbers of riverine reaches?, Computers & Geosciences 119 (November 2017) (2018) 80–91, ISSN 00983004, doi:10.1016/j.cageo.2018.07.001.
  * Costabile and Macchione [2015] P. Costabile, F. Macchione, Enhancing river model set-up for 2-D dynamic flood modelling, Environmental Modelling & Software 67 (2015) 89–107, ISSN 13648152, doi:10.1016/j.envsoft.2015.01.009.
  * Dazzi et al. [2020] S. Dazzi, R. Vacondio, P. Mignosa, Internal boundary conditions for a GPU-accelerated 2D shallow water model: Implementation and applications, Advances in Water Resources 137 (2020) 103525.
  * Sanders [2008] B. F. Sanders, Integration of a shallow water model with a local time step, Journal of Hydraulic Research 46 (4) (2008) 466–475, doi:10.3826/jhr.2008.3243.
  * Dazzi et al. [2018] S. Dazzi, R. Vacondio, A. D. Palù, P. Mignosa, A local time stepping algorithm for GPU-accelerated 2D shallow water models, Advances in Water Resources 111 (2018) 274–288, doi:10.1016/j.advwatres.2017.11.023.
  * Powell et al. [1993] K. G. Powell, P. L. Roe, J. Quirk, Adaptive-Mesh Algorithms for Computational Fluid Dynamics, in: Algorithmic Trends in Computational Fluid Dynamics, Springer New York, 303–337, doi:10.1007/978-1-4612-2708-3_18, 1993\.
  * Carraro et al. [2018] F. Carraro, D. Vanzo, V. Caleffi, A. Valiani, A. Siviglia, Mathematical study of linear morphodynamic acceleration and derivation of the MASSPEED approach, Advances in Water Resources 117 (2018) 40–52.
  * Morgan et al. [2020] J. A. Morgan, N. Kumar, A. R. Horner-Devine, S. Ahrendt, E. Istanbullouglu, C. Bandaragoda, The use of a morphological acceleration factor in the simulation of large-scale fluvial morphodynamics, Geomorphology 356 (2020) 107088, doi:10.1016/j.geomorph.2020.107088.
  * Vanzo et al. [2016] D. Vanzo, A. Siviglia, E. F. Toro, Pollutant transport by shallow water equations on unstructured meshes: Hyperbolization of the model and numerical solution via a novel flux splitting scheme, Journal of Computational Physics 321 (2016) 1–20.
  * Afzal et al. [2016] A. Afzal, Z. Ansari, A. R. Faizabadi, M. K. Ramis, Parallelization Strategies for Computational Fluid Dynamics Software: State of the Art Review, Archives of Computational Methods in Engineering 24 (2) (2016) 337–363, doi:10.1007/s11831-016-9165-4.
  * Owens et al. [2008] J. Owens, M. Houston, D. Luebke, S. Green, J. Stone, J. Phillips, GPU Computing, Proceedings of the IEEE 96 (5) (2008) 879–899, doi:10.1109/jproc.2008.917757.
  * Mudalige et al. [2012] G. Mudalige, M. Giles, I. Reguly, C. Bertolli, P. Kelly, OP2: An active library framework for solving unstructured mesh-based applications on multi-core and many-core architectures, in: 2012 Innovative Parallel Computing (InPar), IEEE, ISBN 978-1-4673-2633-9, 1–12, doi:10.1109/InPar.2012.6339594, 2012\.
  * Brodtkorb et al. [2012] A. R. Brodtkorb, M. L. Sætra, M. Altinakar, Efficient shallow water simulations on GPUs: Implementation, visualization, verification, and validation, Computers and Fluids 55 (2012) 1–12, ISSN 00457930, doi:10.1016/j.compfluid.2011.10.012.
  * Smith and Liang [2013] L. S. Smith, Q. Liang, Towards a generalised GPU/CPU shallow-flow modelling tool, Computers & Fluids 88 (2013) 334–343, ISSN 00457930, doi:10.1016/j.compfluid.2013.09.018.
  * Vacondio et al. [2014] R. Vacondio, A. D. Palù, P. Mignosa, GPU-enhanced Finite Volume Shallow Water solver for fast flood simulations, Environmental Modelling & Software 57 (2014) 60–75, ISSN 1364-8152, doi:10.1016/j.envsoft.2014.02.003.
  * Vacondio et al. [2017] R. Vacondio, A. Dal Palù, A. Ferrari, P. Mignosa, F. Aureli, S. Dazzi, A non-uniform efficient grid type for GPU-parallel Shallow Water Equations models, Environmental Modelling & Software 88 (2017) 119–137, ISSN 13648152, doi:10.1016/j.envsoft.2016.11.012.
  * Hou et al. [2020] J. Hou, Y. Kang, C. Hu, Y. Tong, B. Pan, J. Xia, A GPU-based numerical model coupling hydrodynamical and morphological processes, International Journal of Sediment Research 35 (4) (2020) 386–394, doi:10.1016/j.ijsrc.2020.02.005.
  * Castro et al. [2011] M. J. Castro, S. Ortega, M. de la Asunción, J. M. Mantas, J. M. Gallardo, GPU computing for shallow water flow simulation based on finite volume schemes, Comptes Rendus Mécanique 339 (2-3) (2011) 165–184, ISSN 16310721, doi:10.1016/j.crme.2010.12.004.
  * Lacasta et al. [2014] A. Lacasta, M. Morales-Hernández, J. Murillo, P. García-Navarro, An optimized GPU implementation of a 2D free surface simulation model on unstructured meshes, Advances in Engineering Software 78 (2014) 1–15, ISSN 09659978, doi:10.1016/j.advengsoft.2014.08.007.
  * Lacasta et al. [2015] A. Lacasta, C. Juez, J. Murillo, P. García-Navarro, An efficient solution for hazardous geophysical flows simulation using GPUs, Computers & Geosciences 78 (2015) 63–72, doi:10.1016/j.cageo.2015.02.010.
  * Petaccia et al. [2016] G. Petaccia, F. Leporati, E. Torti, OpenMP and CUDA simulations of Sella Zerbino Dam break on unstructured grids, Computational Geosciences 20 (5) (2016) 1123–1132, doi:10.1007/s10596-016-9580-5.
  * Juez et al. [2016] C. Juez, A. Lacasta, J. Murillo, P. García-Navarro, An efficient GPU implementation for a faster simulation of unsteady bed-load transport, Journal of Hydraulic Research 54 (3) (2016) 275–288, ISSN 0022-1686, doi:10.1080/00221686.2016.1143042.
  * García-Feal et al. [2018] O. García-Feal, J. González-Cao, M. Gómez-Gesteira, L. Cea, J. Domínguez, A. Formella, An Accelerated Tool for Flood Modelling Based on Iber, Water 10 (10) (2018) 1459, doi:10.3390/w10101459.
  * Reguly et al. [2016] I. Z. Reguly, G. R. Mudalige, C. Bertolli, M. B. Giles, A. Betts, P. H. Kelly, D. Radford, Acceleration of a Full-Scale Industrial CFD Application with OP2, IEEE Transactions on Parallel and Distributed Systems 27 (5) (2016) 1265–1278, ISSN 1045-9219, doi:10.1109/TPDS.2015.2453972.
  * Giles et al. [2012] M. B. Giles, G. R. Mudalige, Z. Sharif, G. Markall, P. H. J. Kelly, Performance Analysis and Optimization of the OP2 Framework on Many-Core Architectures, The Computer Journal 55 (2) (2012) 168–180, ISSN 0010-4620, doi:10.1093/comjnl/bxr062.
  * Reguly et al. [2018] I. Z. Reguly, D. Giles, D. Gopinathan, L. Quivy, J. H. Beck, M. B. Giles, S. Guillas, F. Dias, The VOLNA-OP2 tsunami code (version 1.5), Geoscientific Model Development 11 (11) (2018) 4621–4635.
  * Giles et al. [2020] D. Giles, E. Kashdan, D. M. Salmanidou, S. Guillas, F. Dias, Performance analysis of Volna-OP2–massively parallel code for tsunami modelling, Computers & Fluids 209 (2020) 104649.
  * Beckers et al. [2020] F. Beckers, A. Heredia, M. Noack, W. Nowak, S. Wieprecht, S. Oladyshkin, Bayesian Calibration and Validation of a Large-Scale and Time-Demanding Sediment Transport Model, Water Resources Research 56 (7) (2020) e2019WR026966.
  * Jung and Merwade [2015] Y. Jung, V. Merwade, Estimation of uncertainty propagation in flood inundation mapping using a 1-D hydraulic model, Hydrological Processes 29 (4) (2015) 624–640.
  * Peter [2017] S. J. Peter, Dam Break Analysis under Uncertainty, Ph.D. thesis, ETH Zurich, Zurich, Switzerland, doi:10.3929/ethz-b-000209879, 2017.
  * Armanini [2018] A. Armanini, Principles of river hydraulics, Springer, 2018.
  * Graf [1966] W. H. Graf, On the determination of the roughness coefficient in natural and artificial waterways, Hydrological Sciences Journal 11 (1) (1966) 59–68.
  * Bezzola [2002] G. R. Bezzola, Fliesswiederstand und Sohlenstabilität natürlicher Gerinne, Ph.D. thesis, Eidgenössische Technische Hochschule Zürich, 2002\.
  * Exner [1925] F. M. Exner, Ueber die Wechselwirkung zwischen Wasser und Geschiebe in Fluessen, Tech. Rep., Akademie der Wissenschaften, Mathematische Naturwissenschaft Abt. IIa, Wien, Austria, 1925.
  * Parker [2004] G. Parker, 1D sediment transport morphodynamics with applications to rivers and turbidity currents: E-book, Minneapolis, MN, available from: http://hydrolab.illinois.edu/people/parkerg/morphodynamics_e-book.htm .
  * Engelund and Hansen [1972] F. Engelund, E. Hansen, A monograph on sediment transport in alluvial streams, Teknisk Forlag, Copenhagen, 1972.
  * Meyer-Peter and Müller [1948] E. Meyer-Peter, R. Müller, Formulas for bed-load transport, in: IAHSR 2nd meeting, Stockholm, appendix 2, IAHR, –, 1948.
  * Wong and Parker [2006] M. Wong, G. Parker, Reanalysis and correction of bed-load relation of Meyer-Peter and Müller using their own database, Journal of Hydraulic Engineering 132 (11) (2006) 1159–1168.
  * Grass [1981] A. J. Grass, Sediment transport by waves and currents, University College, London, Dept. of Civil Engineering, 1981.
  * Smart and Jaeggi [1983] G. M. Smart, M. N. R. Jaeggi, Sediment Transport on Steep Slopes, VAW-Mitteilung 64, Versuchsanstalt für Wasserbau,Hydrologie und Glaziologie (VAW). Zürich, ETH Zürich., 1983.
  * Shields [1936] A. Shields, Anwendungen der Ähnlichkeitsmechanik und der Turbulenzforschung auf die Geschiebebewegungen, Tech. Rep., Mitteilung der Preussischen Versuchsanstalt für Wasserbau und Schiffbau. Berlin, Deutschland, 1936.
  * van Rijn [1989] L. C. van Rijn, Handbook Sediment Transport by Current and Waves, Delft Hydraulics Laboratory, Delft, The Netherlands, 1989.
  * Chen et al. [2010] X. Chen, J. Ma, S. Dey, Sediment Transport on Arbitrary Slopes: Simplified Model, Journal of Hydraulic Engineering 136 (5) (2010) 311–317, ISSN 0733-9429, doi:10.1061/(ASCE)HY.1943-7900.0000175.
  * Ikeda [1982] S. Ikeda, Lateral bed-load transport on side slopes, Journal of the Hydraulics Division, ASCE 108 (11) (1982) 1369–1373.
  * Talmon et al. [1995] A. M. Talmon, N. Struiksma, M. Van Mierlo, Laboratory measurements of the direction of sediment transport on transverse alluvial-bed slopes, Journal of Hydraulic Research 33 (4) (1995) 495–517, ISSN 0022-1686, doi:10.1080/00221689509498657.
  * Engelund [1974] F. Engelund, Flow and bed topography in channel bends., Journal of the Hydraulics Division ASCE 100 (11) (1974) 1631–1648.
  * Rozovskii [1961] I. L. Rozovskii, Flow of Water in Bends of Open Channels, Academy of Science of the Ukrainian S.S.R, Institute of Hydrology and Hydraulic Engineering, 1961\.
  * Toro et al. [1994] E. F. Toro, M. Spruce, W. Speares, Restoration of the contact surface in the HLL-Riemann solver, Shock waves 4 (1) (1994) 25–34.
  * Toro [2009] E. F. Toro, Riemann solvers and numerical methods for fluid dynamics., Springer-Verlag GmbH, ISBN 978-3-540-49834-6, doi:10.1007/978-3-540-49834-6, 2009.
  * Vanzo [2015] D. Vanzo, Eco-hydraulic quantification of hydropeaking and thermopeaking: development of modeling and assessment tools, Ph.D. thesis, University of Trento, Trento, Italy, 2015.
  * Duran et al. [2013] A. Duran, Q. Liang, F. Marche, On the well-balanced numerical discretization of shallow water equations on unstructured meshes, Journal of Computational Physics 235 (2013) 565–586, ISSN 00219991, doi:10.1016/j.jcp.2012.10.033.
  * Soares-Frazão and Zech [2011] S. Soares-Frazão, Y. Zech, HLLC scheme with novel wave-speed estimators appropriate for two-dimensional shallow-water flow on erodible bed, International Journal for Numerical Methods in Fluids 66 (8) (2011) 1019–1036, ISSN 02712091, doi:10.1002/fld.2300.
  * Hervouet and Petitjean [1999] J.-M. Hervouet, A. Petitjean, Malpasset dam-break revisited with two-dimensional computations, Journal of hydraulic research 37 (6) (1999) 777–788.
  * Valiani et al. [2002] A. Valiani, V. Caleffi, A. Zanni, Case study: Malpasset dam-break simulation using a two-dimensional finite volume method, Journal of Hydraulic Engineering 128 (5) (2002) 460–472.
  * Singh et al. [2011] J. Singh, M. S. Altinakar, Y. Ding, Two-dimensional numerical modeling of dam-break flows over natural terrain using a central explicit scheme, Advances in Water Resources 34 (10) (2011) 1366–1375.
  * Cordier et al. [2011] S. Cordier, M. H. Le, T. M. De Luna, Bedload transport in shallow water models: Why splitting (may) fail, how hyperbolicity (can) help, Advances in Water Resources 34 (8) (2011) 980–989.
  * Bellal et al. [2003] M. Bellal, B. Spinewine, C. Savary, Y. Zech, Morphological evolution of steep-sloped river beds in the presence of a hydraulic jump: Experimental study, in: XXX IAHR Congress, Citeseer, 133–140, 2003.
  * Goutiere et al. [2011] L. Goutiere, S. Soares-Frazão, Y. Zech, Dam-break flow on mobile bed in abruptly widening channel: experimental data, Journal of Hydraulic Research 49 (3) (2011) 367–371, doi:10.1080/00221686.2010.548969.
  * Juez et al. [2014] C. Juez, J. Murillo, P. García-Navarro, A 2D weakly-coupled and efficient numerical model for transient shallow flow and movable bed, Advances in Water Resources 71 (2014) 93–109.
  * Xia et al. [2010] J. Xia, B. Lin, R. A. Falconer, G. Wang, Modelling dam-break flows over mobile beds using a 2D coupled approach, Advances in Water Resources 33 (2) (2010) 171–183, doi:10.1016/j.advwatres.2009.11.004.
  * Yen and Lee [1995] C.-l. Yen, K. T. Lee, Bed topography and sediment sorting in channel bend with unsteady flow, Journal of Hydraulic Engineering 121 (8) (1995) 591–599.
  * Kaveh et al. [2019] K. Kaveh, M. Reisenbüchler, S. Lamichhane, T. Liepert, N. D. Nguyen, M. D. Bui, P. Rutschmann, A Comparative Study of Comprehensive Modeling Systems for Sediment Transport in a Curved Open Channel, Water 11 (9) (2019) 1779.
  * Brufau et al. [2002] P. Brufau, M. E. Vázquez-Cendón, P. García-Navarro, A numerical model for the flooding and drying of irregular domains, International Journal for Numerical Methods in Fluids 39 (3) (2002) 247–275, doi:10.1002/fld.285.