Integral porosity shallow water model at district scale – Case study in Nice

. After three hours of intense rainfall, the city of Nice was ﬂash ﬂooded on October 3, 2015, resulting in casualties and severe damages in property. This study presents a porous shallow water-model based numerical simulation of the ﬂash ﬂood event in a district of Nice, and compares the results with a high-resolution conventional shallow water model. This contribution aims to discuss practical aspects of applying a porous shallow water model to a real world case. The porous shallow water model is an integral porosity-type shallow water model. It uses unstructured triangular meshes. The conventional shallow water model is a distributed memory parallelized high-performance computing code, that uses a uniform Cartesian grid. The study site is an approximately 5km 2 spanning district of the city of Nice, France. Topography information is available in a 1m resolution and in addition, the available digital elevation model includes inframetric structures such as walls and small bridges. In the presentation of the case study, challenges of the pre-processing step of the integral porosity shallow water model are addressed. Notably, a method to semi-automatically generate “good” triangular meshes using the open-source geoinformation system QGIS and the mesh generator Gmsh is presented. During the post-processing step, the results of the porous model are mapped back onto the high-resolution topography to make the results more meaningful. The agreement between the high-resolution reference solution and the porous model results are poor. A speed up of about 10 to 15 was observed for the present case.


Introduction
Shallow water models with porosity aim to reduce the computational cost of urban flood simulations by moving the model point of view from the classical mesoscale to the macroscale. The motivation for this is that a macroscopic point of view enables the use of coarser resolutions. Cell sizes typically exceed the scale of the buildings. The resulting mesh has a lower number of cells and the computational effort is reduced. This is achieved by using the geometric information of the urban district to calculate porosity terms in each cell, that represent the fraction of the cell available to the flow. Pioneering work in this field is the single porosity shallow water model presented in [1,2]. An alternative shallow water model with additional porosity terms at the cell edges is presented in [3]. The porosity terms at the edge represent the fraction of the edge available to flow. This model is called the integral porosity shallow water model (IP). Research on the IP model is published in [4][5][6][7][8][9][10][11][12].
This contribution applies the IP model to simulate the pluvial flood event of October 3, 2015 in Nice, France. The novelty of this contribution is its application scale, as the IP model is not often applied at a district scale. An exception is a fluvial flood case studied in [10]. However, the domain in this fluvial flood example is smaller than the domain of the present case study. Test cases in literature often comprise analytical benchmarks for verification purposes [1,3,6], and laboratory experiments with rather simple urban layouts [1,7]. Most of these experiments are from the IMPACT project [13], and feature a dam-break flow in simplified urban districts.
Different challenges arise when the IP model is applied to a real world case at a large scale. The urban layout becomes complex, the terrain becomes rather rough. As a consequence, we identify two challenges that need to be addressed: (1) designing the computational mesh to optimally account for the buildings becomes non-trivial, and (2) interpreting the results obtained at the macroscale requires more effort.
In this contribution, we report pre-and post-processing methods that we applied to address these issues. Significant amount of pre-processing is required to collect the data required by the IP model and to design the computational mesh. We put effort into automating these steps. When the IP model results have been obtained at the macroscale, mapping the macroscale results back to the high-resolution topography data set recovers some of the flow characteristics. The interpretation of the results becomes easier.
In the following sections, we briefly show the governing equations of the IP model (section 2) and then discuss the case study of the flood event in Nice, France (section 3). Finally, we draw conclusions (section 4).

