Discrete Boltzmann Numerical Simulation of Simplified Urban Flooding Configurations Caused by Dam Break

The Discrete Boltzmann Equation (DBE) is a versatile simulation method, consisting of linear advection equations, which can be applied to the Shallow Water Equations (SWE). The aim of this study is to assess the accuracy of the DBE to simulate dam-break of shallow water flows in presence of obstacles. Dam-break flows in the presence of obstacles can be regarded as simplified models for floods in urban areas. In this study, three cases of dam-break flows in the presence of obstacles are considered: two with an isolated obstacle and the third in presence of an idealized city. A comparison between DBE and benchmark results shows a fair agreement, confirming the validity of the DBE as simulation method for the SWE.


INTRODUCTION
Flooding is considered one of the most significant hazards society is facing. Huge losses in terms of human lives, damages to buildings, houses, and civil infrastructures are caused every year by hurricanes and flooding. Moreover, the impact of flooding in the future is expected to become increasingly important due to deforestation and depopulation of rural zones (Sanders and Schubert, 2019).
Urban flooding are particularly fearsome events, considering that more than half of the world's population lives in urban areas and that the urban population is continuously increasing (Song et al., 2019). The prediction of expected consequences of urban flooding as well as the design and the adoption of mitigation measures is a crucial issue that nowadays takes advantages of numerical simulations. The main mathematical model for flood simulation is based on the two-dimensional Shallow Water Equations (SWE), which is obtained by averaging mass and momentum balance equations along the vertical direction under the classical assumption of hydrostatic pressure distribution (Valiani and Begnudelli, 2006). In the last decades, great attention has been paid to the development of numerical solvers for the SWE (Toro and Garcia-Navarro, 2007;Toro, 2009). The latter tackle directly the SWE by means of the finite volumes discretization and have reached a sufficient level of complexity so that they are able to account for crucial issues such as to cite just a few, the treatment of topography source terms (Duran and Marche, 2014;Hou et al., 2018), the use of unstructured grids (Zhao et al., 2019), and the evolution of wet-dry fronts (Ferrari et al., 2019).
An alternative option is represented by the Smoothed Particle Hydrodynamics (SPH) formulation of the SWE (Chang et al., 2018). The latter give rise to strongly non-linear equations whose numerical treatment needs particular care. Moreover, the numerical algorithms of SPH-based methods are intrinsically non-local, thus implying typically a high computational burden.
Solution methods for the SWE have also been developed in the Lattice Boltzmann Equation (LBE) framework (Zhou, 2004); these methods take advantage of the simplicity and versatility of the stream and collide LBE algorithm (Succi, 2001). The LBE proposes a mesoscopic representation of the flows based on the concept of the probability density function associated to a fluid particle traveling along a given direction with a given particle velocity (Succi, 2001). The finite set of particle velocities is defined in such a way that a symmetric regular spatial lattice is generated (Succi, 2001). Macroscopic quantities, such as vertically averaged velocity and water depth in shallow water flows, are calculated as statistical moments of the probability density functions (Zhou, 2004). However, no matter how appealing the application of the LBE to shallow water flows may be, the LBE is limited to the simulation of subcritical flows unless an ad hoc-defined lattice is adopted instead of the usual one, as in Hedjripour et al. (2016), where only one-dimensional transcritical shallow water flows have been considered. This is an unacceptable limitation, as the simulation of flooding problems cannot rule out the occurrence of transcritical flows. This fact addressed the research toward the use of multispeed particle velocity sets together with a finite difference discretization of the Boltzmann equation. The  Kleefsman et al. (2005). resulting model, the multispeed Discrete Boltzmann Equation (DBE), is able to deal with transcritical flows, both 1D and 2D, preserving the linearity of the LBE numerical algorithm (La Rocca et al., 2015). The DBE has been also successfully applied to the investigation of the dynamics of two-phase (La Rocca et al., 2018) and polydisperse (La Rocca et al., 2019) shallow granular flows.
The aim of this paper is to assess the ability of the DBE in simulating dam-break flows impacting on obstacles as a preliminary step for the simulation of urban flooding events on real topography. Indeed dam-break flows impacting on obstacles can be considered as simplified models of urban flooding, and the scientific literature abounds with experiments and numerical simulations of dam-break flows impacting on obstacles. In this paper, the experimental and numerical dam-break flows of Kleefsman et al. (2005), Soares-Frazão and Zech (2007), and Soares-Frazão and Zech (2008) are considered as benchmarks. In Kleefsman et al. (2005) and Soares-Frazão and Zech (2007), the dam-break flow interacts with an isolated obstacle, while in Soares-Frazão and Zech (2008) the dam-break flow interacts with an idealized city, which realized by means of a regular array of parallelepipeds.
In these works, water depth and velocity are measured at given measurement points and/or are numerically calculated by means of previously validated finite volume-based numerical solvers either applied to the SWE Zech, 2007, 2008), or to the Navier-Stokes equations (Kleefsman et al., 2005). A comparison between numerical results obtained solving the SWE by means of the DBE and reference results is encouraging and proves the ability of the DBE as a simulation tool of dam break flows interacting with obstacles.
The structure of the paper is as follows: first, the mathematical model is briefly presented; second, the considered benchmarks are described; third, results are reported and commented on; and fourth, conclusions are drawn.

