Numerical study of building drag dissipation formulations in the integral porosity shallow water model

. The integral porosity shallow water model is a type of porous shallow water model for urban ﬂood modeling, that deﬁnes two types of porosity, namely a volumetric porosity inside the computational cell and a conveyance porosity at each edge. Porosity terms are determined directly from the under-lying building geometry, hence buildings do not need to be discretized exactly. This enables simulations with signiﬁcantly reduced CPU time on meshes with cell sizes larger than the building size. Here, the macroscopic model view leads to an additional source term at the unresolved building-ﬂuid interface, yielding a building drag dissipation source term. In literature, several formulations for this term can be found. The integral porosity shallow water model is sensitive to the building drag dissipation, and using the drag parameters as a calibration parameter enhances the accuracy of model results. However, the ideal way to achieve this is still an open research question. In this contribution, we present a simple technique to estimate building drag dissipation that uses the conveyance porosity conﬁguration to estimate the projected area inside the cell, which is then used in a drag force equation. The advantage of this approach is that it is computationally inexpensive, no additional parameters need to be stored, and only a single parameter has to be calibrated. The proposed approach is compared with drag dissipation formulations from existing literature in a laboratory experiment that features a dam-break against an isolated obstacle. The aim of the comparison is to evaluate present existing building drag dissipation models with regard to accuracy and computational cost.


Introduction
The shallow water equations with porosity are most commonly used for flood simulations in urban areas, e.g. [1][2][3][4][5][6][7][8][9][10][11]. Another type of shallow water equations with porosity are applied to rainfall-runoff simulations in natural catchments, [12][13][14]. The essential idea of these * e-mail: ilhan.oezgen@wahyd.tu-berlin.de equations is to apply a homogenization technique to approach the multiscale nature of the flow processes. This allows to move from the microscale to the macroscale and reduces the computational burden. In porous shallow water equations porosity parameters are used to maintain a detailed description of obstacles (buildings in urban areas or microtopography in natural catchments) at a coarse scale, cf. [1,3] and [12,15] for urban areas and natural catchments, respectively. In the following, our discussion will be limited to a specific porous shallow water model, called the integral porosity shallow water model (IP).
The IP model is derived in [3]. Here, a binary phase function is used both inside of the control volume as well as along its boundary to obtain a volumetric porosity inside the control volume and a conveyance porosity at each edge of the control volume. The porosity terms depend on the geometric properties of the urban area and the computational mesh, as they express the space available for the flow inside a computational cell or at a cell edge. There are two points to note about the IP model: (1) the porosity terms can only be obtained via integration over a control volume, that means that the resulting governing equations are in integro-differential form; (2) the dependency of the porosity terms on the computational mesh causes the IP model results to be very mesh-dependent. As a consequence of the first point, numerical solutions for the IP model can only be obtained with methods that solve equations in integral form. Obviously, a good candidate is the finite volume method. Attempts to reduce the mesh dependency of IP model results are made in [11,16], where the conveyance porosity of an edge is related to the smallest length parallel to the edge that is available to the flow.
This paper studies formulations for the building drag dissipation source term in the IP model. This term appears during the homogenization procedure, because the building-fluid interface can not be resolved anymore but its influence on the flow still has to be accounted for [3]. Several strategies to discretize the source term have been proposed in literature, such as isotropic formulations in [1,2,4,15] and anisotropic drag formulations in [7,17]. Current research agrees that anisotropic drag formulations should be preferred over isotropic ones [9,15,17]. The discretization of the drag dissipation source term is still an open issue. The difficulty comes from the fact that this process in the urban area is not well understood. Another difficulty is the reproduction of a complex process at a coarse scale, using averaged variables and porosity parameters. In addition, the lack of laboratory experiments specifically designed for this purpose makes it difficult to validate building drag dissipation models. In this paper, we introduce an anisotropic drag source term discretization technique, that is in its philosophy close to the one in [4]. Since all existing drag dissipation models are essentially wrong [9], our aim is to provide a simple yet useful anisotropic calculation of the drag dissipation source that avoids overparameterization.
First, we briefly present our mathematical models. Then, we compare this proposed formulation to the equations in [4] and [9] in a test case against measurement data and a high resolution shallow water model. Finally, we draw conclusions.