Integral porosity shallow water model
We briefly present the governing equations. The IP model equations can be written as where q is the vector of conserved variables, Ω is the area of the cell, F is the flux tensor, Γ is the counter-clockwise path along the boundary of Ω, n is the outwards pointing unit normal vector along Γ, s Ω is the source term inside the cell, s Γ is the source term projected to the cell edge, φ is the volumetric porosity and ψ is the conveyance porosity. The notation ∂ t indicates a temporal derivative and ∂ x and ∂ y indicate a spatial derivative in the Cartesian coordinate system in x-and y-direction, respectively.
The vectors are defined as where h is the water depth, q x and q y are the unit discharges in x-and y-direction, respectively, g is the gravitational acceleration, s 0,x = −gh∂ x z and s 0,y = −gh∂ y z, with z being the bottom elevation, are the bottom slope source terms in x-and y-direction, respectively, s f,x and s f,y are the bed friction source terms in x-and y-direction, respectively, s d,x and s d,y are the building drag dissipation source terms in x-and y-direction, respectively, and h 0 is the water depth evaluated for a quiescent state inside the control volume. i is a mass source term, e.g. rainfall intensity. Bed friction can be calculated by using one of the several existing friction laws, e.g. Manning's law with n being Manning's coefficient, and ||q m || being the Euclidian norm of the vector q m . The building drag dissipation source term can be calculated via a drag force equation, cf. e.g. [11]: Here, c d, j is a drag coefficient of proportionality that can be calculated in several ways. It usually depends on the area of the building that is facing the flow and the flow depth, cf. [3].
In this work, we calibrate c d,x = c d,y directly.

Pluvial flooding in a district in Nice, France
On October 3, 2015, an intense rainfall with a duration of about three hours fell on the French Riviera. The rainfall caused several district of the city Nice to be flash flooded, resulting in casualties and property damage. We study a heavy rainfall event in a district in the heavily urbanized city of Nice, that spans about 5 km 2 . The urban topography was provided by the municipality of Nice in form of a high-resolution digital elevation model (DEM), that has a spacing of 1 m. The DEM was further modified to include inframetric information such as vertical walls and other structures that are not captured by the resolution of the DEM [14]. The available DEM is visualized in Fig. 1. All building heights have been artificially increased to 190 m, so they exceed any terrain elevation signifiantly. The buildings are indicated with white color in Fig. 1 (right).

Pre-processing
We used the geoinformation system QGIS [15] and the meshing software Gmsh [16]. In addition, several self-written small scripts have been used to process data. Using QGIS, we extracted a two-dimensional polygon layer of the buildings inside the domain. This layer was used later to compute the porosity terms of the model and to guide the meshing process. Further, using QGIS we have extracted the boundary of the district as a polygon layer. This layer was used to build the mesh. The extracted boundary polygon contained too many points and meshing this domain resulted in very small cell sizes near the boundary. Thus, we simplified the boundary polygon by manually removing points from it.
We further extracted the vertices of the simplified polygon using the Gmsh-plugin of QGIS. This allowed an instant formulation of points using Gmsh syntax. In order to better represent the blocking effects of the buildings, we extracted the building centroids from the building layer and forced them to be included as vertices in the mesh. The process is illustrated in Fig. 2.
The bottom elevation and its gradient have a large influence on the rainfall-runoff. This is important to note, because due to the coarse resolution of the computational mesh, several data points are located inside one cell. Further, the accuracy of the model degrades not only because the effect of buildings is not accounted for, but also because the bottom elevation is not sufficiently resolved and wrong gradients are computed. The presence of buildings in a cell also poses a difficulty, because we are only interested in elevations that are not buildings. Indeed, in our modeling approach, building information is considered via the porosity terms, and thus the elevation of the buildings is not required. We apply an algorithm that averages the bottom elevation of the vertices of a cell that are not located in a building.
Finally, we compute the porosity parameters using equidistant points in each cell using a recursive search. At each point, we check if it coincides with a building or not, and then calculate the porosity based on this evaluation. The resulting cell porosity is shown in Fig. 3.
This concludes our overview of the applied pre-processing methods. The outcome of this pre-processing is a computational mesh that includes the building centroids as vertices, a bottom elevation data set that does not include buildings and porosity parameters in each cell.