THE MATHEMATICAL MODEL
The SWE, obtained by averaging mass and momentum balance equations along the vertical direction under the classical assumption of hydrostatic pressure distribution (Valiani and Begnudelli, 2006), assume the following form: where h is the water depth and U, V the vertically averaged velocity components along the horizontal directions x, y respectively. n m is the Manning friction coefficient. The first term at right hand side of the second and third Equation (1) accounts for the bottom slope and is calculated as in Valiani and Begnudelli (2006), i.e., taking as a constant the free surface elevation η, defined as the following: η = z f + h, being z f the bottom elevation. This formulation of the bottom slope term ensures the correct balance between flux gradients and source terms (Valiani and Begnudelli, 2006). The free surface η is considered instantaneously and locally constant: i.e., it is constant within the computational cell during the time step.
The DBE is given by (La Rocca et al., 2015): where f k is the probability density function relative to the kth particle velocity, whose cartesian components along x, y are respectively c xk , c yk . The particle velocity set as the kth equilibrium probability density function f eq k is defined in La Rocca et al. (2015) and reported concisely in the Appendix for the reader's convenience. τ * is the dimensional relaxation time, F e x , F e y are the cartesian components of the external forces acting on the flow. In this case, they are defined by the following: The DBE (2) is equivalent to the SWE (1) in the sense that the water depth h and the unit discharges Uh, Vh, obtained as the zero coincide with the water depth h and the unit discharges Uh, Vh that are obtained by solving the SWE (1) provided that the ratio Fr/Re of the Froude and the Reynolds number of the flow is small (La Rocca et al., 2015). It is worth observing that the DBE (2) that consists of a set of uncoupled, linear advection equations with constant advection velocity; this is the most appealing aspect of the DBE, when compared with the highly non-linear SWE, and it is an advantage per se. Indeed, the numerical solution of the DBE can be tackled by means of computational methods that are rather simple in comparison with those needed to treat the complex nonlinearities of the SWE. Such an advantage may involve a gain of efficiency, whose quantification is, however, very difficult, as it depends on a variety of factors, such as e.g., the level of coding optimization and parallelization, the performance of the processors, etc., and this is beyond the scope of this work.

Computational Details
The DBE (2) is discretized by means of an explicit first order approximation of the time derivative and a first order upwind approximation of the space derivative (La Rocca et al., 2015). Details of the numerical discretization, which is rather standard, can be found in Fletcher (2012) and are not reported here for the sake of concision. Boundary conditions are imposed on the macroscopic variables (water depth and vertically averaged velocity); then, the probability density functions at equilibrium are calculated and the values assigned as boundary values to the probability density functions (Ubertini et al., 2003). Stability is ensured by setting the order of magnitude of the Courant number C (C = U max × t/ s) to O(10 −1 ), U max , t, and s being the highest celerity of flow perturbations, the time step, and the space step, respectively. The adopted numerical method introduces numerical viscosity proportional to s, whose effect decreases by increasing the number of computational grid points. Finally, the value of the dimensionless relaxation time τ , defined as τ = τ * / t has been chosen in the range 0.5 < τ < 0.6.

