Modelling urban ﬂoods using a ﬁnite element staggered scheme with porosity and anisotropic resistance

. Artiﬁcial porosity models for urban ﬂooding use porosity as a statistical descriptor for the presence of buildings, which are then treated as subgrid-scale features. Computational e ﬃ ciency makes porosity models attractive for large-scale applications. These models are typically implemented in the framework of two-dimensional (2D) ﬁnite volume collocated schemes. The most e ﬀ ective schemes, falling under the category of Integral Porosity models, allow accounting for a wealth of sub-grid processes, but they are known to su ﬀ er from oversensitivity to mesh design in the case of anisotropic porosity ﬁelds. In the present exploratory study, a dual porosity approach is implemented into a staggered ﬁnite element numerical model. The free surface elevation is deﬁned at grid nodes, where continuity equation is solved; ﬂuxes are conveyed by triangular cells, which act as 2D-links between adjacent grid nodes. The presence of building is modelled using an isotropic porosity in the continuity equation to account for the reduced water storage, and an anisotropic conveyance porosity in the momentum equations to compute bottom shear stress. Both porosities are deﬁned on an element-by-element basis, thus avoiding mesh-dependency. Although su ﬀ ering a number of limitations, the model shows promising results.


Introduction
Artificial porosity models for urban flooding, inspired by the concept of porous medium, use the porosity as a statistical descriptor for the presence of buildings, which are then treated as subgrid-scale features. Such models are attracting increasing attention in the field of largescale simulation of urban floods.
Since the first appearance of porosity models more than twenty years ago [1], different formulations have been proposed in order to capture salient features of the complex flow fields occurring in urban areas using coarse grids. The Single Porosity (SP) approach, expressed in differential form, has been the subject of numerous works pointing out skills but also limitations [2][3][4][5][6][7][8]. The Multiple Porosity (MP) model proposed by Guinot [9], in which regions characterized by different porosity, flow depth and velocity coexist within the same control volume, accounts for head losses due to momentum exchanges and for anisotropic effects. The Integral Porosity (IP) model introduced in [10] uses an integral formulation to properly account for discontinuities in the urban medium and distinguishes between a storage and a connectivity porosity that represent, respectively, the volume fraction available for mass and momentum storage and the areal fraction of cell-sides acting on the computation of the inter-cell fluxes. IP model was further improved with the introduction of a depth-dependent porosity [11]. Finally, the Dual Integral Porosity (DIP) model [12] definitely formalizes the distinction between porosity and flux variables, with the former referring to control volumes, and the latter to boundaries; simple closure relations relate these variables to each other. Laboratory and numerical experiments [13] showed that the IP approach outperforms the SP one. Wave propagation properties are further improved in the DIP model.
The interest in porosity models is primarily due to the great advantage in terms of computational time and resources with respect to high-resolution applications based on the classical Shallow Water Equations (SWEs). Besides, the latest contributions on this topic [14,15] are devoted to shed light on intriguing issues inherent in such models. Specifically, the IP and DIP approaches are known to be markedly mesh-dependent. As already observed in [10], IP approach is inherently mesh-dependent as the conveyance porosity is defined locally at cellsides. While the use of suitable unstructured grids can limit this problem [10,16], Guinot [14] showed that oversensitivity to grid design is due both to a polarization of fluxes and to cell edge orientation that emerge when the flux porosity field is anisotropic. Mesh design guidelines are provided, but their fulfilment is not trivial in real-world applications [14].
Porosity models for urban floods have been typically implemented in the framework of finite volume, collocated schemes (i.e., conserved variables refer to the barycentre of the cell and fluxes are evaluated at the cell sides). In the present exploratory study, a staggered finite element (FE) numerical model is used [17,18]. The model was recently enhanced to account for anisotropic resistances and applied to model overland and shallow inundation flows in agricultural landscapes, where small-scale linear features (e.g., ditches) force anisotropy in the flow [19]. The scheme was shown to be very slightly mesh-dependent.
In order to apply this model to large-scale urban floods, a simple method is proposed in which, rather than a drag coefficient [8], the bottom shear stress accounts for anisotropy as a function of the cross-sectional area actually available to conveyance in a given direction. This area, suitably normalized, closely resembles the conveyance porosity introduced in the IP/DIP approaches. However, in this case, the conveyance porosity is not defined locally at the cell interfaces [10,12,14], but averaged over each cell [20], and it is defined tensorially along principal, mutually orthogonal directions.
The numerical model is first described (Section 2) and then applied to schematic case studies (Section 3) in order to highlight both strengths and limitations of the proposed approach, which are finally summarized in Section 4.

