SW2D-Lemon: A New Software for Upscaled Shallow Water Modeling

J. G. Caldas Steinstraesser, et al., Advances in hydroinformatics, 10.1007/978-981-19-1600-7_2 SW2D-LEMON: A NEW SOFTWARE FOR UPSCALED SHALLOW WATER MODELING Joao Guilherme Caldas Steinstraesser Inria, IMAG, Univ Montpellier, CNRS, Montpellier, France Carole Delenne Univ Montpellier, HSM, CNRS, IRD, Inria, Montpellier, France Pascal Finaud-Guyot Univ Montpellier, HSM, CNRS, IRD, Inria, Montpellier, France Vincent Guinot Univ Montpellier, HSM, CNRS, IRD, Inria, Montpellier, France Joseph Luis Kahn Casapia Inria, IMAG, Univ Montpellier, CNRS, Montpellier, France Antoine Rousseau Inria, IMAG, Univ Montpellier, CNRS, Montpellier, France ABSTRACT We present a new multi-OS platform named SW2D-LEMON (SW2D for Shallow Water 2D) developed by the LEMON research team in Montpellier. SW2D-LEMON is a multi-model software focusing on shallow waterbased models. It includes an unprecedented collection of upscaled (porosity) models used for shallow water equations and transport-reaction processes. Porosity models are obtained by averaging the twodimensional shallow water equations over large areas containing both a water and a solid phase. The size of a computational cell can be increased by a factor 10 to 50 compared to a 2D shallow water model, with CPU times reduced by 2 to 3 orders of magnitude. Applications include urban flood simulations as well as flows over complex topography. Besides the standard shallow water equations (the default model), several porosity models are included in the platform: (i) Single Porosity, (ii) Dual Integral Porosity, and others are currently under development such as (iii) Depth-dependent Porosity model. Various flow processes (friction, head losses, wind, momentum diffusion, precipitation/infiltration) can be included in a modular way by activating specific execution flags. We recall here the governing equations as well as numerical aspects and present the software features. Several examples are presented to illustrate the potential of SW2D. 1. INTRODUCTION In this paper, we present a new multi-OS platform named SW2D-LEMON (SW2D for Shallow Water 2D) developed by the LEMON research team in Montpellier. Simulating urban floods and free surface flows in wetlands requires considerable computational power. Twodimensional shallow water models are needed. Capturing the relevant hydraulic detail often requires computational cell sizes smaller than one meter. For instance, meshing a complete urban area with a sufficient accuracy would require 106 to 108 cells, and simulating one second often requires several CPU seconds. This makes the use of such model for crisis management impossible. Similar issues arise when modelling wetlands and coastal lagoons, where large areas are often connected by an overwhelming number of narrow channels, obstructed by vegetation and a strongly variable bathymetry. Describing such channels with the level of detail required in a 2D model is impracticable. A new generation of models overcoming this issue has emerged over the last 20 years: porosity-based shallow water models. They are obtained by averaging the two-dimensional shallow water equations over large areas containing both water and a solid phase [8]. The size of a computational cell can be increased by a factor 10 to 50 compared to a 2D shallow water model, with CPU times reduced by 2 to 3 orders of magnitude. Research on porosity-based shallow water models has accelerated over the past 15 years, with significant contributions led by V. Guinot and our research team [1], [2], [3], [4], [13], [15], [20]. SW2D-LEMON is the software continuation of this 15 years research production. It follows the original research code first developed by V. Guinot in Fortran 90. SW2D-LEMON is multi-platform (Linux, MacOS, Windows) and includes a convenient graphical user interface (GUI) together with an unprecedented collection of upscaled (porosity) models used for shallow water equations and transport-reaction processes. Applications include urban flood simulations as well as flows over complex topography with unstructured grids. Besides the standard shallow water equations (the default model), several porosity models are already included in the platform such as Single Porosity and Dual Integral Porosity. Various flow processes (friction, head losses, wind, momentum diffusion, precipitation/infiltration) can be included in a modular way by activating specific execution flags. The paper is organized as follows: in Section 2, we recall the traditional shallow water equations as obtained by Saint-Venant, before defining the notion of porosity, presenting the single and dual porosity models and then the (standard) features of the numerical schemes used in the code. Section 3 describes the SW2D-Lemon workflow by presenting the inputs and outputs as well as the different options available to the modeler. Section 4 presents the three different licenses provided for the software before illustrating the installation process under Windows. Finally, we end this article with a presentation of the test cases that are provided in version 1.0 of the code, that will be made available to the Polytech'Montpellier engineering school from September 2021. For the sake of simplicity, the SW2D acronym will hold for SW2D-Lemon hereafter. For more information on SW2D-LEMON, please visit http://sw2d.inria.fr 2. MODELLING 2.1 Equations and source terms The current version of the SW2D simulation engine incorporates the Shallow Water Equations (SWEs), the Single Porosity (SP) [1] and the Dual Integral Porosity (DIP) [2] models. Other models published by members of the LEMON team, such as the Multiple Porosity (MP) [3] and Depth-Dependent Porosity (DDP) [4] models, are currently under integration and testing and should be released soon. The governing equations for these three models can be written in the form 𝜕!𝐮+𝐌%∇.𝐅=𝐬 (1) where t is time, u and s are respectively the conserved variable and source term vectors, and F and M are respectively the flux and inertia tensors. The SWE model, an extension of Saint Venant’s one-dimensional equations [5], uses the following definitions for u, s, F and M: 𝐮= + ℎ ℎ𝑢 ℎ𝑣 / ,𝐬= 1 𝑅−𝐼 𝑔ℎ 6 𝑆",$ −𝑆%,$ 8 −𝐼𝑢+𝑊$ 𝑔ℎ 6 𝑆",& −𝑆%,& 8 −𝐼𝑣+𝑊& : %% (2a) 𝐅= 1 ℎ𝑢 ℎ𝑣 ℎ𝑢'+𝑔ℎ'/2 ℎ𝑢𝑣 ℎ𝑢𝑣 ℎ𝑣'+𝑔ℎ'/2 : ,𝐌=𝐈𝐝 (2b) where 𝑔 is the gravitational acceleration, ℎ is the water depth, 𝑢 and 𝑣 are respectively the 𝑥 − and 𝑦 − components of the flow velocity, 𝐼 is the infiltration rate, 𝑅 is the rainfall intensity, 𝑆!,# and 𝑆$,# ( 𝑋 = 𝑥, 𝑦 ) are respectively the bottom and friction slopes in the 𝑋 − direction, and 𝑊 # ( 𝑋 = 𝑥, 𝑦 ) is the wind drag specific force in the 𝑋 − direction. An eigenvalue analysis of this model is provided in [6]. The wind drag force is computed using Smith and Banke’s model [7]. The SP model was first introduced in a depth-dependent version [7] to model the effects of subgrid scale topography. A depth-independent version for the modelling of urban floods was presented in nonconservation form [8] and later adapted for shock-capturing finite volume methods and discontinuous urban feature properties [1]. In this model, u, s, F and M are defined as 𝐮=𝜙 + ℎ ℎ𝑢 ℎ𝑣 / ,𝐬= 1 𝑅−𝜙𝐼 𝜙𝑔ℎ 6 𝑆",$ −𝑆%,$ 8 +𝑔ℎ𝜕$𝜙−𝜙𝐼𝑢+𝜙𝑊$ 𝜙𝑔ℎ 6 𝑆",& −𝑆%,& 8 +𝑔ℎ𝜕&𝜙−𝜙𝐼𝑣+𝜙𝑊& : %% (2a) 𝐅=ϕ 1 ℎ𝑢 ℎ𝑣 ℎ𝑢'+𝑔ℎ'/2 ℎ𝑢𝑣 ℎ𝑢𝑣 ℎ𝑣'+𝑔ℎ'/2 : ,𝐌=𝐈𝐝 (2b) where 𝜙 is the plan view fraction of space available to water. Porosity models are derived from the SWEs (1, 2a-b) by performing a volume averaging [9] over a control volume containing a water and a solid phase. The underlying assumption of the SP model is that the statistics of the water-solid partition converge to the same values when computed over a horizontal 2D domain and over its boundary. As noticed in [10], this is equivalent to assuming the existence of the Reference Elementary Volume (REV) [11], an assumption that does not hold in typical urban areas [3]. The DIP model [2] was introduced as an improved version of the Integral Porosity model [10]. It uses the following definitions for the terms in Equation (1) 𝐮=𝜙( + ℎ ℎ𝑢 ℎ𝑣 / ,𝐬= 1 𝑅−𝜙(𝐼 𝜙(𝑔ℎ 6 𝑆",$ −𝑆%,$ 8 +𝑔ℎ𝜕$ ( 𝜙)−𝜙( ) −𝜙(𝐼𝑢+𝜙(𝑊$ 𝜙(𝑔ℎ 6 𝑆",& −𝑆%,& 8 +𝑔ℎ𝜕& ( 𝜙)−𝜙( ) −𝜙(𝐼𝑣+𝜙(𝑊& : %% (3a) 𝐅= ⎣ ⎢ ⎢ ⎢ ⎡ 𝜙(ℎ𝑢 𝜙(ℎ𝑣 *! " *#ℎ𝑢'+𝜙(𝑔ℎ'/2 *! " *#ℎ𝑢𝑣 *! " *#ℎ𝑢𝑣 *! " *#ℎ𝑣'+𝜙(𝑔ℎ'/2 ⎦ ⎥ ⎥ ⎥ ⎤ ,𝐌= I 1 0 0 𝐃 M 𝐈𝐝 (3b) 𝐃=𝜖𝐑 P 𝜇+0 0 𝜇' R 𝐑,+ , 𝜖= S 1if%𝜕!ℎ>0 0if%𝜕!ℎ≤0 (3c) where 𝜙% and 𝜙& are respectively connectivity and storage porosities, 𝜇'3(𝑘 = 1,2) is a momentum dissipation coefficient (between 0 and 1) in the kth principal direction, and R is the rotation matrix that transforms the x-direction into the first principal direction of D. The purpose of the dissipation tensor D is to model the momentum dissipation arising from the multiple wave reflections against obstacles on the subgrid scale. This term is active only in the presence of positive waves, hence the 𝜖 switch ( 𝜖 = 03or3𝜖 = 1) in Eq. (3c). In contrast with a classical source term in Eq. (1), the inertia-modifying tensor M preserves the selfsimilar properties of the solution when Riemann problems are dealt with. This self-similarity property was identified in [3], verified in [2], [4] and confirmed in [12] by numerical experiments. So far, no model has been proposed for the dissipation coefficients 𝜇' as functions of the hydraulic conditions and urban geometry, and the momentum dissipation tensor D must be calibrated. For all three models (1), (2) and (3), arbitrary initial conditions may be provided. The standard (timedependent) boundary condition types handled by the models are the following: – prescribed unit discharge (inflowing or outflowing), – prescribed water depth, – prescribed free surface elevation, – prescribed Froude number. In all four cases, if the water is flowing into the domain, the flow velocity vector is assumed normal to the boundary. If the water is flowing out of the domain, the lateral component of the momentum is advected with the flow. An impervious boundary is considered as a particular case of a prescribed unit discharge boundary, with the normal discharge set to zero. 2.2 Numerical aspects Equation (1) is solved numerically using a finite volume approach on unstructured grids. Explicit, Godunovtype shock-capturing methods are used for the solution of the hyperbolic part. The algorithm is designed so as to allow for elements with an arbitrary number of edges. Therefore, the mesh may be a mixture of triangular, quadrangular or any type of polygonal elements. Besides, the cells are not required to be convex. The solution procedure involves a first-order time splitting [13]. Within a computational time step, the sequence is the same for all models: 1) Hyperbolic part 1.1) Determine the maximum permissible time step based on CFL requirements for solution stability on arbitrary-shaped grids [14]. 1.2) Reconstruct the flow variables in case a higher-order reconstruction is used (MUSCL-EVR scheme [15][14]). 1.3) Compute the fluxes and the source terms induced by porosity gradients and topography gradients across the interfaces (including boundary interfaces) between the computational cells using approximate, HLLC-type [15] [16] Riemann solvers. In order to speed up the computational process, the user is allowed to specify a threshold water depth. When the water in two neighbouring cells is smaller than this threshold, the two cells are assumed dry and zero fluxes are set automatically across the interface. The details of the Riemann solvers used can be found in the references for the SP and DIP models [1] [2]. 1.4) Carry out a mass and momentum balance over the non-dry cells. 1.5) Apply the momentum dissipation tensor D to the cells where the water level has been identified as rising as a result of Step 1.4. 1.6) Revise the balance using a divergence correction procedure in case negative water depths are obtained. Two options are available when the water depth becomes negative within a given cell : (i) all the fluxes across the interfaces of the cell under concern are multiplied by the same factor so as to obtain a zero water depth, or (ii) the computational time step is decreased to the maximum possible value that yields a zero water depth. 2) Compute the effect of source terms. All source terms are local functions of the average flow variables within a given computational cell, which involves solving only local, first-order ordinary differential equations with respect to time within each cell. The contributions of the various source terms are computed sequentially 2.1) Compute the effects of friction. The user can specify the friction model formulation (Manning, Strickler, Chezy, etc.) 2.2) Compute the effects of precipitation/infiltration 2.3) Compute the effects of wind shear stress It is worth noting that the MUSCL-EVR approach [14] used in the hyperbolic part is significantly computationally cheaper than the classical second-order time stepping MUSCL approach (also called MUSCLHancock [18]), because it involves a single step time integration and the reconstruction of a single variable (the free surface elevation) instead of three in the usual procedure. 3. SOFTWARE FEATURES SW2D is multi-platform and runs on Linux, MacOS and Windows. It includes all models presented in Section 2 above: classical Shallow Water models, SP and DIP porosity models (other porosity models are currently being implemented and will be available in next software releases). It can be used with any initial data and with the four boundary condition types mentioned above. Several types of forcing (precipitation, infiltration, wind, friction) are also implemented. 3.1 Input files SW2D requires input files providing geometry information, initial and boundary conditions, forcing terms, parameter values, etc. Here we briefly describe the files needed for a simulation (for more details, see the online documentation on the SW2D site): - Mesh file: SW2D is compatible with the 2DM format (ascii files) provided by the Surface-water Modeling System (SMS) Aquaveo(c) software. SMS can be downloaded from the Aquaveo website. Alternatively, users may generate 2DM file from your own data using Python modules such as Py2DM. It is also possible to convert 2D Selafin/Seraphin files (TELEMAC-MASCARET format) using a dedicated executable included in the SW2D distribution. The mesh is unstructured. The computational cells may have an arbitrary number of edges. It is for example possible to mix triangles with quadrangles, as long as these meshes are correctly defined in the 2DM file. The 2DM format incorporates the elevation of the mesh nodes, so that no additional elevation data file is needed. Input files are processed by the sw2dConverter program. This converter produces a binary file containing all the necessary information for the main SW2D program. The format of this binary file is specific to SW2D. As stressed above, sw2dConverter can read files of various formats, such as those from SMS, BlueKenue or Telemac-Mascaret. - Initial condition file: this text file specifies the initial state of the simulation. The well-posedness of the shallow water problem requires that the initial water level/depth and the components of the flow velocity/unit discharge vector be specified for all computational cells. Two options are made available to the user. If the initial state is uniform, only one line is needed (Figure 1). The values read within this line are applied uniformly to all cells. Otherwise, the initial values are specified on a cellby-cell basis, with one line containing the 3 initial state variables for each cell in the mesh. Figure 1: Content of an IC file for water at rest (surface elevation = 1m) - Boundary condition time series file: this file provides the time series for all the boundary conditions used in the simulation. The time series are placed after a header specifying the number n and the types of boundary conditions. They are arranged in the form of n + 1 columns. For each line, the first column is the time value (in seconds) and the remaining n + 1 records indicate the numerical value of the boundary condition at this specific time (Figure 2). During the simulation, the numerical value of the boundary condition is interpolated using the nearest two times in the file. When the simulated time is larger than the largest time in the time series file, the last read boundary condition numerical value is used. Figure 2: Content of a BC file with 3 types of boundary conditions, all constant in time: zero discharge, free surface at 8m and Froude number equals to 0.7. Changes in the last line would lead to time interpolation for q (resp. z,c) at boundary 1 (resp. 2, 3) between t=0 and t=86400 - Input SW2D file: this file specifies the simulation settings: choice of the physical model (shallow water, SP, DIP, etc.), physical and numerical settings (simulation time time, type of boundary conditions, CFL constraint, type of finite volume scheme) and output parameters (type and frequency of simulation result writing). These parameters are assigned a default value (see the SW2D documentation for more details) that may be overwritten by the user (Figure 3). a supprimé: Figure 1 a mis en forme : Non souligné a supprimé: Figure 3 Figure 3: Sample SW2D input file. The dual integral porosity (DIP) model is solved using the first-order Godunov scheme. The simulated time is 100 s, the maximum permissible computational time step (dtmax) is 1 second and the simulation results are stored in the form of maps every second. - [Optional]: depending on the physical models (porosity, friction, source terms, etc), additional files are needed for the specification of the corresponding model parameters. The format for such files is provided in the software documentation, along with sample files. 3.2. Output files The model outputs are stored every dtmap seconds. A unique (ascii) output file is created for each stored time. The simulation results are stored on a cell-by-cell basis, with one line for each cell (Figure 4). The verbosity of the storage is controlled in the SW2D input file by a flag that may take 3 values: - ‘none’: nothing is stored - ‘standard’: main variables are stored (x, y, zb, h, q, r) - ‘all’: all variables are stored (main variables above, plus with additional ones such as velocities, free surface elevation, Froude number, etc.) Figure 4: Sample output file (verbosity flag set to ‘standard’). 3.3 SW2D workflow The simulation workflow (Figure 5) consists of two main steps. The first step is geometry preprocessing. It consists in converting the mesh and boundary condition input files to the specific format required by the computational engine. The second step is the hydrodynamic computation itself. These two steps are distinct because the geometry processing step can be time-consuming and needs to be done only once. This allows the same model mesh and geometry to be used in numerous simulation scenarios without having to process the mesh again. SW2D is provided with a Graphical User Interface (GUI) that makes it possible to launch and visualize the simulation progress without the need for external tools (see Section 5). Experienced users may also choose to bypass the GUI and use the software in command line mode. In both cases, the output files are stored in a dedicated directory, at a frequency chosen by the user thanks to the dtmap parameter. During the program execution, a number of log messages can be displayed by the user in a dedicated window (or in the console): information, warnings and errors. a mis en forme : Non souligné a supprimé: Figure 4 a mis en forme : Non souligné a supprimé: Figure 5 Figure 5: Workflow for SW2D execution 4. LICENSING AND INSTALLING 4.1 Licensing SW2D is the property of Inria (Institut National de Recherche en Informatique et Automatique) and UM (Université de Montpellier). It has been developed in C++ by researchers of the LEMON project-team (common to Inria, UM, CNRS and IRD). Although it is still under development, it is already distributed under 3 different types of licences: Public Research License Inria and UM grant to the academic user, a free of charge, without right to sublicense, nonexclusive right to use the software for research purposes for a period of one year. Private Licence Inria and UM provide SW2D to any private user (including companies) within the framework of a contract between the parties. Educational Licence Inria and UM have formed partnerships with engineering schools where SW2D is distributed free of charge, strictly within the framework of the training of engineering students. If results obtained through the use of SW2D were to be published, authors should cite corresponding articles, taking information from the SW2D publication page. The software is provided only as a compiled library file. Any decompiling process is strictly unauthorized. Upon request (in the framework of research partnerships), any partner may become a member of the SW2D development team and be granted access to the SW2D C++ source code. 4.2 Installing SW2D This section presents the example of a binary installation file for Windows, similar to the one provided to educational partners. The process is similar to any Windows program installation (see Figure 6 and Figure 7). 1. Request the sw2d-install-windows.exe file from the SW2D consortium and download it. 2. Locate and double-click the downloaded .exe file. It will usually be in the ‘Downloads’ folder. 3. A dialog box will appear. Follow the instructions to install the software. a mis en forme : Non souligné a supprimé: Figure 6 4. The software will be installed. For more convenience, it is recommended to add shortcut icons on the user desktop. Figure 6: Installation process under windows: first steps Figure 7: Installation process under windows: last steps Sample files, which allow testing the software, are provided with the installation along with a user's guide. 5. TEST CASES The present section focuses on different test cases built for teaching. They exhibit the different features of SW2D as the principal points where the 2D flow modeler should pay attention to. 5.1 Gardon test case This test case aims to model a flood propagation in the Gardon river near Ners (downstream of Alès, France). It is part of the teaching curriculum at Polytech Montpellier. This site corresponds to the typical configuration where a 2D model should be used, with a sharp meander and a complex channel geometry inducing drying beds, flow divertion into parallel branches and highly non-uniform velocities across the channel section. The channel is also subjected to widening and narrowing, inducing significant transverse flow velocities. Hydraulic structures are also present: a strong narrowing of the floodplain near the downstream section due to a road embankment and a sill (with a 45 degrees angle with the main flow direction) with elevation variation of nearly 2m potentially leading to transcritical flow on the steep downstream slope. This test case is therefore rather synthetic of typical river hydraulics complexities, with mixed and competing influences of geometry, hydraulic parameters and boundary conditions. Figure 8 shows the test case situation and the mesh used in the 2D model. a mis en forme : Non souligné a supprimé: Figure 8 Figure 8: Gardon test case situation and mesh (basemap ©IGN http://geoportail.gouv.fr). The hydraulic model is the classical 2D shallow water model. The simulation input files are specified in Table 1. Table 1: Input files for the Gardon test case. File name Object friction_strickler_map.txt Distribution of the friction coefficient (uniform with 𝐾 = 50𝑚(/*. 𝑠+( ) Gardon.2dm Mesh file (5178 unstructured cells) Gardon.bc Boundary condition file hydro_boundary_time_series.tx t Boundary condition time series (prescribed hydrograph upstream, time-dependent Froude number downstream) initial_conditions_map.txt Initial conditions map (uniform water level and zero flow assigned uniformly to all cells) input.sw2d SW2D input file The propagation of a propagation is simulated with a prescribed upstream unit discharge rising from 0.1m2/s to 15m2/s. Figure 9 shows the results obtained at the peak time of the inflowing hydrograph. As shown in the screen capture, the GUI allows for a split view to visualise several variables at the same time. a mis en forme : Non souligné a supprimé: Table 11 a mis en forme : Non souligné a supprimé: Figure 9