Dam-Break Flow Impacting on an Isolated Obstacle-I
This case has been considered in Kleefsman et al. (2005) and many subsequent works, as it is widely used as benchmark to validate free-surface software (Kees et al., 2011). The reader can refer to Kees et al. (2011) for a detailed description of the By lifting the gate, the water column collapses under the action of gravity and impinges on a fixed parallelepiped (0.40 × 0.16 × 0.16 m) placed at 2.40 m from the right end of the tank, giving rise to a complex dam-break flow traveling back and forth within the tank. The Manning coefficient is estimated to the following: n m ≈ 0.01 s/m 1/3 . The experimental setup is equipped with four heights (H 1 , H 2 , H 3 , and H 4 ) and four pressure probes (P 1 − P 4 ) (Kees et al., 2011).
In Figure 1, the plots of the DBE numerical profiles of the free surface at t = 0.4s and t = 0.56 are compared to the corresponding experimental profiles of Kleefsman et al. (2005). At t = 0.4s the flow has almost reached the obstacle, while at t = 0.56 the impact has occurred, and the formation of a water splash with a height of around 0.4 m is evident both in the experiment and in the numerical simulation. The position of the height (H 1 , H 2 , H 3 , and H 4 ) and pressure (P 3 ) probes is also reported in Figure 1.
In Figure 2, the time history of the water height at probes H 1 , H 2 , H 3 , and H 4 is plotted. Dots represent the experimental results obtained at MARIN, while the solid line represents the DBE numerical results. The latter are those obtained with the finest grid (800 × 280). The space step s is equal to 4 × 10 −3 m, while the order of magnitude of the time step t is O(10 −4 )s. The general agreement is fairly good. The main discrepancy between experimental and numerical results is evident for the probe H2, where the DBE results overestimate the peak of water height due to the impact of the flow on the obstacle, and this is attributed to the fact that the shallow water model does not describe correctly the strong vertical dynamics involved during impact; the duration is, however, very short in comparison to the duration of the considered experiment. However, it is interesting to observe that the instant of impact of the flow on the obstacle (t ≈ 0.5s) is correctly predicted by the DBE numerical results. The experimental and numerical time history of the pressure measured by the probe P 3 is shown in Figure 3. In the framework of the shallow water equations, the numerical pressure is calculated assuming a hydrostatic pressure distribution: p = ρg(h − z), z being the elevation of the pressure gauge from the bottom of the tank. The experimental pressure shows an isolated peak, due to the sudden impact of the water against the obstacle, which is a typical dynamic feature and of course cannot be reproduced by the numerical hydrostatic pressure distribution. However, as time goes by, dynamic effect loose importance and the hydrostatic pressure distribution becomes a fair representation of the pressure field.

Dam-Break Flow Impacting on an Isolated Obstacle-II
This case has been considered in Soares-Frazão and Zech (2007), where the reader can find a detailed description of the experimental setup. The experiments have been realized in a tank 36 m long, 3.6 m wide, with a trapezoidal cross section. A gate (width 1 m), placed 7.70 m from the left end of the tank, separates the reservoir from the part of the tank representing a valley with an isolated building. The latter is a parallelepipedal obstacle, with a 0.8 × 0.4 m base area. The 0.8 m long side is inclined at an angle of 64 • with respect to the y axis (Soares-Frazão and Zech, 2007). The water level is set to 0.4 m and 0.02 m within the reservoir and the rest of the tank, respectively.
Finally, the Manning coefficient is set to n m ≈ 0.01 s/m 1/3 as in Soares-Frazão and Zech (2007). Water height and velocity measurements have been performed at several points (Soares-Frazão and Zech, 2007). Velocity values are measured at 0.036 m from the bottom and cannot be compared with the DBE numerical vertically averaged velocity, as the flow presents intense re-circulation zones with high values of the vertical velocity component. In Figure 4 are plotted the time histories of the free surface at measurement points G 1 , G 3 , G 4 , and G 6 . Solid lines represent experimental values, while dotted and dashed lines represent the DBE numerical values obtained with a coarse (600 × 67) and a fine (1500 × 165) grid, respectively. The difference between the DBE numerical results obtained with the two grids is not meaningful, mainly at measurement point G 6 . The differences with the experimental results observed at measurement points G 1 , G 3 , G 4 , and G 6 are ascribed to the fact that there the vertical motions are not negligible and that the shallow water model loses consistency. Nevertheless, the DBE numerical results grab the main features of the flow, and the general agreement with experiments can be considered fairly good. The DBE numerical results are in satisfying agreement with the numerical results of Noël et al. (2003), reported in Soares-Frazão and Zech (2007) and obtained by directly integrating the SWE. This fact is depicted in Figure 5 where the snapshots of the free surface at t = 1s, t = 3s, and t = 10s reproduce the same flow features observed in the numerical calculation of Noël et al. (2003).