Integral porosity shallow water model
We introduce the notation ∂ t to indicate a temporal derivative and ∂ x and ∂ y to indicate a spatial derivative in the Cartesian coordinate system in x-and y-direction, respectively.
The governing equations of the IP model can be written as a balance law 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 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 bed elevation, are the bed 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. 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 = � q x , q y � T . The equations are solved with a first order Godunov-type finite volume method suitable for both structured and unstructured meshes, using a Harten, Lax and van Leer approximate Riemann solver with the contact wave restored (HLLC) [18] for numerical flux calculation.

Building drag dissipation source term
The building drag dissipation source term results from the macroscopic model point of view. In [4], an isotropic estimation of the term is given as where c denotes a coefficient of proportionality andw p denotes the width of the obstruction projected in flow direction, which is calculated as with n b being the number of edges and w p,k being the projected width at the edge k, calculated as Based on the first three letters of the name of the first author of the corresponding publication, we will refer to this approach as the SCH approach hereinafter. A concept sketch is given in Fig. 1 (top). The second approach from literature is presented in [9]. The building drag dissipation is calculated as where θ is the angle between the unit discharge vector and the street axis, and c d,x and c d,y are two free calibration parameters that represent drag coefficients of proportionality in x-and y-direction, respectively. This formulation is a generalization of the drag dissipation source term formula in [17]. We refer to this approach as the GUI approach.  Finally, we propose to calculate the drag dissipation source term by using which is very similar to the SCH approach. But in contrast, we calculate two projected areas w p,x and w p,y in x-and y-direction, respectively, as where ∆k is the characteristic length of the cell andψ is an average conveyance porosity defined asψ The porosity terms ψ k and ψ m are the porosity terms at the two edges that intersect a line drawn through the cell center in j-direction ( j = x, y). We illustrate the concept in Fig. 1 (bottom). Only the edges intersected by the dashed line are taken into account. The edges that are not intersected are neglected. The number of free calibration parameters is reduced to one, i.e. c. This approach will be referred to as the WHD approach. For triangular cells, one conveyance porosity term will go into the drag dissipation source term calculation twice. If square cells are used, this is not the case. If quads are used, one porosity term might go into the calculation twice, depending on the distortion of the mesh.
Finally, both the SCH and the WHD approach can be seen as special cases of the GUI approach with θ = π/2 and a specific way of calculating c d,x and c d,y .