The numerical model
The dual porosity scheme for urban floods is implemented in the 2DEF finite element hydrodynamic model [17,18,[21][22][23] . The 2DEF model solves a modified version of the SWEs, which is obtained through a phase averaging procedure over a Representative Elementary Area (REA) in order to deal with flooding and drying over irregular topography [2], reading where h is the free-surface elevation, Y is the equivalent water depth defined as the volume of water per unit area (thus accounting for porosity), q = (q x , q y ) is the discharge per unit width, τ = (τ x , τ y ) is the bottom shear stress, ρ is the water density, g is gravity, Re are the horizontal components of Reynolds and dispersion stresses. The term ϑ(h) in Eq. (3) is an h-dependent storativity coefficient [2,22], defined as the ratio between the wet and the total area of a cell. In Eqs. (1) and (2), local and advective accelerations are lumped into the total time derivative of the depth-averaged velocity components q x /Y and q y /Y [17,24,25]. Unlike other well-known schemes that use staggered grids [e.g. [26][27][28], in the present scheme the free surface elevation is defined at grid nodes; triangular cells link together adjacent nodes and convey fluxes as a function of quantities defined on an element-by-element basis (free-surface gradient, stress terms, etc.). After being conveniently linearised, the SWEs are solved using an Eulerian, semi-implicit, Galerkin FE scheme for continuity; convection is treated explicitly by means of a Lagrangian approach [17]. The resulting staggered scheme naturally handles discontinuities in bottom elevations between adjacent elements, but is not suitable to deal with supercritical flows and shock waves, nor with rapidly varying flows [29][30][31].

Porosity treatment of urbanized areas
It is worth noting that the model equations (1)-(3) naturally handle the effects due to the presence of buildings in both storage reduction (by means of the parameter ϑ) and momentum equations (through the equivalent water depth Y, which depends on ϑ as well).
Defina [2] enforced a statistically-based model for subgrid macro-roughness by setting where erf(·) is the error function, z b the mean bottom elevation, and a r the macro-roughness height. In [22], with the aim of modelling hydrological response of agricultural catchments, ϑ(h) was modified to account for the porosity of topsoil layer, n e , as ϑ(h) = n e + (1 − n e ) · η(h); n e acts as the lower limit for ϑ (Fig. 1). Similarly, for a urban layout where the fraction of space available to water storage (i.e., the fraction of area not occupied by buildings) is the porosity φ (of course, φ ≥ n e ), the storage parameter in the continuity equation becomes where φ is now the upper limit of ϑ(h) (solid line in Fig. 1). The single porosity approach is isotropic by definition; yet, many urban layouts are characterized by markedly anisotropic patterns [10,32]. Moreover, a fraction of the domain that is free of buildings does not contribute to transport due to the sheltering effects of buildings that results in dead zones. The SP approach can thus be enhanced with the introduction of a conveyance (or connectivity) porosity, Ψ , that scales the transport capacity, as it was done in the IP and DIP approaches [10,12]. Of course, to account for dead zones, the conveyance porosity should be taken such that Ψ ≤ φ; the difference φ -Ψ thus represents an estimate of the fraction of the domain that does not contribute to transport [10]. Importantly, the conveyance porosity must account for directional dependence.
The usefulness of distinguishing between storage and conveyance porosities was clearly demonstrated in many literature studies [10,12,13,16]. However, IP and DIP schemes were implemented using collocated FV schemes, in which conserved variables (defined on an element-by-element basis) are "duplicated" at cell interfaces in order to compute fluxes; accordingly, while storage porosity used to define the conserved variables refers to grid cells [2,12], the conveyance porosity used to define fluxes only refers to cell interfaces [10,12]. Unfortunately, the point-based assessment of the connectivity porosity makes these models suffering from over-sensitivity to mesh design [14].
In the FE staggered model used in this study, while continuity is balanced at grid nodes (storage porosity at nodes is weighted-averaged from that of surrounding elements), fluxes are defined on an element-by-element basis. Accordingly, conveyance porosity must refer to the entire cell and not to cell interfaces. This need can be seen as an opportunity to better describe the conveyance porosity field by better reflecting the connectivity properties of the urban medium within the cells and not only at the cell interfaces [14]. The estimation of an element-by-element conveyance porosity is in general not trivial for complex urban layouts, and it is far beyond the scope of the present study. Here, the conveyance porosity is defined and introduced in the model under strong simplifying hypotheses. With reference to the one-dimensional (1D) schematic case depicted in Fig. 2a, Q denotes the total discharge; q = Q/W and U = q/Y are, respectively, the discharge per unit width and the velocity, phase-averaged over the area of the cell [2]. The presence of buildings results in a channel contraction, where the actual velocity is In Eq. (6), Ψ = w/W is the conveyance porosity, which is assumed equal to the width ratio of the channel contraction [33].
Similarly to the previous 1D case, in the two-dimensional (2D) schematic case of Fig. 2b, conveyance porosity can be estimated along two principal, mutually orthogonal directions (say L − T ), for which the conveyance porosity is respectively maximum (Ψ L ) and minimum (Ψ T ). Conveyance porosities are evaluated, along each principal direction, with reference to the narrowest cross-section (Fig. 2b); this allows to write the associated flow velocities as In the model, conveyance porosities are introduced only to evaluate the energy dissipation due to bottom shear stress, which is thus computed with reference to the velocities (u c L , u c T ). The use of the velocity at the contraction is equivalent to assume that water flows within the narrowest cross-section for the length of the entire computational cell. This is a simplifying hypothesis justified by the fact that energy dissipation varies quadratically with the velocity and linearly with the length of the path; most of the dissipation occurs where velocity is greater, and the jet developing downstream of the contraction causes the velocity to decrease slowly along the flow direction (i.e., the effective length of the narrowest section is longest than its purely geometrical length).
As the conveyance porosity is directionally dependent, resistance terms τ = (τ x , τ y ) are introduced in the model in tensorial form [19] τ ρgY = Rq (8) with q the modulus of the phase-averaged discharge per unit width and R a second-order tensor. By denoting with n L and n T the Manning roughness coefficients along, respectively, the L and T directions, the principal components of the tensor R are where H is a parameter derived from Y/φ according to Appendix A of [2], in order to account for the effect of sub-grid scale macro-roughness on bottom shear stress. If the principal directions are rotated by some angle α relative to the model frame (Fig. 2b), the components of R descend from the transformation law of a second-order tensor From here on, the model implementation follows exactly that presented in [19]. It's worth noting that the assumption of mutually orthogonal principal direction can be removed by considering different angles for each principal axes [12]; nonetheless, the effectiveness of the tensorial representation of conveyance porosity, the existence of (and the possibility of detecting) principal directions, and the computation of hydraulically significant conveyance porosities in real urban layouts are questionable tasks that deserve further investigation. At the moment, the model does not account for building drag (this can be included according to [12]), nor for transient momentum dissipation due to shock wave reflection [12,34], since the model is not intended to deal with supercritical flows nor with shock waves.

Model results
The dual porosity model presented above is tested considering different schematic applications and comparing the results from dual porosity models with those of a "reference simulation", i.e., a refined solution in which the buildings are explicitly resolved using the "building hole" technique [16]. As a rule, the porosity models are run over a mesh that is much coarser than that of the reference simulation, thus leading to sensible saving of computational resources and time. Additional simulation, not shown in this contribution, confirmed that numerical solutions are insensible to the mesh design, apart from the resolution of numerical solutions which obviously depends on the resolution of the mesh.
In the first application (Test 1), a schematic urban layout (Fig. 3) is placed at the centre of a tilted plane (bottom slope of 0.065% southward, Manning roughness coefficients n L = n T = 0.0286 m −1/3 s, macro-roughness height a r = 0.1 m). Based on the geometry of the layout, the urban porous medium is characterized by a storage porosity φ = 0.3 and, following Fig. 2b, by conveyance porosities Ψ L = 0.25 and Ψ T = 0.067, with α = 50 • . A constant inflow discharge of 50 m 3 /s is prescribed at the centre of the northern boundary, and a suitable rating curve is imposed at the southern boundary. Simulations are protracted until steady state conditions are attained.
The model results are shown in Fig. 3. The solution from the dual porosity model (central panels) closely resembles that obtained from the reference simulation (left panels). On the contrary, the SP approach (right panels) is unable to correctly characterize the marked anisotropy in the flow field forced by the building layout. In the second application (Test 2), a hypothetical flood wave (Fig. 4) is routed over an initially dry, tilted area, with the street grid of a real urban layout placed at the centre of the plane (Fig. 5, bottom slope of 0.065% southward, Manning roughness coefficients n L = n T = 0.0154 m −1/3 s along the streets and n L = n T = 0.0286 m −1/3 s elsewhere, macro-roughness height a r = 0.1 m). The storage and the conveyance porosities are estimated geometrically as φ = 0.22, Ψ L = 0.14, and Ψ T = 0.08; principal axes are oriented at α = 65 • to the model frame. The discharge hydrograph (Fig. 4) is prescribed at the centre of the northern boundary, and a suitable rating curve is imposed at the southern boundary.
The results from the dual porosity model (lower panels in Fig. 5) compares well with those obtained by solving the classical 2D SWEs over a high-resolution grid (upper panels in Fig. 5), both in terms of water depth and of timing of flooding. It has to be stressed that all the additional parameters needed by the dual porosity model (i.e., φ, Ψ L , Ψ T , and α) have been derived geometrically, and not by calibration.

Discussion and conclusions
First, it must be noted that the dual porosity model presented in this contribution is still in its embryonic stage. The treatment of the urban medium is obviously too simplistic to capture the complexity of real urban settlements, and the model suffers a number of severe limitations as well. For example, the use of two principal directions to model the conveyance porosity is likely an insufficient framework in order to characterize the actual directional dependence of resistance due to buildings and obstacles; the evaluation of effective values of conveyance porosity, averaged over the cell, is undoubtedly a hard task from both the geometric and hydraulic points of view. Furthermore, the Finite Element numerical scheme used in this study is neither suitable to deal with shock waves, nor with rapidly varying flows.
Notwithstanding these limitations, the preliminary results provided by the dual porosity model when applied to simplified applications are promising: with respect to high resolution simulations based on the classical 2D SWEs, the resistance exerted by urban layouts on the flow field is well reproduced by the porosity model without the need of calibration, and mesh-independence is preserved.