Model results
Rainfall is applied in the model as a mass source term according to a time series. The model results are observed at six measurement points. These measurement points have been selected from a larger set of about 200 points that coincide with real points of interest such as sensitive All boundaries of the domain are closed boundaries. This is reasonable because the boundaries coincide with real obstructions to the flow. The IP model results are compared with a high-resolution model (FullSWOF [17]). This model uses a mesh with square cells with an edge length of 1 m. The high-resolution simulation runs in parallel on 128 cores. The computational cost expressed as wall time is about 16 h. In contrast, the IP model runs on 32 cores and takes about 1.5 h to complete. We can calculate a rough estimation of the speed up as ξ = n CPU,1 n CPU,2 where ξ is the speed up, n CPU,1 and n CPU,2 are the number of CPUs utilized for the highresolution and porous shallow water model run, respectively, and t 1 and t 2 are the wall times for the high-resolution and porous shallow water model run, respectively. The comparison is plotted in Fig. 4. For the IP model, a peak of 3.5 · 10 −3 m is observed at about 70 min at point 1. This is approximately twice the water depth predicted by the high-resolution model. At points 2 and 3, the high-resolution model does not predict any water depths. However, the IP model predicts maximum water depths of 4.5 · 10 −3 m and 19.5 · 10 −3 m for point 2 and point 3, respectively. At point 5, both model results at least resemble each other for the first 55 min of the simulation, both reaching a water depth of about 5 · 10 −3 m. After this point, the model results diverge from each other. The highresolution model predicts a maximum water depth of 0.1 m at about 105 min. After that, the water depth decreases until it reaches 0. The IP model predicts a maximum water depth of 0.05 m at the end of the simulation. At point 6, both model results deviate from each other. As at point 1, the IP model overpredicts the water depths at point 6. The high-resolution model predicts a maximum water depth of about 0.05 m at about 85 min. In contrast, the IP model predicts the maximum water depth to be 0.15 m at the same time.
In addition to point measurements, we compare flood extension areas. We map the maximum water depth in each cell. A direct comparison is given in Fig. 5. We make following observations: (1) the north-western and south-eastern regions in the high-resolution model have a higher water depth than most areas, which was well reproduced by the IP model; (2) smaller areas and spots with high water depths like the south-west and the central area as well  as the north-east are not well replicated by the IP model. Overall, narrow spaces and streets can not be well accounted for in the IP model.

Post-processing
We map the maximum water depths that have been computed by the IP model to the highresolution DEM. This was done by suming up the IP model bottom elevation and its max-imum water depths. From this value the bottom elevation of the high resolution DEM was substracted. This way small and narrow spaces like streets can be incorporated into the visual presentation of the results, cf. Fig. 6.
In Fig. 6 the maximum water depth distribution in the IP model is more similar to the one in the high-resolution model. By mapping the water depth on the high resolution DEM, small local water depths are better represented. For instance the bridges in the south and south-west of the district now display a higher depth than in Fig. 5.

Conclusions
We have presented a case study of a real world pluvial flood event in an urban district. In our discussion, we focused on pre-processing and post-processing methods that we used to prepare and interpret the data. In summary, in order to use the IP model for flood simulations, it is necessary to calculate the porosity terms based on the urban geometry. This can be achieved by extracting a polygon layer that describes the buildings by means of image analysis. Any common GIS should be able to do this.
Following [3], we found in preliminary studies, that meshes which incorporate the building centroids as vertices in the computational mesh yield the best agreement between IP and high-resolution model results. In order to automate this process, we extract the centroids from the building polygon layer. Many mesh generators support a constrained meshing, e.g. Gmsh [16] and Triangle [18].
After the model run, we map the predicted maximum water depths on the high-resolution DEM. This makes the model results more meaningful and helps understanding and communicating the results. We must say that the IP model results do not agree with the high-resolution prediction. Reasons for the large deviation are that the coarse resolution does not represent the topography of the catchment well. However, rainfall-runoff depends on the correct representation of bottom slope. Further, rainfall-runoff is dominated by small water depths, usually smaller than 10 cm. This is the range of accuracy of the IP model.
Unfortunately, the presented case and its preliminary studies have been disheartening. It must be concluded that the IP model in its current state is not very suitable for rainfall-induced flood simulations. An accurate representation of the topography and a model accuracy in the range of a few millimeters can not be provided by the IP model. Therefore, we must recommend applying the IP model to flood events that do not depend on rainfall concentration in the urban area. Fluvial flood events might be a more suitable application area.
As an outlook, we list two points that might be worth future investigation: (1) Improving porosity term calculation by taking into account orientation of obstacles; (2) an even more refined meshing that forces more critical points in the triangulation. The drawback of the first point is that it requires more pre-processing. A limitation of the second point is that it might yield severely distorted meshes.