Dam-break flow against an isolated obstacle
This test case consists of a laboratory experiment that has been reported in [19]. In its initial state, a reservoir with a water elevation ofη r = 0.4 m is separated from a downstream channel by a gate. At t = 0, this gate is quickly opened. Inside the downstream channel, a single rectangular block is located that represents a building. The block is impermeable and gets never submerged during the experiment. Opening the gate generates a dam-break wave propagating in the downstream channel. Water depths are measured at four measurement gauges around the building. The geometry of the domain is sketched in Fig. 2 (top). The roughness of the channel is reported with a Manning's coefficient of n = 0.01 sm −1/3 . A thin layer of water of aboutη c = 0.02 m is reported in the channel [19]. The computational mesh is an unstructured quad mesh with an element size of 0.5 m around the obstacle, an element size of 0.075 m in the reservoir and 0.0375 m in the channel at places that are not in the close surrounding of the obstacle. This gives a mesh with 44796 elements. The mesh is shown in Fig. 2 (bottom). In our model, the downstream boundary is an open Neumann boundary with ∂ x h = ∂ y h = 0. All other boundaries are closed Dirichlet-type boundaries.
We compare the three model results obtained by different drag dissipation approaches in a sensitivity analysis-like manner. For each approach, we vary the free calibration parameters and study the deviation between the results. We compare IP model results with a highresolution classical shallow water model (SWE) that uses an unstructured quad mesh with element size 0.02 m in the whole domain. The mesh consists of 177522 elements (about 4 times more elements as in the IP model). Same initial and boundary conditions apply.
The results are shown in Fig. 3, 4 and 5 for the WHD, SCH and GUI approach, respectively. The SCH approach appears to be more sensitive to c than the WHD approach. The same values of c yield a higher water depth at gauge 1 for the SCH approach. Gauge 1 is the most sensitive gauge, other gauges are not influenced as significantly by the drag coefficient of proportionality.
For the GUI approach, we report that the value of c D,y does not influence the results significantly. Varying it from c D,y = 0.1 to c D,y = 2 changes the model results only insignificantly. We omit plotting these curves, because they are visually indistinguishable. Instead, Fig. 5 shows results for varying c D,x -values with a fixed c D,y = 0.5.
At gauges 3, 4 and 5, results for all approaches are similar. Increasing the drag coefficient of proportionality slightly decreases the water depth at gauge 3 and 5. At gauge 4, it decreases the water depth at the beginning of the simulation but increases it at later stages.
The difference of the influence of the drag coefficient of proportionality comes from the location of the gauges. The gauge located in front of the obstacle, namely gauge 1, is influenced the most by the drag dissipation source term. Gauge 5 is located behind the obstacle, and therefore is not sensitive to the drag dissipation source term. We think that at this location, the influence of the porosity terms outweighs the influence of the drag dissipation  mechanism. Gauge 3 and gauge 4 are located at the sides of the obstacle. The water reaching gauge 3 is marginally less blocked by the obstacle than the water at gauge 4, because it is located slightly more in the upstream direction. We observe that gauge 4 shows slightly more sensitivity to the drag dissipation.
Overall, fair agreement between IP model, SWE model and experimental data was obtained. In terms of computational efficiency, the WHD approach requires slightly more floating point operations (FLOP) than the SCH approach to determine the conveyance porosity terms in each direction via intersection. This is the cost for introducing anisotropy and can be approached in two ways: (1) always re-determine the intersection points from scratch for each time step; or (2) store the porosity terms and look them up afterwards. The first approach requires less memory but more FLOP. The GUI approach requires determining the angle between vector of the unit discharge and street axis. Consequently, a number of FLOP has to be carried out at each time step to calculate this angle. If the buildings are not aligned but distributed in the domain in a more random way, which is expected in real world cases, determining the street axis may require postprocessing with a GIS.

Conclusions
We presented a novel building drag dissipation approach for the integral porosity shallow water model. We compared this approach to two recent approaches from literature in a case study considering dam-break against an isolated obstacle. Results indicate that the presented approach is at least as useful as existing drag dissipation models.
We report that for the investigated test case, the presented approach shows less sensitivity to the free calibration parameter than the approach in [4] (SCH). Comparison with the gener-alized drag formula in [9] (GUI) is not as straight-forward, as the free calibration parameters represent different processes.
It is interesting to note that the WHD approach is close to the GUI approach in terms of accuracy, while it uses only a single calibration parameter. This implies that similar results can be obtained with significantly less effort.
The aim of deriving the presented drag dissipation formula is not a deeper process understanding. Further, we do not claim that this approach is any more correct than existing drag dissipation models. The benefit of the presented approach is that it estimates the drag dissipation term based on values that are already calculated, i.e. the conveyance porosity terms. This gives a useful anisotropic estimation of the drag dissipation source term, that can easily be implemented into existing models. All drag dissipation models yield similar model results for water depth for the investigated case.
A comparison of the approaches by means of automatic calibration of the free calibration parameters in each approach is currently in preparation.