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
BBasement v3: a modular freeware for river process modelling over multiplecomputational backends
Davide Vanzo a,b, ∗ , Samuel Peter e , Lukas Vonwiller f , Matthias Buergler a , Manuel Weberndorfer d ,Annunziato Siviglia c , Daniel Conde a , David F. Vetsch a a Laboratory of Hydraulics, Hydrology and Glaciology. ETH, Swiss Federal Institute of Technology, Zürich, Switzerland. b Dept. of Surface Waters - Research and Management. Eawag, Swiss Federal Institute of Aquatic Science and Technology,Kastanienbaum, Switzerland. c Dept. Civil, Environmental and Mechanical Engineering. University of Trento, Trento, Italy. d ID Scientific IT Services. ETH, Swiss Federal Institute of Technology, Zürich, Switzerland. e Axpo Group, Baden, Switzerland. f TK CONSULT AG , Zürich, Switzerland.
Abstract
Modelling river physical processes is of critical importance for flood protection, river management andrestoration of riverine environments. Developments in algorithms and computational power have led to awider spread of river simulation tools. However, the use of two-dimensional models can still be hindered bycomplexity 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, morphodynamicand scalar advection-diffusion processes.
Basement leverages different computational platforms (multi-coreCPUs and graphics processing units GPUs) to enable the simulation of large domains and long-term riverprocesses. The adoption of a fully costless worflow and a light GUI facilitate its broad utilization. We testits robustness and efficiency in a selection of benchmarks. Results confirm that
Basement could be anefficient and versatile tool for research, engineering practice and education in river modelling.
Keywords:
GPU-CUDA, river modelling, unstructured grid, shallow water, sediment transport, pollutanttransport
Software availabilityName of software:
Basement , version 3 (v3)
Website:
E-mail: [email protected]
Developer:
Numerical Modelling division at the Laboratory of Hydraulics, Hydrology and Glaziology(VAW), ETH Zürich ∗ corresponding author: [email protected] Preprint submitted to Environmental Modelling & Software February 26, 2021 a r X i v : . [ phy s i c s . f l u - dyn ] F e b anguage: 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 environmentalsciences. River science [ sensu
1] is no exception in this trend, with a number of tools been developed toaddress 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 thelast decades allowed for the spread of more and more sophisticated numerical tools. In the context of riverscience, two-dimensional depth-averaged (hereinafter 2D) models are nowadays of common use in researchand engineering practice. This is particularly true for some applications such as flood modelling and rivermorphodynamics [e.g. 3, 9]. The increasing usage of 2D river models is also closely bonded with the growingavailability of high-resolution river datasets. In particular, advances in LiDAR, UAV-Photography andothers 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 ad-hoc implementations of GPU-based models for 2D hydrodynamic [e.g. 31, 32, 33, 34], and occasionally mor-phodynamic simulations [e.g. 35] have become available. The vast majority of these models are based onstructured 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 theneed for mappings and indirections. To the best of the Authors’ knowledge, few hydrodynamic models im-plement GPU-acceleration on unstructured meshes [36, 37, 38, 39], with very limited ad-hoc implementationsfor 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]. Anincrease in availability of freeware GPU-based river models would be beneficial for environmental modelersin academic research and education, but also in consultancy and engineering offices.3n this paper we introduce the
Basement software (version 3), a freeware application developed atthe 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 prac-tical interest. It can seamlessly run on GPU-enabled workstations, as well as on more standard multi-coreCPUs. This flexibility in the choice of the computational backend is achieved by integrating the OP2framework [30, 42, 43]. This framework provides an additional abstract layer for the acceleration of nu-merical models on unstructured computational meshes, and has been successfully implemented in similarmodelling context [44]. The obtained parallelization performance alleviates the computational limitationswhen 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] ofdeterministic 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 alight Graphical User Interface (GUI). The
Basement software aims to enable a broad range of potentialusers 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 documentationmaterial and tutorials, the software is designed to be a valuable didactic tool for engineering and river sciencestudents.The paper is structured as follows: §2 provides the software application context that justifies the adoptedmathematical and numerical strategies. §3 to §5 report the mathematical basis, the numerical strategiesand main features of the basic modules of
Basement . The software design, the modelling workflow and theparallelization solutions are presented in §6. A selection of benchmarks are reported in §7, whilst conclusionsand 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 atdifferent spatial and temporal scales. For example,
Basement can be used to simulate large scale (i.e. basinscale) flood propagation, but also reach scale morphodynamic processes such as formation and evolutionof fluvial bars. Moreover, it can be applied together with high-resolution topographies (in the order ofcentimeters) to simulate ecohydraulic processes at different ecological scales (e.g. habitat modelling). Thisrange 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 river4ows in Alpine contexts, the propagation of natural flood waves as well as hydropeaking events. This isensured by the adoption of a robust and accurate shock-capturing explicit solver for the hydrodynamicproblem (§4.2); accurate front propagation: it is possible to simulate extreme events such as dam-break induced floods, butalso ecologically-relevant processes such as the wetting-drying of riparian areas and in-channel mor-phologies due to artificial flow alterations. This is achieved by an implemented shock-wave capturingnumerical 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 enablesfor an accurate description of complex river morphologies and riverrine structures (§4.1). The adoptionof an unstructured mesh also reduces the strong anisotropy of structured meshes, which can be crucialfor particular applications. large problems: the software adopts a parallelization strategy tailored to the acceleration of problems onunstructured meshes (§6.4). Moreover,
Basement simulations can efficiently be executed on differentcomputational backends. Those backends inlude GPU cards, therefore allowing for the simulation oflarge 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 ashydrodynamics, sediment or advection-diffusion of a scalar (e.g. a non-reactive pollutant) can besimulated 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 goodparallelization 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 describingthe 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 aso-called Initial-Boundary Value Problem [18], where process-specific initial and boundary conditions arerequired 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 basicmodules. The main module features are also listed in Table 1 of the Supplementary Material.5 . Mathematical formulation
The hydrodynamic module solves the so-called shallow water equations (hereinafter SWE) [e.g. 18]. Thetwo-dimensional SWE are of practical interest with regard to water flows with a free surface under theinfluence of gravity.Considering a Cartesian reference system ( x, y, z ) where the z axis is vertical and the ( x, y ) plane ishorizontal (Fig. 1a), the system of governing equations can be expressed as: ∂ t H + ∂ x q x + ∂ y q y = S h ∂ t q x + ∂ x (cid:18) q x h + 12 gH − gHz B (cid:19) + ∂ y (cid:16) q x q y h (cid:17) = − gH∂ x z B − ghS fx ∂ t q y + ∂ x (cid:16) q x q y h (cid:17) + ∂ y (cid:32) q y h + 12 gH − gHz B (cid:33) = − gH∂ y z B − ghS fy (1)where the system unknowns are the water surface elevation H [m], and the two directional components of q = ( q x , q y ) [ m / s ], representing the flow discharge per unit width. With z B [m] we indicate the bottomelevation, whilst h = ( H − z B ) [m] is the water depth, and g [ m / s ] the acceleration due to gravity. Notethat the depth-averaged velocity vector can be consequently expressed as u = ( u, v ) = ( q x /h, q y /h ) [ m / s ].Finally S fx and S fy [-] represent the dimensionless friction terms in x and y direction, whilst S h [ m / s ]represents potential external contribution/subtraction of flow discharge to the mass conservation equation. To solve the system (1), closure relationships for the friction terms S fx , S fy and the contribution ofexternal inflow/outflow discharge S h must be provided. Friction terms.
Under the hypothesis of turbulent flow, hence under the assumption that the energy lineslope is proportional to the square of the flow velocity, the friction terms S fx , S fy can be written as: S fx = u (cid:107) u (cid:107) ghc f ; S fy = v (cid:107) u (cid:107) ghc f , (2)where c f is the dimensionless friction coefficient. Several formulae are available for c f . Basement imple-ments four well know formulations of power or logarithmic type, given in Table 1.
External inflow/outflow discharge.
The term S h [ m / s ] is used to represent additional sources of water likerainfall and springs, or water abstraction (sink), and can be defined over subsets of the computationaldomain. The external water source can be provided by the user as total discharge [ m / s ] or as intensity[ mm / h ], per squared meter. Different behaviour can be imposed for each external source/sink area:6 xy H hz B ϕ c ϕ b qq B sn 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 ( q B ) deviation angle ϕ b fromthe flow direction q due to gravitational effects caused by the local lateral slope s . (c) Bed load transport ( q B ) deviation angle ϕ c from the flow direction q due to the spiral flow motion caused by the curvature of radius Rc .Table 1: Friction closure relationships for the hydrodynamic problem . Formulations for the dimensionless frictioncoefficient c f ; for both Chézy and Bezzola entries, d is the 90 th percentile of the sediment grain size distribution. Closure Expression Parameters Range Ref
Strickler c f = k str h / / √ g k str [ m / s − ] 7-40 [49]Manning c f = h / / ( n √ g ) n [ m − / s ] 0.025-0.143 [49]Chézy c f = 6 .
25 + log ( h/K s ) for h > K s c f = 6 . for h ≤ K s K s [m] K s = n k d , [50]with n k = 2 ÷ Bezzola c f = 2 . (cid:112) − y R /h ln (10 . h/y R ) , for h/y R > c f = 1 . (cid:112) h/y R ln (10 . h/y R ) , for . ≤ h/y R ≤ c f = 1 . , for h/y R < . y R [m] y R = nd , n ≈ [51] The basic morphodynamic module solves the so-called Exner equation [52]. It describes the bed evolutiondue to erosion or deposition, which results in the elevation change of the actual bed level z B . Assuming thesame coordinate reference of Fig. 1a, it reads (1 − p ) ∂ t z B + ∂ x q B x + ∂ y q B y = S b , (3)where p [-] is the bed sediment porosity, assumed constant in space and time, S b [ m / s ] is an externalsource term specifying local inputs or outputs of sediment material (e.g. slope collapse or excavation) and q B = ( q B x , q B y ) [ m / s ] is the specific sediment transport flux.The morphodynamic module in its basic form accounts only for sediment transport occurring in the formof bed-load or total-load [49, 53]. The simulation of sediment transport as suspended load is delegated to aspecific module, planned in future versions of the software. Two closure relationships are needed to numerically solve the problem of (3): a sediment transportformula and the external source/sink of sediments.
Sediment transport formula.
Basement implements three different types of sediment transport formulaeas 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 andHansen 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 (cid:107) q B (cid:107) are given in Table 2. In the MPM-like formulation, θ is the dimensionless bed shear stress (i.e. Shields8 able 2: Sediment transport closure relationships for the morphodynamic problem.
Expressions provide an esti-mation of the specific sediment transport magnitude (cid:107) q B (cid:107) [ m / s ]. Type Expression Parameters Ref
MPM-like α ( θ − θ cr ) m (cid:112) ( s − gd m α = 8 , m = 1 . , θ cr = 0 . [55] α = 4 . , m = 1 . , θ cr = 0 . [56]Grass-like α ( (cid:107) u (cid:107) − u cr ) m α ≈ O ( − − , m = 3 , u cr = 0 . [57]Engelund and Hansen . c f θ / (cid:112) ( s − gd m − [54]Smart and Jäggi α (cid:16) d d (cid:17) . J . (cid:107) u (cid:107) ( θ − θ cr ) d m α = 8 , θ cr = 0 . [58]parameter [49]), θ cr is the critical dimensionless bed shear stress, d m is the representative grain diameter, s = ρ s /ρ is the relative density of the sediment with respect to water. The coefficients α , m and the criticalthreshold θ cr can be assigned by the user or adopted from literature (see Table 2). The Grass-like modelproposes a simple bedload transport formula, where (cid:107) q B (cid:107) is a function of the flow velocity magnitude, with u cr as critical threshold velocity. The coefficients α , m and the critical threshold u cr can be assigned by theuser or adopted from literature (Table 2). Engelund and Hansen [54] proposed a total transport formula foruniform bed material without a threshold condition for incipient motion. Local corrections of the sediment transport.
The morphodynamic module implements three corrections tothe basic Exner equation (3) to account for the influence of local characteristics of the flow and the bottomon the sediment transport. Namely, i) the influence of local slope on incipient motion, ii) the effect of lateralbed 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 horizontalbed. In case of a sloped bed in flow direction or transverse to it, the stability of grains is either increasedor reduced due to the gravity. The critical shear stress value can be adapted consequently to account forthe influence of local longitudinal and transversal slopes. A common approach is to scale the critical shearstress for almost horizontal bed θ cr with a correction factor k : θ ∗ cr = kθ cr . (4) Basement implements the correction factor k as proposed in [60] and [61]. Implementation details are9iven in the official documentation.The bedload direction can be corrected to account for two relevant morphodynamic processes linked tothe slope of the bed and the curvature of the flow. The deviation of the bedload direction from the flowdirection can thus be modelled as a deviation angle ϕ = ϕ b + ϕ c , sum of the correction angle for bed slope( ϕ b ) and curvature ( ϕ c ), as depicted in Fig. 1b and c. The bedload vector is then rotated with the rotationmatrix T ( ϕ ) , being T = cosϕ − sinϕsinϕ cosϕ , (5)where the angle is positive counterclockwise.The angle ϕ b is estimated with the approach proposed in [62] and [63] for the effect of the local transversalbed slope. In particular, the bedload direction deviates from the flow direction in presence of a localtransversal bed slope, due to the gravity acting on the bedload sediment particles (Fig. 1b). The bed loaddeviation ϕ b with respect to the flow is therefore evaluated as tan ϕ b = − N l (cid:114) θ cr θ · s · n q , for s · n q < , (6)where N l is an experimental lateral transport factor ( . ≤ N l ≤ . ), s = ( ∂ x z B , ∂ y z B ) is the local bedslope and n q is the unit vector perpendicular to q in downhill direction (Fig. 1b).The angle ϕ c accounts for the effect of a marked flow curvature. Due to three dimensional spiral flowmotion that establishes in curved flows, the bed load direction tends to point towards the inner side of thecurve, while the flow direction points towards the outer side (Fig. 1c). This curvature effect is taken intoaccount according to an approach proposed in [64], where the deviation angle ϕ c is determined as tan ϕ c = − N ∗ hR , (7)where h is the water depth, N ∗ is a curvature factor, and R denotes the radius of the river bend, positivefor curvature in counterclockwise direction. The curvature factor N ∗ mainly depends on bed roughness andassumes values N ∗ ≈ for natural streams [64], and values up N ∗ ≈ 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 totalvolume flux including porosity [ m / s ]. Similarly to the hydrodynamic case (§3.1.1), different approaches areadopted for the sediment sink, namely exact , available and infinity . 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-erodible10iver bottom, as in case of bedrock or concrete cover, a non-erodible fixed bed depth z rel (Fig. 1a) can beset. This threshold also determines the volume of sediment available for transport. The fixed bed elevationis defined relative to the initial bottom elevation z B with z rel ≤ . A number of environmental processes, such as pollutant, temperature or nutrient transport, can bemodelled assuming the passive advection and diffusion of a scalar quantity, in the form of dissolved orparticulated particles [e.g. 27]. The scalar advection-diffusion module allows for the simultaneous simulationof multiple passive species, up to a maximum of 5. The transport of a generic species c can be described bythe following advection-diffusion equation: ∂ t q c + ∂ x (cid:104) q x q c h − h ( K xx ∂ x φ c + K xy ∂ y φ c ) (cid:105) + ∂ y (cid:104) q y q c h − h ( K yx ∂ x φ c + K yy ∂ y φ c ) (cid:105) = S φ c , with c = [1 , , (8)where the unknown is q c , the specific mass of the species c . It can be expressed as q c = hφ c , with φ c thevolumetric concentration and h the water depth, as usual. The term S φ c is a net source of c and K ij [ m / s ]are the components of the 2D diffusion tensor. For the scalar advection-diffusion module, the closure relationships are used to model the contributionof external scalar input and output. In particular, the term S φ c represents an additional scalar mass fluxthat is added within a set of elements defined by regions. The source can be specified either as an imposedconcentration value or a total volumetric flux [ m / s ]. The behavior is analogous to the case of hydro- andmorphodynamics, §3.1.1 and §3.2.1.The terms K ij of the diffusion tensor vary considerably with respect to the physical nature of thetransported species. Diffusive transport is modelled in terms of both molecular diffusion K m and turbulentdispersion K tij , such that K ij = K m I ij + K tij , with I ij the identity matrix. The molecular diffusion isassumed as an isotropic Fickian process with constant coefficient K m . Turbulent dispersion is non-isotropic( K tij ) and scale with the friction velocity u ∗ = (cid:107) u (cid:107) /c f and water depth via a longitudinal α L and transversal α T non-dimensional coefficient. Suitable values for open channel flows in natural environments are α L =13and α T =1.2 [27].
4. Numerical solution
The numerical solution of the governing equations (1), (3) and (8) is sought in a finite volume frame-work, with a spatial discretization based on unstructured meshes (§4.1). For the temporal integration, anexplicit first order Euler scheme is used. In its basic configuration, the temporal integration proceeds in asynchronous-decoupled way for all the modules, meaning that the modules are independently integrated in11 i Γ i j Ω j l i j n i j xy x-y- 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). time with the same timestep (§4.5). The following sections detail the domain discretization strategy and theadopted numerical solver for the fluxes calculation of the three basic modules. The interested reader shouldrefer to the provided references for specific implementation details. The problem is discretised adopting a finite volume approach over unstructured triangular meshes. Aconforming triangulation T Ω of the computational domain Ω ⊂ R by elements Ω i such that T Ω = (cid:83) Ω i , isassumed. Given a finite volume element Ω i (Fig. 2), j = 1 , , is the set of indexes such that Ω j is a neighbourof Ω i . Γ ij is the common edge of two neighbour cells Ω i and Ω j , and l ij its length. n ij = ( n ij,x , n ij,y ) is theunit vector which is normal to the edge Γ ij and points toward the cell Ω j . The system of governing equations (1) can be cast in vectorial form as ∂ t U + ∂ x F x + ∂ y F y = S , (9)where left-handside terms of (9) are U = Hq x q y , F x = q x q x h + 12 gH − gHz B q x q y h , F y = q y q x q y hq y h + 12 gH − gHz B . (10)12he vector of source terms can be written as S ( U ) = S h + S fr ( U ) + S bed ( U ) , where S h = S h , S fr = − ghS fx − ghS fy , S bed = − gH∂ x z B − gH∂ y z B . (11)By integrating the governing system of equations (9) in the control volume V = [Ω i ] × [ t n , t n +1 ] , we obtainthe general update formula for the triangular element i : U n +1 i = U ni − ∆ t | Ω i | (cid:88) j =1 l ij [ F ij ] + ∆ t S i . (12)Problem unknowns at cell i and discrete time n are represented by cell averages U ni ; the numerical solutionsought at time t n +1 = t n + ∆ t is denoted by U n +1 i . In (12), F ij are the hydrodynamic fluxes estimated atthe cell interface ij , as in Fig. 2.To compute the fluxes 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 basicHLL scheme to account for the influence of intermediate contact waves. Further details on the HLLCapproach are available in Chapter 10 of [18]. The solver is proved to be robust and efficient in simulatingunsteady flows and the advection of passive tracers [27].The numerical discretization of the three terms of the source vector S (11) is conducted separately,according to the nature of each term. The external inflow/outflow contribution S h is added explicitly to thecontinuity equation, as it is not a function of the problem unknowns. The stiff friction source terms S fr ( 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 S bed ( U ) arediscretized using the modified-state approach proposed by [69]. This results in an easy and robust treatmentof complex topographies and wetting and drying problems [27]. 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 ∆ t of the hydrodynamic problem. The general update formula for the Exner problem reads: z n +1 Bi = z nBi − − p ∆ t | Ω i | (cid:88) j =1 l ij [ q Bij ] + ∆ tS bi , (13)with the same symbols introduced for (3). The term q Bij represents the normal sediment flux at the cellinterface ij (Fig. 2). 13or the numerical estimation of the term q Bij , a number of approaches are available in literature. In thecurrent version,
Basement implements an Approximate Riemann Solver of HLL-type [ sensu q Bij = λ + s q Bi − λ − s q Bj + λ + s λ − s ( z Bj − z Bi ) λ + s − λ − s , (14)where pedix i ( j ) refers to quantities evaluated at the corresponding cell (Fig. 2), and λ + s , λ − s are speedestimation of the morphological problem. We adopt the following closure [70]: λ − s = min ( λ i , λ j ) , λ + s = max ( λ i , λ j ) . (15)The expression for the terms λ and λ , calculated for both cell i or j reads: λ , = 12 (cid:32) u n − c ± (cid:115) ( u n − c ) + 4 ∂q B,n ∂q n c (cid:33) , (16)where u n is the normal velocity at the cell interface ij (Fig. 2) and c = (cid:112) ( gh ) is the so-called wave celerity. The scalar advection-diffusion problem is solved in a synchronous-decoupled way with respect to theshallow water problem (§4.2). We reformulate the governing equation (8) via a Cattaneo-type relaxationtechnique, as proposed by [27]. Two additional scalar conservation equations are then added to (8): ∂ t ψ cx − ∂ x φ c ζ = − ψ cx ζ , ∂ t ψ cy − ∂ y φ c ζ = − ψ cy ζ (17)where the new symbols are ζ , a positive and small relaxation time, and ψ cx and ψ cy are two auxiliaryvariables that recover ∂ x φ c and ∂ y φ c , respectively, for a sufficiently small ζ [27]. After a trivial substitutionof ψ cx ≈ ∂ x φ c and ψ cy ≈ ∂ y φ c into (8), the system composed by (8) and (17) can be rewritten in vectorialform as ∂ t U + ∂ x A x + ∂ y A y + ∂ x D x + ∂ y D y = S c + S rel (18)where the vectors U , A x , D x , S c and S rel read U = q c ψ cx ψ cy , A x = q c q x h , D x = − h ( K xx ψ cx + K xy ψ cy ) − q c ζh , S c = S φ c , S rel = − ψ cx ζ − ψ cy ζ , (19)with Q representing the conserved scalar quantities, A x is the advective fluxes vector, D x is the diffusive-relaxed fluxes vector, both in x direction. The scalar source terms are S c , whilst the source terms arising14orm the relaxation are S rel . For brevity, we omit the formulation for the y direction ( A y and D y ), whichis 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 aflux-splitting approach combining the advective and diffusive-relaxed fluxes and employs different solvers foreach. The HLLC solver, applied also for the hydrodynamic fluxes (§4.2) provides the advective componentof the scalar fluxes at the cell interface A ij . For the diffusive-relaxed component, the SVT technique derivesthe fluxes at the interface D ij directly from the Riemann invariants of a two non-linear wave Riemannproblem.Similarly to the hydro- and morphodynamic problems, the control volume V = [Ω i ] × [ t n , t n +1 ] is used tointegrate the governing system (18), in order to obtain the following scalar update formula at the element i U n +1 i = U ni − ∆ t | Ω i | (cid:88) j =1 l ij [ A ij + D ij ] + ∆ t ( S c + S rel ) , (20)where the fluxes A ij and 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 thenature of the terms. The simple scalar sources S c are computed with a first-order Euler scheme, while thestiff relaxation source terms S rel are integrated by means of a locally implicit Euler method. Numerical integration proceeds with a dynamic timestep ∆ t , evaluated at each time loop (Fig. 6b) thatfulfills the well-known Courant-Friedrichs-Lewy stability conditions [18]. In the current implementation, thecondition is expressed as: ∆ t = CF L min ≤ i ≤ N (cid:18) min ≤ j ≤ (cid:18) ρ ij λ ij (cid:19)(cid:19) , (21)where ρ ij is twice the distance between the edge j and the centroid of the cell i (Fig. 2), and N is the totalnumber of domain elements. The term λ ij is an estimation of the largest eigenvalue of the hydrodynamicproblem (1), namely λ ij = | u n | + √ gh with the simbology already introduced. The CFL coefficient rangesbetween 0 and 1: by default it is set to 0.9, if not specified otherwise.
5. Initial and boundary conditions
All modules require the user to define the initial conditions of the simulation. Two types of initialconditions are similarly available for all the modules:• region defined: user explicitly defines the initial values of the problem unknowns (e.g. water depthand specific discharge for hydrodynamics). Different values can be assigned to different region of thecomputational domain; 15 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 initialconditions. In this case, the domain will progressively fill with water, in relation to the assigned inflowboundary conditions or internal sources (§3.1.1).
The boundary conditions (hereinafter BCs) have different specifications for each core module (see follow-ing 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 Ω is defined bythe domain boundaries, as Γ , , . An external standard BC is dependent only on the local flow conditionsand on some user-defined rules. This represents the most common case, for example to define impermeablewalls 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 localBCs are defined also with information from a linked boundary. Typical example is a weir, where the flowdischarge 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 Ω , and not at the edges(Fig. 3). This BC type comes in handy in case of very large domain application, because it allows to testdifferent configurations of hydraulic structures (e.g. different locations of a weir or training wall), withoutthe 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 readercan refer to the official documentation for further details. 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:• Standard BCs: in addition to wall
BC, inflows (upstream BCs) and outflows (downstream BCs) can beassigned. As standard inflows, three options are provided (namely uniform, explicit, zhydrograph ). Forall cases, given a total volume discharge Q [ m / s ] or water surface elevation [m] as input, the inflowscondition for the mathematical unknowns of system (1) is set. For the standard outflows a value forthe water depth h must be specified. Possible options are: uniform conditions, hydraulic weir, ratingcurve, hydrograph and zerogradient (i.e. Neumann BC). It is worth remarking that the specific type16 utflow inflow t r a i n i n g w a l l weirweirdomain boundary standardlinkedinternal wall (default) Figure 3:
Example of modelling domain with different types of boundary conditions.
The computational domain Ω can include the river channels but also the surrounding floodplains. The domain is delimited by the external BCs: impermeablewalls (default type, Γ , , ) are depicted in solid black, while standard inflow and outflow are in dashed blue. A weir (dottedbrown) is modelled with an external linked BC. A training wall (thick solid green) is modelled with an internal BC. of upstream and downstream BCs should be selected depending on the local flow conditions (i.e. sub-or super-critical).• Linked BCs: this type of boundaries establish a link between two certain region of the domain wherethe governing equations are not solved. This type of BCs are particularly designed to simulate thebehaviour of hydraulic structures within the river channel, such as weirs, gates, bridges, spillways.• Internal BCs: they are fictitious boundaries defined as segments at the interfaces of some computationalcells. 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 wallconditions are applied until reaching a given threshold value (time or water depth) after which thewall is removed, and the SWE are solved. With the rating curve option (or h-Q relation ) a giveflow relation is applied on one side of the internal boundary, while on the other side, wall conditionsapply (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 theactual free unsteady flow (i.e. governing equations) over the (former) boundary. The sediment flow is defined as a specific bedload flux, which is averaged and evenly distributed at thedomain boundary conditions over the boundary length. In analogy with the hydrodynamic module, themorphodynamic boundaries are of type external standard and linked .17 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 fromthe hydrodynamic conditions under transport capacity conditions or iii) bed equilibrium condition,where the upstream bed elevation is kept constant. Two downstream BCs are available, allowingthe simulation of i) equilibrium condition and ii) check-dam. In this second option an equilibriumboundary condition is activated only if the bed level reaches a given threshold value, otherwise a walltype boundary is assumed.• Linked BCs: one BC is available. It allows for the simulation of sediment transport through givenhydrodynamic linked conditions, hence to ensure sediment continuity in the simulated channel.
Scalar BCs are defined in terms of concentration of total volumetric rate [ m / s ], evenly distributedthroughout the length of the relevant domain boundary. The implemented types are:• Standard BCs: three types are available. i) scalar inflow as a constant value; (ii) scalar inflow as atime-series and (iii) zerogradient (i.e. Neumann BC) outflow.
6. Software design
The standard modelling procedure involves three phases: the pre-processing phase, the numerical simu-lation 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 phasesand provide a short description of the different configuration and results file formats as used by
Basement .1. Pre-processing: in this phase the user is required to define the model domain and the input data. Themesh file (customized 2dm format, "MyMesh.2dm" in Fig. 4) contains the description of the triangularunstructured computational mesh. The file can be generated with BASEmesh, a Python script as wellas a QGIS plugin (see software website for details), or via grid generator software that supports the2dm format. In addition, further input data such as time series of water and sediment discharge (orother quantities) to be used as BCs can be provided (ASCII format, "MyData.txt" in Fig. 4).2. Numerical simulation: the actual simulation can be run via either CLI or GUI. In
Basement
Basement re-processing Numerical simulation Post-processing
Timeseries
External data Input fi les Output fi les Result data Topography Visualization(e.g. hydrology)
Figure 4:
Modelling workflow.
Pre-processing: generation of the computational mesh from topographical data and definitionof input time series. Numerical simulation:
Basement . Post-processing: elaboration and visualization of the results.
3. Post-processing: the XDMF file ("output.xdmf" in Fig. 4) can be opened with the Crayfish plugin forQGIS 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 areprovided at software website).
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 theGraphical User Interface or Command Line Interface. This modular design allows a customization of thesimulation workflow by the user and an efficient batch processing of
Basement steps. For instance, theprograms 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 monitoringthe three simulation steps. In particular, the GUI validates the configuration parameters and automaticallyadds required parameters where default values are available.1. The pre-simulation step focuses on the model definition. In particular, the "model.json" commandfile contains: i) physical properties, ii) initial conditions and iii) boundary conditions of the physicalproblem, and further iv) numerical parameters. The setup executable first reads the computationalmesh "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. The simulation is carried out on a selected computational backend (§6.4). It is driven by the commandfile "simulation.json" that contains the simulation parameters (e.g. execution time, output time anddesired output quantities). The program reads and executes the model "setup.h5" generated in theprevious step. The results of the simulation are stored in a second binary container: "results.h5".19 igure 5:
Software components and simulation steps.
The
Basement software is composed of a set of executables (redrectangles) driven by JSON configuration files (grey labels). Data is stored in HDF5 containers (green cylinders). In dashedarrows: the special actions of simulation re-run and restart .
3. The post-simulation step is configured using a command file "results.json" that contains the se-lected output format (currently only XDMF is supported). The output is then available for thepost-processing phase. (§6.1).When it is necessary to run a new simulation starting from the results of a previous one, two optionsare 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 indicatedan 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 takenfrom 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 optionis particularly useful in the case of large models (i.e. millions of computational cells): in this situation thepre-simulation step can be computationally demanding (i.e. tens of minutes). If the user only needs toextend the simulation duration, for example, then the
Re-run option allows to skip the pre-simulation step.
Basement aims to simulate different river processes with a high level of flexibility and efficiency: withthis in mind, we designed the software adopting a modular approach. The two core concepts of this designare
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 "HYDExternal Source", Fig. 6a). Modules are activated by the user in the pre-simulation step. An activatedmodule triggers the execution of a number of
Kernels throughout the simulation.20 nitialize HYD variables Hydrodynamics (HYD) (a) MODULES
HYD Boundary ConditionsHYD External SourceMorphodynamics (MOR)MOR Boundary Conditions Initialize MOR variables calculate timestepHYD boundary fl uxesMOR boundary fl uxesHYD internal fl uxesMOR internal fl uxesexternal HYD sourcesupdate all variables simulation startsimulation ends time loop (b) KERNEL SEQUENCING ... N Figure 6:
Examples of modules and associated kernels. (a)
At model setup, depending on the simulated physicalproblem, the user triggers the activation of a set of modules. (b)
The active modules trigger a unique kernel sequence tobe executed at simulation time to correctly simulate the requested processes. Each module corresponds to a different set ofkernels. 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 thecomputational domain or a subset of it. Depending on the specific task, Kernels can be scheduled for asingle execution (i.e. initialization Kernels) or for repeated execution in each iteration of the integrationtime loop (Fig. 6b). While the time loop is executed with a timestep ∆ t that satisfies certain stabilityconditions (§4.5), the Kernels can be scheduled for execution at different (larger) time intervals to reflectthe nature of the simulated processes.The architecture based on Modules and Kernels has two main advantages. First, it is flexible in that itallows users (and software developers) to easily add or remove specific modules without interfering with otherexisting 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). The parallelization strategy of the
Basement numerical core addresses two main aspects: i) the useof different technologies (i.e. computational backends) generated from the same, unique software source21 able 3:
Description of available backend types for
Basement v3.
Type Descriptionseq 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 CPUcode. 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. Tothis end the numerical core of
Basement
Basement currently supports multi-core CPUs and GPUs. When starting the simulation, the user canselect 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 thatthe choice of graphics processing units is currently limited to Nvidia (CUDA) cards. The precise requirementsare 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).
7. Results
A set of selected test cases (T1-T6) are proposed here to test the robustness, accuracy and efficiencyof the three basic modules. Table 4 summarizes the key features of each test case. The interested readercan refer to the official software documentation for further examples. Finally, §7.7 focuses on the softwareperformance and scalability. All the test cases are freely available at (link provided after paper acceptance).
The scope of this test is to assess the robustness and accuracy of the hydrodynamic solver when simu-lating a shock-type hydrodyamic wave travelling on a highly irregular and dry domain. The collapse of the22 able 4:
Main features of the adopted benchmarks.
ID Module Comparison Key features
T1 hydrodynamic field data i) highly unsteady shock wave generation and propa-gation, ii) wet-and-dry processes, iii) performance ofthe hydrodynamic moduleT2 morphodynamic lab data i) sediment transport with a transcritical 1D flow, ii)upstream BCs: uniform inflow, iii) downstream BCs:imposed water levelT3 morphodynamic lab data i) 2D dam-break in a complex domain, ii) sedimenttransport with an advancing wet-and-dry front, iii)downstream BCs: free outflow, iv) performance ofthe morphodynamic moduleT4 morphodynamic lab data i) erosion and deposition in a channel bend, ii) sed-iment transport direction correction, iii) upstreamBCs: unsteady hydrographT5 scalar advection-diffusion numerical sol i) 1D strong rarefaction waves, ii) conservation of asteady discontinuity of the scalar quantityT6 scalar advection-diffusion numerical sol i) 2D complex domain, ii) advection and diffusion oftwo different scalar quantities, iii) conservation andmixing of the scalars, iv) performance of the scalaradvection-diffusion module23alpasset dam, in the Reyran River Valley (Fréjus, France), represents a well-established hydrodynamicbenchmark for numerical models [e.g. 71, 72, 73]. In 1959, the 66.5 m high dam collapsed almost instan-taneously, generating an up to 40 m high flood wave that propagated down the Reyran valley, destroyedthe two villages Malpasset and Bozon and reached the Mediterranean Gulf 21 minutes later [72]. Thepropagation 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 for17 survey points, marked as P1 to P17 in Fig. 7 and the flood arrival time is known from three electrictransformer stations which have been destroyed by the flood wave. The locations of the transformer stationsare indicated as A, B and C in Fig. 7. Coordinates and recorded arrival times are listed in Table 5. Wemake 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 areset to walls (§5.2.1), with exception of the downstream boundary located in the Mediterranean Gulf, wherea fixed water level was set to 0 m. The initial conditions are a fixed water surface elevation of 100 m inthe reservoir, and dry conditions in the rest of inland domain. The initial velocity was set to . / s in theentire domain. In accordance with [71], the Manning’s friction coefficient was set to .
033 m − / s for thewhole domain. The CFL number was set to 0.9.The simulated maximum water levels are compared to the 17 field observations in Fig. 8, with overall goodagreement. The average relative error is 7.15%, with the largest observed at P13, with an overestimationof 30.6%. To highlight the effects of the topographical approximations of the Digital Elevation Model (andhence of the computational mesh), we compared recorder and simulated water level values as follows. Foreach punctual maximum water level recorded in field observations (blue triangles in Fig. 8), we comparedthe maximum simulated values of the spatial mean, maximum and minimum among the computational cellcontaining the observation point and its three neighbours (black and red series in Fig. 8). We expect lowerdiscrepancies between recorded and simulated values where the numerical values, hence the topographicalelevations, are spatially homogeneous. As a matter of fact, points P1, P7 and P13 (Fig. 8) have the largediscrepancy between measured and simulated values, but also the largest spatial variability of the numericalvalues (red shaded area). This suggests that such discrepancies relates more to the local topographicalapproximations of the DTM rather than to the numerical model.Observed and simulated times of flood arrival are given in Table 5. Simulated values are in goodagreement with measured ones for all the electrical transformer stations (ET). Simulated arrival times havea maximum relative error of 3.8% for ET B, corresponding to an absolute delay of 47 s. It is worth mentioningthat the friction value influences the simulated arrival times.24 lectr. transformersSurvey pointsBottom elevation [m]-20020406080100Water depth [m]01020304050 Figure 7:
T1: planar view of the Malpasset test case.
Visualization of the computational domain illustrating the bottomelevation (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, asreference.
P1 P2 P3 P4 P5 P6 P7 P8 P9 P10 P11 P12 P13 P14 P15 P16 P17Survey points020406080100 M a x i m u m w a t e r l e v e l [ m ] field datanumerical maxnumerical meannumerical min Figure 8:
T1: Malpasset dam-break wave maximum water level.
Numerical values are compared to the field data atthe survey points P1 to P17: simulated values are given as maximum, minimum (red lines) and mean value (black line) of thecomputational cell containing the survey point and its three neighbours. able 5: T1: Malpasset dam-break wave arrival times.
Observed and simulation time of flood arrival (TFA) and relativeerror (Err) for the three electrical transformer stations (ET) destroyed by the flood wave.
ET x y TFA obs
TFA sim
Err0 [-] [m] [m] [s] [s] [%]A 5550 4400 100 103 3B 11900 3250 1240 1287 3.8C 13000 2700 1420 1435 1
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 morphologicalsolver 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 computa-tional domain is a composed by a straight 6.9 by 0.5 m channel, representing the lower part of the originalexperimental flume, and is discretized with 24,612 triangular elements. The sediment has a characteristicdiameter of 1.65 mm. The water and sediment discharge at the upstream boundary are set to .
012 m s − and .
196 m s − (with porosity) respectively. The flume is at initial uniform flow conditions, characterizedby a supercritical flow. At t =0 s, a fixed water level of 0.2093 m is imposed at the downstream boundaryand sediment transport out of the domain is stopped. This results in the formation of an hydraulic jumpmoving upstream in the flume, and a subsequent downstream propagation of a sediment bore. The CFLnumber is set to 0.9, the bed porosity is assumed constant and equal to 0.42, and the simulation durationis 500 s.Fig. 9 shows the initial and final profiles of the simulated bed and water elevations. The solver reproduceswell the sharp transition between super- and sub-critical flow conditions. The position of the sediment frontin time is shown in Fig. 10, with good agreement between simulated and experimental values. 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 awell-know morphodynamic test for numerical models [e.g. 14, 70, 77]. The domain consists of a flat flumewith a non-symmetrical sudden enlargement (Fig. 11). The bed is composed of a coarse uniform sand with a26 W a t e r a nd b e d e l e v a ti on [ m ] bed elevation t=0 swater elevation t=0 sbed elevation t=500 swater elevation t=500 s Figure 9:
T2: propagation of a sediment bore.
Initial (dashed lines) and final (solid lines) longitudinal profiles of bedelevation (black) and water elevation (blue) for the propagation of a sediment bore test. T i m e [ s ] experimentalnumerical Figure 10:
T2: evolution in time of the sediment front position.
Blue triangles are the experimental values [75], whilethe black line is the numerical solution. nitial water depth = 25cm Dam Free outflow Y [m]
P1 P2 P3 P4P5 P6
X [m] CS1 CS2
Figure 11:
T3: planar view of the dam-break over mobile bed setup (deformed axis). Experimental and numericalresults are compared at survey points P1 to P6 during the simulation, and at cross-sections CS1 and CS2 at the end ofsimulation. median diameter of d m =1.82 mm. The initial conditions are defined by an horizontal layer of fully saturatedsand of thickness 0.1 m over the whole domain and an initial water storage of depth 0.25 m upstream of thedam, located at section x =3.0 m. At time t =0 s, the dam is suddenly removed, resulting in the propagationof a dam break wave with consequent sediment transport.The computational domain is discretized by unstructured triangular cells at different resolutions (followsin Tab. 6). Inviscid wall boundary conditions are set at the upstream and lateral domain boundaries, whilea free-outflow condition is used at the downstream outlet. The Manning coefficient is set to . − / s ,the sediment density and porosity are set to / m and 0.47, respectively. The sediment transport isevaluated with the MPM-like formula (Table 2), setting θ cr = 0.0495 for the critical Shields stress, α = 3.97and m = 1.5 for the remaining parameters [56]. The CFL number is set to 0.9. The numerical simulationslast 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 arrivaltime 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 flowpattern generate by multiple wave reflections while simulation proceeds in time, which potentially generatetri-dimensional flow structures. Nevertheless, obtained series are coherent with the ones of [14], where asecond-order accuracy model was employed.Numerical bed elevations after 12 s are compared to the experimental results in Fig. 13. The simulatedscour and deposition patterns are well reproduced. The magnitude of the scour at cross section CS1 ( y ≈ . m) matches well, with an underestimation of the deposition pattern ( y ≈ . m). At cross section28 w a t e r s u rf ace e l e v a ti on [ m ] P1P1P1P1
P2P2P2P2 w a t e r s u rf ace e l e v a ti on [ m ] P3P3P3P3
P4P4P4P4 w a t e r s u rf ace e l e v a ti on [ m ] P5P5P5P5
P6P6P6P6
Figure 12:
T3: experimental and numerical water elevation at six survey points.
Sub-panels represent the surveypoints P1 to P6 as in Fig. 11: experimental points are given as full circles, whilst numerical results as lines, with four differentmesh sizes (27k, 54k, 108k, 216k computational cells). .0 0.1 0.2 0.3 0.4 0.5transversal coordinate y [m]0.060.080.100.120.14 b e d e l e v a ti on z b [ m ] CS1 at x = 4.1 m y [m]0.060.080.100.120.14 b e d e l e v a ti on z b [ m ] CS2 at x = 4.4 m 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 fourdifferent mesh sizes (27k, 54k, 108k, 216k computational cells).
CS2, the simulated deposition magnitude matches well with the experimental one, but with a small shifttoward the lateral boundary.
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 experimentfrom [79], already adopted as morphodynamic benchmark test [e.g. 80]. The flume is U-shaped, with abend of ° having a costant radius along the center line of R c =4 m. The cross section is rectangular withwidth 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 witha simplified (triangular) flood hydrograph, having a base flow of .
02 m s − , and peak flow of .
053 m s − .The rising and falling limbs last 100 and 200 minutes, respectively. Afterwards a constant baseflow waskept for another 100 minutes. During the experimental run, a steady point bar in the inner side of the beddevelops 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 to0.4, and used the MPM-like formula with parameters from [56], as in Table 2. The lateral slope factor N l isequal to 1.4 and the curvature factor N ∗ is set to 11 (6 and 7). At the numerical domain boundaries uniformflow 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 bedvariation with respect to the initial flat configuration is scaled with the approaching (i.e. upstrem reach)30 ° 90° 180° . . . . - . - . - . Figure 14:
T4: scour and deposition on a channel bend.
Planar view of relative bed change ∆ z b / h , with h approachingflow depth. Numerical results on the left, laboratory results from [79] on the right. outer sideinner side R / R c [-]1.00.50.00.51.0 ∆ z b / h [-] experimentalnumerical Figure 15:
T4: scour and deposition at one cross-section.
Cross-sectional view (at °) of relative bed change ∆ z b / h ,with h approaching flow depth. In the x -axis the radial coordinate R is scaled with the center-line curvature radius R c ; blackline is the numerical solution, blue triangles are from the laboratory experiment of [79]. flow depth h . The magnitude of scours and depositions for the numerical run ranges between -0.75 and0.75, matching fairly well with the experimental values. Also the positioning of the point bar, with themaximum deposition anticipating the middle of the bend ( °) is well reproduced numerically.Fig. 15 shows the cross-sectional profile of the relative bed change for the numerical simulation andthe experimental run in the middle of the flume bend (at °). The numerical profile reproduces well theexperimental trend. This test demonstrates the software capability in simulating an unsteady morphologicalbuilding process, i.e. a point bar development during a flood. Such process can be well reproduced onlyby implementing suitable corrections of the sediment transport direction due to gravity and curvature, aspresented in §3.2.1. 31 .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 chosentest is an idealized one-dimensional problem, but nevertheless it is particularly challenging for a pletora ofnumerical scheme [18]. The domain is a simple, flat-bed, channel 100 m long and 0.1 m wide. The domainis deliberately chosen very narrow, to mimic a 1D setup, given that the exact solution of the problem isavailable 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 φ present adiscontinuity: q x = − . / s if x <
50 m , q x = 3 . / s otherwise ,φ = 1 . if x <
50 m , φ = 0 . otherwise ,h = 1 . ∀ x. (22)The domain is discretized with 1362 computational cells, lateral walls are reflective and inviscid, whilsttransparent boundary conditions are set at beginning and end of the channel. The CFL is set to 0.95 and thesimulation timeout is t =2.5 s. As the simulation starts, two strong rarefaction waves start diverging fromthe 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 contactwave 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 inFig. 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 underlineshow the numerical solution correctly approximates the exact solution in all the domain sections. The scalardiscontinuity is perfectly maintained throughout the simulation, confirming the accurate resolution of thesteady contact wave. With this test we assess the solver capability in correctly preserving the liquid and scalar species masswhen simulating the advection and diffusion of species during a dam-break phenomena. The test is partic-ularly 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]. Thedomain is composed by a rectangle [0; 75] × [ −
15; 15] m. The bottom η ( x, y ) is fixed during the simulation,32 Water depth num exact
Concentration num exact (a) −3−2−10123 −3−2−10123Longitudinal coordinate x [m]20 30 40 50 60 70 8020 30 40 50 60 70 80
Specific discharge num exact (b)Figure 16:
T5: comparison between numerical and exact solution of a steady contact wave.
Solution is givenat timeout t =2 s. Exact solutions are depicted with lines, whereas numerical values are symbols, decimated for the sake ofvisualization. Panel (a) shows the water depth (blue diamonds) and the scalar concentration (red circles), panel (b) shows thespecific discharge on the x direction (green circles). and defined as η ( x, y ) = max (0 , η ( x, y ) , η ( x, y ) , η ( x, y )) , with η ( x, y ) = 1 − (cid:112) ( x − + ( y + 9) ,η ( x, y ) = 1 − (cid:112) ( x − + ( y − ,η ( x, y ) = 3 − (cid:112) ( x − . + y . (23)The initial conditions are given by: h = 1 . if x <
16 m , h = 0 .
125 m otherwise ,φ = 0 . if x <
16 m , φ = 0 otherwise ,φ = 0 if x <
16 m , φ = 0 . otherwise ,q x = q y = 0 m / s everywhere , (24)presenting a virtual dam at x =16 m separating two discontinuous volumes of water and scalar mass. Herethe domain is discretized with 492,277 triangular cells with a maximum of characteristic length of 0.1 m. Thehydrodynamics setup features reflective wall boundaries, a CFL of 0.95 and frictional sources compatiblewith a Manning coefficient of n = .
01 m − / s . The scalar setup features two initially unmixed species, bothwith a constant and isotropic diffusion coefficient K c = .
25 m / 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 thetwo 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 flat33opography, is still lagging by approximately 15 m (Fig. 17b). At this point, the reflected bores propagateupstream and further mix both species, symmetrically around the x axis. By t =50 s these bores overcomethe 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 frictionsources gain relevance and dissipate most of the kinetic energy in the flow, with a near-static solution beingobtained 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 fullyhomogeneous across the domain (Fig. 18).The model is fully conservative, with the total liquid and scalar mass preserved during the entire simu-lation. As the simulation approaches the lake-at-rest conditions, the observed quantities correctly convergeto their resting values of h =0.364 m, φ =0.386 and φ =0.023 (Fig. 18). The performance and scalability of the software depends not only on the implemented parallelizationstrategies but also on the physical model to be reproduced. In general, "simpler" models, i.e. only fewsimulated physical processes, are likely to show higher computational performances. To test
Basement ’scomputational 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 com-putational meshes. The sizes of these meshes are given in Table 6 and have been chosen to cover a broadrange of spatial resolutions, ranging from thousands to hundred of thousands computational cells. Thesimulations have been run with a set of computational backends. In particular, CPU-based simulations havebeen performed on an Intel Xeon Gold 6154 (3.00GHz) workstation equipped with 36 cores (two socketswith 18 physical cores each), whilst GPU-based simulations have been run on three GPUs (GeForce GTX1050 Ti, GeForce GTX 1080 Ti, and Tesla P100; see Table 7 for the main characteristics). Moreover, thesimulations have been benchmarked in both single and double precision mode. The GPUs were integratedin a workstation with a 32-core Intel Xeon Gold 5218 (2.30GHz) processor (two sockets with 16 physicalcores each).For a given a mesh size, the speedup achieved by a parallelized backend p is computed using the formulaspeedup = T s /T p , where T s ( T p ) is the total computational time used by the serial (parallel) backend. Theresults for all the investigated benchmarks are depicted in Fig. 19.As anticipated, the speedup depends on the simulated processes. Comparing the speedup values amongdifferent benchmarks in Fig. 19, cases T1 (hydrodynamics) and T3 (morphodynamics) show higher values onaverage than T6 (advection-diffusion). Such results are expected, given the increased complexity (numberof equations and operations to be solved) of T6. 34
15 30 45 60 75 − x ( m ) y ( m ) a) t =0 s − x ( m ) y ( m ) b) t =20 s − x ( m ) y ( m ) c) t =50 s . . . . H (m) 0.5 0.25 0.0 0.05 0.1 φ (-) φ (-) Figure 17:
T6: dam break over complex topography . 3D visualization of the water surface elevation H (left panels) andplanar view of concentration distribution ( φ and φ ) for both scalar species 1 and 2 (right panes). Subpanels (a,b,c) reportdifferent simulation timeout. . . . . . . t (h) 0.0 0.5 1 1.5 2 2.5 3 . . . . . . t (h) h ( m ) q x ( m /s ) φ (-) φ (-) Figure 18:
T6: time series at two locations.
Evolution of hydrodynamic quantities ( h and q x ) and scalar concentrations φ (red) and φ (blue) in time at position ( x, y ) = (7 . , . m (left) and ( x, y ) = (67 . , . m (right). Dashed lines representtheoretical scalar concentration values, at rest.Table 6: Number of computational cells for the performance and scalability benchmarks.
Mesh ID T1 T3 T6
24 945 27 444 24 388
52 102 47 187 49 155
101 417 109 344 98 163
499 060 218 912 196 829
Table 7:
Characteristics of the benchmarked Nvidia GPU cards.
Type Architecture CUDA Cores
GeForce GTX 1050 Ti Pascal 768GeForce GTX 1080 Ti Pascal 3584Tesla P100 Pascal 358436 S p ee dup Number of cells Tesla P100 STesla P100 DGTX1080Ti SGTX1080Ti D
GTX1050Ti S
GTX1050Ti DOpenMP-32OpenMP-16OpenMP-8OpenMP-4OpenMP-2 (a) T1
100 0 50000 100000 150000 200000 250000 S p ee dup Number of cells Tesla P100 S
Tesla P100 DGTX1080Ti S
GTX1080Ti DGTX1050Ti SGTX1050Ti DOpenMP-32OpenMP-16OpenMP-8
OpenMP-4OpenMP-2 (b) T3 S p ee dup Number of cells
Tesla P100 S
Tesla P100 DGTX1080Ti SGTX1080Ti DGTX1050Ti S
GTX1050Ti D
OpenMP-32OpenMP-16OpenMP-8OpenMP-4OpenMP-2 (c) T6Figure 19:
Speedup of backends for varying mesh size, test cases T1, T3 and T6.
The benchmarks were executedwith different degrees of parallelism on the CPU (using OpenMP) and the GPU. The final "S" and "D" in GPU series denotesingle and double precision, respectively.
Basement ’s parallelization can be evaluated in more detail by comparingthe speedup values along the vertical axis of the plots. In the following we focus on benchmark T3 (morpho-dynamics) which shows an "intermediate" scalability among the three benchmarks (Fig. 19b). Looking atthe mesh with 47k elements as an example, the CPU-based family (i.e. OpenMP on multiple cores) shows aspeedup efficiency (i.e. speedup / ncores · ) of 87% with 2 CPU cores (Speedup=1.7), and of 58% with 32cores (Speedup=18.5), with an average efficiency of 74%. Basement performs even better on some GPUcards. 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 bythe tested GPUs ranges between 7 and 60.The benchmarks in Fig. 19 also show how the maximum speedup changes with mesh size, computationalbackend and simulated processes. For all three cases, the CPU-based parallelizations show a mild speedupincrease with an increasing number of computational cells. The dependency on the mesh size is slightly morepronounced when the number of computational cores is increased. This reflects the fact that CPU-basedsolutions 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), theparallelization efficiency for 2-4 CPU cores is almost constant for all mesh sizes and above 80%. On theother 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 problemsize, but also on the simulated processes. In benchmark T1 (Fig. 19a), the speedup clearly increases withproblem size. This can be explained with the overhead of GPU parallelization, which becomes more andmore negligible with increasing domain size. Conversely, benchmark T6 (Fig. 19c) shows little impact of thedomain 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 greaterthan 10 between the single and double precision versions. Of course, the adequate choice depends on therequirements of the specific application.The analysis above shows
Basement ’s performance on different computational backends and underlinesthe differences when simulating different processes. The results summarized in Fig. 19 can also serve asa guideline for the interested reader/user when choosing an appropriate computational configuration for agiven application. Finally it is worth highlighting that all the tested hardware configurations can be easilyinstalled in standard office workstations. 38 . Conclusions
In this paper we introduced the main features of
Basement version 3, a freeware tool for river simula-tion.
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 scalehydrodynamic processes modelled with several hundreds of thousand elements in good agreement with themeasurements. On the opposite end of the spectrum, the morphological solver is able to handle demandingsediment transport scenarios well, albeit with known limitations. With the scalar advection-diffusion modulea 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 onrequest in
Basement ’s pre-simulation step. The advantage of this approach is that only the required tasksare 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 computationalpower. Such advantages are reflected in the presented benchmarks. Given a large enough domain, thesoftware 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 theexisting basic modules. As an example, the modelling of the sediment transport in presence of non uniformsediment size and the simulation of water temperature dynamics are in implementation phase. On theother hand, efforts are dedicated also to develop novel modelling solutions for river processes such as thebio-morphodynamic feedbacks between vegetation and sediment transport. Table 2 of the SupplementaryMaterial provides an overview of under development features. The modularity of the development frameworkallows also for further refinement of single specific numerical solvers and the implementation of high-orderschemes 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, developmentand implementation was done by DV, SP, LV, MB, MW, with the coordination and supervision of DFVand AS. Implementation and testing of the reported features was done by MB, DV, MW and DC. Themanuscript was conceptualized by DV, AS and DFV, and drafted by DV, with support of MB and MW. AllAuthors contributed to the manuscript review. 39 onflicts 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 theEnvironment (BAFU). 40 eferencesReferences [1] D. J. Gilvear, M. T. Greenwood, M. C. Thoms, P. J. Wood (Eds.), River Science, John Wiley & Sons, Ltd, doi:10.1002/9781118643525, 2016.[2] S. K. Brewer, T. A. Worthington, R. Mollenhauer, D. R. Stewart, R. A. McManamay, L. Guertault, D. Moore, Synthesizingmodels useful for ecohydrology and ecohydraulic approaches: An emphasis on integrating models to address complexresearch questions, Ecohydrology 11 (7) (2018) e1966, ISSN 19360584, doi:10.1002/eco.1966.[3] 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 InternationalRiver Interface Cooperative (iRIC) software, Earth Surface Processes and Landforms 45 (1) (2019) 11–37, doi:10.1002/esp.4653.[4] A. Crosato, M. S. Saleh, Numerical study on the effects of floodplain vegetation on river planform style, Earth SurfaceProcesses and Landforms 36 (6) (2011) 711–720.[5] D. Sharma, A. Kansal, Assessment of river quality models: a review, Reviews in Environmental Science andBio/Technology 12 (3) (2012) 285–311, doi:10.1007/s11157-012-9285-8.[6] R. D. Williams, J. Brasington, D. M. Hicks, Numerical Modelling of Braided River Morphodynamics: Review and FutureChallenges, Geography Compass 10 (3) (2016) 102–127, doi:10.1111/gec3.12260.[7] S. J. Dugdale, D. M. Hannah, I. A. Malcolm, River temperature modelling: A review of process-based approaches andfuture directions, Earth-Science Reviews 175 (2017) 97–113, doi:10.1016/j.earscirev.2017.10.009.[8] 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.[9] A. P. Zischg, M. Mosimann, D. B. Bernet, V. Röthlisberger, Validation of 2D flood models with insurance claims, Journalof hydrology 557 (2018) 350–361.[10] W. A. Marcus, M. A. Fonstad, Remote sensing of rivers: the emergence of a subdiscipline in the river sciences, EarthSurface Processes and Landforms 35 (15) (2010) 1867–1872, doi:10.1002/esp.2094.[11] J. Thomas Steven Savage, F. Pianosi, P. Bates, J. Freer, T. Wagener, Quantifying the importance of spatial resolution andother factors through global sensitivity analysis of a flood inundation model, Water Resources Research 52 (11) (2016)9146–9163.[12] G. Pasternack, 2D modeling and ecohydraulic analysis, University of California at Davis, California, ISBN 9781466320093,2011.[13] 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.[14] A. Siviglia, G. Stecca, D. Vanzo, G. Zolezzi, E. F. Toro, M. Tubino, Numerical modelling of two-dimensional morphody-namics 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.[15] J. Wyrick, A. Senter, G. Pasternack, Revealing the natural complexity of fluvial morphology through 2D hydrodynamicdelineation of river landforms, Geomorphology 210 (2014) 14–22, doi:10.1016/j.geomorph.2013.12.013.[16] M. Guan, N. Wright, P. Sleigh, S. Ahilan, R. Lamb, Physical complexity to model morphological changes at a naturalchannel bend, Water resources research 52 (8) (2016) 6348–6364, doi:10.1002/2015WR017917.[17] 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.
18] E. Toro, Shock-capturing methods for free-surface shallow flows, Wiley and Sons Ltd, ISBN 0471987662, doi:10.1080/00221680309499935, 2001.[19] M. Nahorniak, J. Wheaton, C. Volk, P. Bailey, M. Reimer, E. Wall, K. Whitehead, C. Jordan, How do we efficientlygenerate high-resolution hydraulic models at large numbers of riverine reaches?, Computers & Geosciences 119 (November2017) (2018) 80–91, ISSN 00983004, doi:10.1016/j.cageo.2018.07.001.[20] 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.[21] S. Dazzi, R. Vacondio, P. Mignosa, Internal boundary conditions for a GPU-accelerated 2D shallow water model: Imple-mentation and applications, Advances in Water Resources 137 (2020) 103525.[22] 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.[23] S. Dazzi, R. Vacondio, A. D. Palù, P. Mignosa, A local time stepping algorithm for GPU-accelerated 2D shallow watermodels, Advances in Water Resources 111 (2018) 274–288, doi:10.1016/j.advwatres.2017.11.023.[24] K. G. Powell, P. L. Roe, J. Quirk, Adaptive-Mesh Algorithms for Computational Fluid Dynamics, in: Algorithmic Trendsin Computational Fluid Dynamics, Springer New York, 303–337, doi:10.1007/978-1-4612-2708-3_18, 1993.[25] F. Carraro, D. Vanzo, V. Caleffi, A. Valiani, A. Siviglia, Mathematical study of linear morphodynamic acceleration andderivation of the MASSPEED approach, Advances in Water Resources 117 (2018) 40–52.[26] J. A. Morgan, N. Kumar, A. R. Horner-Devine, S. Ahrendt, E. Istanbullouglu, C. Bandaragoda, The use of a morphologicalacceleration factor in the simulation of large-scale fluvial morphodynamics, Geomorphology 356 (2020) 107088, doi:10.1016/j.geomorph.2020.107088.[27] D. Vanzo, A. Siviglia, E. F. Toro, Pollutant transport by shallow water equations on unstructured meshes: Hyperbolizationof the model and numerical solution via a novel flux splitting scheme, Journal of Computational Physics 321 (2016) 1–20.[28] A. Afzal, Z. Ansari, A. R. Faizabadi, M. K. Ramis, Parallelization Strategies for Computational Fluid DynamicsSoftware: State of the Art Review, Archives of Computational Methods in Engineering 24 (2) (2016) 337–363, doi:10.1007/s11831-016-9165-4.[29] 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.[30] 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.[31] 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.[32] 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.[33] 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.[34] 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.[35] J. Hou, Y. Kang, C. Hu, Y. Tong, B. Pan, J. Xia, A GPU-based numerical model coupling hydrodynamical and morpho-logical processes, International Journal of Sediment Research 35 (4) (2020) 386–394, doi:10.1016/j.ijsrc.2020.02.005.[36] M. J. Castro, S. Ortega, M. de la Asunción, J. M. Mantas, J. M. Gallardo, GPU computing for shallow water flowsimulation based on finite volume schemes, Comptes Rendus Mécanique 339 (2-3) (2011) 165–184, ISSN 16310721, doi:
57] A. J. Grass, Sediment transport by waves and currents, University College, London, Dept. of Civil Engineering, 1981.[58] G. M. Smart, M. N. R. Jaeggi, Sediment Transport on Steep Slopes, VAW-Mitteilung 64, Versuchsanstalt für Wasser-bau,Hydrologie und Glaziologie (VAW). Zürich, ETH Zürich., 1983.[59] 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.[60] L. C. van Rijn, Handbook Sediment Transport by Current and Waves, Delft Hydraulics Laboratory, Delft, The Netherlands,1989.[61] X. Chen, J. Ma, S. Dey, Sediment Transport on Arbitrary Slopes: Simplified Model, Journal of Hydraulic Engineering136 (5) (2010) 311–317, ISSN 0733-9429, doi:10.1061/(ASCE)HY.1943-7900.0000175.[62] S. Ikeda, Lateral bed-load transport on side slopes, Journal of the Hydraulics Division, ASCE 108 (11) (1982) 1369–1373.[63] A. M. Talmon, N. Struiksma, M. Van Mierlo, Laboratory measurements of the direction of sediment transport on transversealluvial-bed slopes, Journal of Hydraulic Research 33 (4) (1995) 495–517, ISSN 0022-1686, doi:10.1080/00221689509498657.[64] F. Engelund, Flow and bed topography in channel bends., Journal of the Hydraulics Division ASCE 100 (11) (1974)1631–1648.[65] I. L. Rozovskii, Flow of Water in Bends of Open Channels, Academy of Science of the Ukrainian S.S.R, Institute ofHydrology and Hydraulic Engineering, 1961.[66] E. F. Toro, M. Spruce, W. Speares, Restoration of the contact surface in the HLL-Riemann solver, Shock waves 4 (1)(1994) 25–34.[67] 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.[68] D. Vanzo, Eco-hydraulic quantification of hydropeaking and thermopeaking: development of modeling and assessmenttools, Ph.D. thesis, University of Trento, Trento, Italy, 2015.[69] A. Duran, Q. Liang, F. Marche, On the well-balanced numerical discretization of shallow water equations on unstructuredmeshes, Journal of Computational Physics 235 (2013) 565–586, ISSN 00219991, doi:10.1016/j.jcp.2012.10.033.[70] S. Soares-Frazão, Y. Zech, HLLC scheme with novel wave-speed estimators appropriate for two-dimensional shallow-waterflow on erodible bed, International Journal for Numerical Methods in Fluids 66 (8) (2011) 1019–1036, ISSN 02712091,doi:10.1002/fld.2300.[71] J.-M. Hervouet, A. Petitjean, Malpasset dam-break revisited with two-dimensional computations, Journal of hydraulicresearch 37 (6) (1999) 777–788.[72] A. Valiani, V. Caleffi, A. Zanni, Case study: Malpasset dam-break simulation using a two-dimensional finite volumemethod, Journal of Hydraulic Engineering 128 (5) (2002) 460–472.[73] J. Singh, M. S. Altinakar, Y. Ding, Two-dimensional numerical modeling of dam-break flows over natural terrain using acentral explicit scheme, Advances in Water Resources 34 (10) (2011) 1366–1375.[74] S. Cordier, M. H. Le, T. M. De Luna, Bedload transport in shallow water models: Why splitting (may) fail, how hyper-bolicity (can) help, Advances in Water Resources 34 (8) (2011) 980–989.[75] M. Bellal, B. Spinewine, C. Savary, Y. Zech, Morphological evolution of steep-sloped river beds in the presence of ahydraulic jump: Experimental study, in: XXX IAHR Congress, Citeseer, 133–140, 2003.[76] 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.[77] C. Juez, J. Murillo, P. García-Navarro, A 2D weakly-coupled and efficient numerical model for transient shallow flow andmovable bed, Advances in Water Resources 71 (2014) 93–109.[78] 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.
79] C.-l. Yen, K. T. Lee, Bed topography and sediment sorting in channel bend with unsteady flow, Journal of HydraulicEngineering 121 (8) (1995) 591–599.[80] K. Kaveh, M. Reisenbüchler, S. Lamichhane, T. Liepert, N. D. Nguyen, M. D. Bui, P. Rutschmann, A Comparative Studyof Comprehensive Modeling Systems for Sediment Transport in a Curved Open Channel, Water 11 (9) (2019) 1779.[81] 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.79] C.-l. Yen, K. T. Lee, Bed topography and sediment sorting in channel bend with unsteady flow, Journal of HydraulicEngineering 121 (8) (1995) 591–599.[80] K. Kaveh, M. Reisenbüchler, S. Lamichhane, T. Liepert, N. D. Nguyen, M. D. Bui, P. Rutschmann, A Comparative Studyof Comprehensive Modeling Systems for Sediment Transport in a Curved Open Channel, Water 11 (9) (2019) 1779.[81] 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.