Dam-Break Flow Impacting on an Idealized City
This case is reported in Soares-Frazão and Zech (2008), where the experimental setup is described with great detail, and has been performed with the same experimental setup tank used in Soares-Frazão and Zech (2007). The idealized city has been realized by means of 25 impervious wooden blocks of 0.30 m × 0.30 m, representing buildings, separated by streets 0.10 m wide and arranged in a regular array of 5 × 5. According to Soares-Frazão and Zech (2008), the ratio of the street to building width is realistic despite the strong simplification of the representation of the city.
A gate (width 1 m), placed 7.70 m from the left end of the tank, separates the reservoir from the part of the tank containing the idealized city. The water level is set to 0.4 m and 0.011 m within the reservoir and the rest of the tank, respectively. Finally, the Manning coefficient is set to n m ≈ 0.01 s/m 1/3 as in Soares-Frazão and Zech (2008). Water height has been measured by means of several resistive level gauges, placed along straight lines along the streets of the city, while velocity has been measured both on the free surface, by means of a digitalimaging technique, and close to the bottom by means of an acoustic doppler velocimeter (ADV) (Soares-Frazão and Zech, 2008).
In Figure 6, the DBE numerical results and the experimental and numerical results of Soares-Frazão and Zech (2008) relative to the free surface profile at 4s, 5s, and 6s are shown. The free surface profile is plotted along the third longitudinal street. The agreement is satisfying. Minor discrepancies at t = 5s can be attributed to numerical oscillations in the central part of the obstacles' array. In Figure 7, the DBE numerical Froude number profile at t = 4s along the third longitudinal street, as indicated by the sketch, is shown. The high values attained by the Froude number allow to state that the DBE is a robust tool to investigate super-and transcritical flows. It is worth observing the fairly good agreement among the DBE numerical and experimental velocity values and that the position of the main hydraulic jump (x ∼ 4.3 ÷ 4.4 m), shown in Figure 6, matches with the sudden change from Fr > 1 to Fr < 1 of the Froude number shown in Figure 7. Finally, in Figure 8 the computed water depth in the neighborhood of the city area is shown at t = 5s, t = 10s, respectively. The complex features of the flow comply with those shown in Soares-Frazão and Zech (2008).
In particular, at t = 5s (Figures 8A,B), the flow between the gate and the city is mainly supercritical; it forms a sort of hat-like shaped subcritical zone downstream of a hydraulic jump which front is placed at x H = 3.95 m along the centerline of the tank. From the edges of the subcritical zone, at x E = 4 m, y EN = 0.7 m, and x E = 4 m, y ES = −0.7 m, two symmetrical oblique waves detach and end at the wall of the tank (Figures 8A,B). Points x H , x E , y EN , and y ES coincide in Figures 8A,B. Within the longitudinal street of the city the flow is mainly supercritical, as in Soares-Frazão and Zech (2008). Finally, the wake zone immediately downstream of the city, where a triangular wave is observable, is qualitatively similar to that observed in Soares-Frazão and Zech (2008).
At t = 10s (Figures 8C,D), the structure of the flow is similar to that observed at t = 5s. The subrictical zone is almost rectangular, the front being located at

CONCLUDING REMARKS
The multispeed Discrete Boltzmann Equation (DBE) is a method recently developed to simulate transcritical shallow water flows, and the main advantage is the simplicity of the mathematical model: a set of linear, uncoupled, purely advective equations. In this paper, the DBE has been assessed as solution method for the Shallow Water Equation (SWE) when dam-break flows impacting on obstacles are considered. These flows are strongly unsteady and present frequent sub-supercritical and supersubcritical (hydraulic jumps) transitions.
The consideration of dam-break flows impacting on obstacles is important as they represent simplified models of urban flooding, thus constituting a benchmark for numerical methods, thanks to the wide variety of reference numerical and experimental data published in literature.
In this paper, three dam-break flows impacting obstacles have been considered: two on an isolated obstacle and one on an idealized city that has realized by means of a regular array of parallelepipeds.
The experimental configurations were realized within tanks with a smooth, though sloping, bottom. The latter has been treated in order to account for the correct balance between flux gradients and source terms.
The comparison between DBE and reference results (both numerical and experimental) shows that the DBE can be considered as a reliable solution method for dam break flows impacting on obstacles, even in the presence of complex flow configurations, such as e.g., the case of the idealized city. Discrepancies between DBE and reference results can be attributed both to the use of coarse grids and to the accuracy of the numerical method and to the occurrence of violent vertical motions, which undermine the consistency of the shallow water model. It is worth observing that, for what regards the specific performance of the proposed models in urban areas, there are no restrictions for the DBE on the attainable grid resolutions. Any grid convergence rate can be achieved by just adopting a higher order numerical scheme. Future works will be addressed to the consideration of realistic topography as well as to the implementation of more efficient numerical methods, such as e.g., the finite volume method with unstructured, adaptive grids.

DATA AVAILABILITY STATEMENT
The datasets generated for this study are available on request to the corresponding author.

AUTHOR CONTRIBUTIONS
ML has formulated the Discrete Boltzmann Equation. PP performed the numerical simulations. SM analyzed the results and prepared graphics. All the authors have discussed together about the outcomes and their interpretation. All authors contributed to the article and approved the submitted version.

ACKNOWLEDGMENTS
Valuable discussions with Dr. A. Montessori are kindly acknowledged.