Vincent Guinot

Brett F Sanders

Jochen E Schubert

Dual integral porosity shallow water model for urban ood modelling

With CPU times 2 to 3 orders of magnitude smaller than classical shallow water-based models, the shallow water equations with porosity are a promising tool for large-scale modelling of urban oods. In this paper, a new model formulation called the Dual Integral Porosity (DIP) model is presented and examined analytically and computationally with a series of benchmark tests. The DIP model is established from an integral mass and momentum balance whereby both porosity and ow variables are dened separately for control volumes and boundaries, and a closure scheme is introduced to link control volume-and boundary-based ow variables. Previously developed Integral Porosity (IP) models were limited to a single set of ow variables. A new transient momentum dissipation model is also introduced to account for the eects of sub-grid scale wave action on porosity model solutions, eects which are validated by ne-grid solutions of the classical shallow-water equations and shown to be important for achieving similarity in dam-break solutions. One-dimensional numerical test cases show that the proposed DIP model outperforms the IP model, with signicantly improved wave propagation speeds, water depths and discharge calculations. A two-dimensional eld scale test case shows that the DIP model performs better than the IP model in mapping the oods extent and is slightly better in reproducing the anisotropy of the ow eld when momentum dissipation parameters are calibrated.

Introduction

Since their inception in the nineteen-nineties [START_REF] Dena | A new set of equations for very shallow water and partially dry areas suitable to 2D numerical models[END_REF], shallow water models with porosity have arisen as a promising formalism for the inclusion of subgrid-scale topographical and geometric features [START_REF]Two-dimensional shallow ow equations for partially dry areas[END_REF][START_REF] Bates | Development and testing of a sub-grid scale model for moving boundary hydrodynamic problems in shallow water[END_REF] in two-dimensional ow simulations. Their usefulness to tackle urban ood modelling has been illustrated by a number of applications, with reported CPU acceleration factors of two to three orders of magnitude compared to usual two-dimensional models solving the shallow water equations [START_REF] Guinot | Flux and source term discretisation in twodimensional shallow water models with porosity on unstructured grids[END_REF][START_REF] Guinot | Multiple porosity shallow water models for macroscopic modelling of urban oods[END_REF][START_REF] Kim | Urban ood modeling with porous shallow-water equations: A case study of model errors in the presence of anisotropic porosity[END_REF][START_REF] Lhomme | One-dimensional, two-dimensional and macroscopic approaches to urban ood modelling[END_REF][START_REF] Sanders | Integral formulation of shallow water equations with anisotropic porosity for urban ood modelling[END_REF][START_REF] Soares-Frazao | Twodimensional shallow water models with porosity for urban ood modelling[END_REF][START_REF] Özgen | Shallow water equations with depth-dependent anisotropic porosity for subgrid-scale topography[END_REF]. The shallow water equations with porosity have also motivated a number of numerical developments [START_REF] Benkhaldoun | A non-homogeneous Riemann solver for shallow water equations in porous media[END_REF][START_REF] Cea | Unstructured nite volume discretization of two-dimensional depthaveraged shallow water equations with porosity[END_REF][START_REF] Finaud-Guyot | An approximate-state Riemann solver for the two-dimensional shallow water equations with porosity[END_REF][START_REF] Guinot | Flux and source term discretisation in twodimensional shallow water models with porosity on unstructured grids[END_REF][START_REF] Lhomme | One-dimensional, two-dimensional and macroscopic approaches to urban ood modelling[END_REF].

The rst versions of shallow water models with porosity for simulating urban oods incorporated a single porosity [START_REF] Guinot | Flux and source term discretisation in twodimensional shallow water models with porosity on unstructured grids[END_REF][START_REF] Hervouët | Modelling urban areas in dam-break oodwave numerical simulations[END_REF]. Although the discretization of geometry-induced source terms in the Single Porosity (SP) model is more complex than for the usual shallow water equations [START_REF] Cea | Unstructured nite volume discretization of two-dimensional depthaveraged shallow water equations with porosity[END_REF][START_REF] Finaud-Guyot | An approximate-state Riemann solver for the two-dimensional shallow water equations with porosity[END_REF][START_REF] Guinot | Flux and source term discretisation in twodimensional shallow water models with porosity on unstructured grids[END_REF][START_REF] Lhomme | One-dimensional, two-dimensional and macroscopic approaches to urban ood modelling[END_REF], the wave propagation speeds are identical. At the time the rst nite volume discretization of the SP model was published [START_REF] Guinot | Flux and source term discretisation in twodimensional shallow water models with porosity on unstructured grids[END_REF], it did not fully account for the eect of ow path obstruction by the buildings in addition to the storage eect. It was thus proposed that two dierent porosities should be used: a domain-based porosity to account for mass and momentum storage and a boundary-based, connectivity porosity to account for mass and momentum transport [START_REF] Lhomme | One-dimensional, two-dimensional and macroscopic approaches to urban ood modelling[END_REF]. An eigenvalue analysis of the resulting system showed that wave propagation speeds [START_REF] Lhomme | One-dimensional, two-dimensional and macroscopic approaches to urban ood modelling[END_REF] (and the behaviour of the ow solutions) should be signicantly dierent from those of the SP model. However, no practical guidelines could be provided to infer the connectivity porosity from the urban geometry. Independently from this research, the Integral Porosity (IP) model was proposed almost at the same time and published shortly afterwards [START_REF] Sanders | Integral formulation of shallow water equations with anisotropic porosity for urban ood modelling[END_REF]. The two key issues brought to the fore by [START_REF] Sanders | Integral formulation of shallow water equations with anisotropic porosity for urban ood modelling[END_REF] are that (i) the SP model is based on the implicit assumption that the Reference Elementary Volume (REV) [START_REF] Bear | Dynamics of uids in porous media[END_REF] exists in the urban environment, (ii) the SP formalism is essentially isotropic. The IP model allows these two drawbacks to be eliminated because the integral formulation makes the connectivity porosity a local, directional, and deterministic descriptor based on the intersection of the computational mesh with ow obstructions, and the storage porosity is also given deterministically. This limited the need for calibration to directional building drag parameters to achieve closure in eld-scale applications, whereas SP modelling studies had also examined the calibration of porosity [START_REF] Soares-Frazao | Twodimensional shallow water models with porosity for urban ood modelling[END_REF]. Recently, an IP model that accounts for depth-dependent porosity has been developed as a framework for both subgrid topographic variability and blockage features such as buildings [START_REF] Özgen | Shallow water equations with depth-dependent anisotropic porosity for subgrid-scale topography[END_REF].

A one-dimensional analysis of the wave propagation speeds of the dierential form of the IP model equations [START_REF] Guinot | Macroscopic modelling of urban oods[END_REF] conrms that the IP model structure is similar to the dierential formulation proposed in [START_REF] Lhomme | One-dimensional, two-dimensional and macroscopic approaches to urban ood modelling[END_REF]. Nevertheless, numerical experiments reported in a recent publication [START_REF] Özgen | Urban ood modeling using shallow water equations with depth-dependent anisotropic porosity[END_REF] and in the present paper show that the wave propagation speeds and discharge predictions in the IP model are signicantly dierent from those in rened wave propagation simulations. The wave propagation speeds are the eigenvalues of the Jacobian matrix of the ux vector with respect to the conserved variable vector.

Consequently, erroneous wave speeds usually reect erroneous mass and/or momentum ux models.

In this paper, the introduction of a dual set of ow variables corresponding to integrals over cells and boundaries, and a closure scheme which links the two solutions through the dual porosity denition, is shown to resolve this problem. The resulting Dual Integral Porosity (DIP) model supports a signicant improvement in accuracy in verication test cases, and shows advantages compared to the IP model in a eld-scale test problem.

The remainder of the paper is organised as follows. Section 2 presents the DIP formalism along with the domain-boundary closure model and a modication of the momentum uxes to account for transient momentum dissipation which cannot be modelled like friction with a sink term. Section 3 presents the verication of the DIP model in idealized congurations where analytical solutions are known, as well as comparisons to SP and IP models. Section 4 presents a two-dimensional eld scale test where the practical utility and versatility of the model is validated, i.e., ability to control anisotropic solutions through calibrated parameters. Sections 5 and 6 are devoted to a discussion and conclusions, respectively.

2 Integral porosity shallow water models

Original formulation

The conservation form of the model is obtained from a mass and momentum balance applied to a shallow water domain Ω with boundary Γ. The domain is partitioned into two regions: a solid region occupied by buildings, and a water region. Mass and momentum conservation leads to [START_REF] Sanders | Integral formulation of shallow water equations with anisotropic porosity for urban ood modelling[END_REF]:

∂ t ˆΩ εhdΩ + ˆΓ εq.ndΓ = 0 (1a) ∂ t ˆΩ εqdΩ + ˆΓ ε (u.n) q + g 2 h 2 n dΓ = ˆΩ εs Ω dΩ + ˆΓ εs Γ dΓ (1b) q = hu = hu hv , s Γ = g 2 h 2 η0 n (1c) s Ω = -g (h -h η0 ) ∇z b -c b D + c f D |u| u (1d)
where n is the normal unit vector, q and u are respectively the unit discharge and ow velocity vectors, s Ω and s Γ are respectively the domain and boundary source terms, h is the water depth, u and v are respectively the xand y-velocities, z b is the bottom level, ε is the porosity indicator (equal to 0 in the solid areas and 1 in the water region), η 0 is the average free surface elevation over the domain Ω, h η0 = η 0 -z b , c b D and c f D are respectively the drag coecients for the building and bottom friction-induced stresses. The porosity over the sets Ω, Γ is dened as

φ S ≡ ´Sw dS ´S dS = ˆS ε (x, y) dS, S = Ω, Γ (2)
As stressed in [START_REF] Sanders | Integral formulation of shallow water equations with anisotropic porosity for urban ood modelling[END_REF], the water storage porosity φ Ω and the water conveyance porosity φ Γ are dierent in the general case. By denition, φ Γ scales transport capacity and should be taken smaller than φ Ω . The dierence φ Ω -φ Γ thus represents an estimate for the fraction of the domain that does not contribute to transport. The constraint φ Γ ≤ φ Ω is further substantiated by a wave propagation speed analysis, see Section 2.5.

Although questionable when a Reference Elementary Volume does not exist [START_REF] Sanders | Integral formulation of shallow water equations with anisotropic porosity for urban ood modelling[END_REF], the dierential form of the equations for this model allows the behaviour of the solutions to be analysed. In the limit of an innitesimal control volume Ω, using the divergence theorem leads to a system in the form [START_REF] Lhomme | One-dimensional, two-dimensional and macroscopic approaches to urban ood modelling[END_REF] ∂

t (φ Ω h) + ∇ • (φ Γ q) = 0 (3a) ∂ t (φ Ω q) + ∇ • (φ Γ F) = φ Ω s Ω + ∇ • (φ Γ s Γ ) (3b)
F = hu 2 + g 2 h 2 huv huv hv 2 + g 2 h 2 = u ⊗ u + g 2 h 2 I ( 3c
)
where I is the identity matrix. The wave propagation speeds in the one-dimensional conguration are the eigenvalues of the system [START_REF] Guinot | Macroscopic modelling of urban oods[END_REF][START_REF] Lhomme | One-dimensional, two-dimensional and macroscopic approaches to urban ood modelling[END_REF]:

λ 1 = φ Γ φ Ω (u -c) , λ 2 = φ Γ φ Ω u, λ 3 = φ Γ φ Ω (u + c) , c = (gh) 1 2 (4)
This system is hyperbolic. Its eigenvalues are those of the classical shallow water equations, multiplied by the ratio φΓ φΩ . The particular case φ Γ = φ Ω yields the wave propagation speeds of the classical Saint Venant equations. This property is used in nite volume discretizations of the single porosity equations to derive approximate Riemann solvers [START_REF] Cea | Unstructured nite volume discretization of two-dimensional depthaveraged shallow water equations with porosity[END_REF][START_REF] Finaud-Guyot | An approximate-state Riemann solver for the two-dimensional shallow water equations with porosity[END_REF][START_REF] Guinot | Flux and source term discretisation in twodimensional shallow water models with porosity on unstructured grids[END_REF]. If the mesh is designed such that φ Γ > φ Ω , the wave propagation speeds are larger than those of the classical shallow water equations. Such a behaviour is non-physical.

Dual integral porosity model formulation

The DIP model overcomes ambiguity in the denition of the uid velocity, u, at domain boundaries (versus within domain boundaries) with the introduction of dual ow variables (u Ω , u Γ ) and substitution into equations (1a-b) as follows,

∂ t ˆΩ εh Ω dΩ + ˆΓ εq Γ .ndΓ = 0 (5a) ∂ t ˆΩ εq Ω dΩ + ˆΓ ε (u Γ .n) q Γ + g 2 h 2 Γ n dΓ = ˆΩ εs Ω dΩ + ˆΓ εs Γ dΓ (5b)
For the number of equations to match the number of unknowns, a closure model is needed between the domain-averaged variable vector u Ω and the boundary-averaged variable vector u Γ .

Closure model for the unit discharge. Consider an elemental boundary segment dΓ = [AB]

(Figure 1). Typically, [AB] is the edge of a computational cell in a nite volume or nite element discretization. The total discharge across [AB] is

Q AB = φ Γ q Γ .n dΓ (6)
Consider now the segment [A B ], obtained by shifting [AB] into the domain Ω. For an innitesimal shift, the discharge across [A B ] is equal to the discharge across [AB]. It must be computed using the areal porosity φ Ω because the segment is located within the domain

Q A B = φ Ω q Ω .n dΓ (7)
Figure 1: Closure model for the unit discharge. Denition sketch.

Continuity requires that Q AB = Q A B . This leads to the sucient condition

q Γ = φ Ω φ Γ q Ω ⇔ u Γ = φ Ω φ Γ h h Γ u Ω (8)
Note that modifying only the normal component of the ow velocity would suce to satisfy continuity.

However, doing so would make the direction of q Γ articially dependent on the direction of the interface.

Therefore, both the normal and tangent unit discharges are modied by the same factor. With this correction, the mass ux across the boundary veries

ˆΓ φ Γ q Γ .n dΓ = ˆΓ φ Ω q Ω .n dΓ (9)
Closure model for the water depth. Substituting the closure model ( 8) into the expression for the momentum ux leads to

ˆΓ φ Γ (u Γ .n) q Γ dΓ = φ 2 Ω φ Γ ˆΓ h Ω h Γ (u Ω .n) q Ω dΓ (10)
A model must be proposed for the ratio hΩ hΓ in the momentum equation. The simplest possible model, h Γ = h Ω , is retained. It is based on the assumption that building-induced ow obstruction has negligible eects on the free surface elevation and that the reduction in the ow cross-sectional area is fully balanced with an increase in the ow velocity. This is a reasonable approximation for small values of the Froude number, based on gradually varied ow theory, and experiments have shown that supercritical (high Froude) ow impinging on clusters of buildings causes upstream shock waves that surround buildings by locally subcritical (low Froude) ow [START_REF] Testa | Flash ood ow experiment ni a simplied urban district[END_REF]. The expression for the momentum ux becomes

ˆΓ φ Γ (u Γ .n) q Γ dΓ = φ 2 Ω φ Γ ˆΓ (u Ω .n) q Ω dΓ (11)

Accounting for transient momentum dissipation

Solving the frictionless shallow water equations over street networks with a constant street width in the presence of piecewise constant initial and/or boundary conditions is shown in [START_REF] Guinot | Multiple porosity shallow water models for macroscopic modelling of urban oods[END_REF] to yield solutions with an almost perfectly self-similar (x, t) behaviour. The waves in the solution are seen to propagate at a constant speed, making the solution at the scale of the street spacinperiod dependen on the ratio

x/t. The self-similarity of the solutions is not exact in that a street network with xed width and spacing is not invariant when the x-coordinate is transformed, but self-similarity is observed for the solutions averaged on the scale of a building block. Two such examples are presented in [START_REF] Guinot | Multiple porosity shallow water models for macroscopic modelling of urban oods[END_REF]. The rst deals with the propagation of a ood bore into an idealized urban network. The second consists of an urban dambreak problem. As illustrated by Section 3, the numerical solutions obtained by solving the frictionless shallow water equations on a ne grid cannot be reproduced by a porosity model that does not incorporate a momentum dissipation mechanism of some kind. Using the proposed DIP model without friction will lead to wrongly estimated wave propagation speeds, water depths and ow velocities, especially in the case of positive waves (rising water levels).

An explanation for this behaviour was introduced in [START_REF] Guinot | Multiple porosity shallow water models for macroscopic modelling of urban oods[END_REF]: when a positive wave propagates into a built area (Figure 2), multiple wave reections occur against the building walls and generate moving bores (Figure 2, top). The wave reections generate reactions to the hydrostatic pressure forces from the building walls. In Figure 2, the reaction from wall W 2 is larger than that from wall W 1 because the reected bore is deeper (Figure 2, bottom). In the same way, the reaction from wall W 4 is larger than that from wall W 3 . The resulting force is opposed to the average ow velocity and thus contributes to dissipate momentum. Steady state ow and decreasing water levels do not generate bores and no signicant dissipation occurs due to wave reection. This mechanism was conrmed in [START_REF] Guinot | Multiple porosity shallow water models for macroscopic modelling of urban oods[END_REF] by a detailed inspection of the water depth, velocity and hydraulic head elds in the rened shallow water simulations. Although the dissipation is not instantaneous in reality, it is assumed to occur on a time scale smaller than or of the same order of magnitude as the time scale needed by the wave to travel between two successive building rows. The assumption of instantaneous dissipation is thus a convenient approximation. The tests presented in Section 3 indicate that the proposed model is reasonably accurate even when the ow is supercritical. Classical head loss models that are direct functions of the ow velocity are well-known to not allow for self-similar solutions in the case of piecewise constant initial/boundary value problems, and thus cannot reproduce the losses from subgrid-scale bores described above. The only possibility to retain the self-similar character of the solutions is for the momentum dissipation model to be divergencefree. In other words, the governing equations must consist only in a modication of the momentum ux and/or the conserved variable (e.g. by introducing a so-called rigidity, or inertia matrix). The proposed momentum dissipation model satises this condition. It is derived from that proposed in the multiple porosity model [START_REF] Guinot | Multiple porosity shallow water models for macroscopic modelling of urban oods[END_REF] under the so-called Local Equilibrium Assumption (LEA). Consider again a positive wave entering the averaging volume shown in Figure 2. Within an innitesimal time interval dt, the total momentum in the averaging volume is subjected to a positive variation φ Ω dq Ω , and if the fraction of transported momentum that is dissipated in dead zones is dened by µ ∈ [0, 1], the net variation in the average momentum over the averaging volume is thus (1 -µ) φ Ω dq Ω . This is equivalent to multiplying the divergence of the momentum uxes by a factor (1 -µ) ≤ 1. When the averaging volume is subjected to a zero or negative variation (momentum and volume are withdrawn), no dissipation occurs. Consequently, the above correction is applied only for a rising water depth, that is, when ∂ t h > 0, or equivalently when the divergence of the ow eld is negative ∇ • q < 0.

A B W 1 W 2 W 3 W 4 A B
The anisotropic character of most urban layouts calls for generalizing the proposed momentum dissipation term in the form of a second-order tensor:

D = µ xx µ xy µ yx µ yy = B µ 1 0 0 µ 2 B -1 dq Ω , B = cos θ 1 cos θ 2 sin θ 1 sin θ 2 (12)
where θ 1 and θ 2 are the angles of the principal directions of the urban layout with the x-axis. Note that this does not necessarily result in a symmetric tensor. The innitesimal momentum variation becomes

(I -D) φ Ω dq Ω (13)

Anisotropic building drag model

The proposed model also incorporates an improvement to the building drag force model. Although deemed of secondary importance compared to the closure and the momentum dissipation models presented in Subsections 2.2-3, the revised building drag model allows the accuracy of integral porosity models to be increased further. This is illustrated by the computational examples in Section 4. The building drag model originally implemented in the IP formulation [START_REF] Sanders | Integral formulation of shallow water equations with anisotropic porosity for urban ood modelling[END_REF] uses the following momentum source term (equation (1d)

s b D = -C b D |u| u (14)
Firstly, the source term be made depth-dependent. For a given ow velocity, building drag force can be expected to be proportional to the contact surface between the water and the buildings. The contact surface being proportional to the water depth, it is proposed that the source term be scaled by the water depth. Secondly, many urban street networks are observed to exhibit anisotropy, thus motivating a tensor-based description of the building drag terms. The following formula is thus proposed

s b D = -C b D h h ref |u| u, C b D = C b D,xx C b D,xy C b D,yx C b D,yy (15)
where C b D is a drag coecient tensor and h ref is a reference water depth used for scaling purposes. In what follows, h ref is set to 1 m.

2.5

Final formulation and numerical procedure

Incorporating the new closure model, the momentum dissipation terms and the revised building drag model into the mass and momentum balances lead to the following set of governing equations

∂ t v + ˆΓ MFn dΓ = ˆΩ εs Ω dΩ + ˆΓ εs Γ dΓ (16a) v ≡ φ Ω   h hu hv   , F ≡ φ Γ   βhu βhv β 2 hu 2 + g 2 h 2 β 2 huv β 2 huv β 2 hv 2 + g 2 h 2   , M ≡   1 0 0 0 1 -µxx -µxy 0 -µyx 1 -µyy   (16b) β ≡ φ Ω φ Γ (16c) s Ω =   0 g (h -hη 0 ) S 0,x -c f D |u| u g (h -hη 0 ) S 0,y -c f D |u| v   -   0 0 C b D,xx C b D,xy C b D,yx C b D,yy   h h ref |u| u (16d) s Γ = g 2 h 2 η 0   0 nx ny   (16e)
where n x and n y are respectively the xand y-components of the normal unit vector to the boundary, S 0,x = -∂ x z b and S 0,y = -∂ y z b are respectively the xand y-bottom slopes. Note that the bottom friction term originally proposed in [START_REF] Sanders | Integral formulation of shallow water equations with anisotropic porosity for urban ood modelling[END_REF] is used in the DIP model too for the sake of model comparison.

The equations are discretized using an explicit nite volume approach [START_REF] Cea | Unstructured nite volume discretization of two-dimensional depthaveraged shallow water equations with porosity[END_REF][START_REF] Finaud-Guyot | An approximate-state Riemann solver for the two-dimensional shallow water equations with porosity[END_REF][START_REF] Guinot | Multiple porosity shallow water models for macroscopic modelling of urban oods[END_REF][START_REF] Guinot | Flux and source term discretisation in twodimensional shallow water models with porosity on unstructured grids[END_REF][START_REF] Sanders | Integral formulation of shallow water equations with anisotropic porosity for urban ood modelling[END_REF]. The frictionless part of the equations and the friction-induced source terms are solved sequentially using a time splitting procedure.

Frictionless terms. The frictionless part of the equations is rst solved as

v n+1/2 i = v n i + ∆t A i j∈N (i) L ij M i F ij .n ij + ∆t j∈N (i) S Γ,ij n ij (17a) S = L ij   0 0 s Γ,ij 0 0 s Γ,ij   (17b) s Γ,ij ≡ -g φ Ω,i -φ Γ,ij 2 h 2 i + φ Γ,ij max 0, η i -z b,j + η i -z b,i 2 max 0, min η i , z b,j -z b,i (17c)
where the subscript i is the cell number, the superscript n denotes the time level, N (i) is the set of neighbouring cells for the cell i, A i is the plan view area of the cell i, L ij and n ij are respectively the length and normal unit vector of the interface between the cells i and j, ∆t is the computational time step. The ow variables in equations (17a-17b) are cell-averaged quantities and should bear the subscript Ω. It is omitted in the notation for the sake of clarity. The discretized source term S Γ,ij accounts for the variations in the porosity and bottom levels across the interface. It thus incorporates the discretization of the source term s Γ and the contribution of the bottom slope in s Ω , as detailed in [START_REF] Finaud-Guyot | An approximate-state Riemann solver for the two-dimensional shallow water equations with porosity[END_REF].

Friction terms. The solution is updated to account for the friction term using time splitting with an exact solution [START_REF] Guinot | Multiple porosity shallow water models for macroscopic modelling of urban oods[END_REF]. The ow variables are rst rotated into the coordinates formed by the principal directions of the building drag tensor. In this coordinate system, the equations become

∂ t φ Ω   h hu 1 hu 2   = -φ Ω     0 c f D + C b D,1 h h ref |u| u 1 c f D + C b D,2 h h ref |u| u 2     (18)
where the subscripts 1 and 2 denote the component of the velocity along the wo principal directions.

Simplifying by φ Ω leads to

∂ t hu p = -c f D + C b D,p h h ref |u| u p , p = 1, 2 (19)
Under the assumption of small time steps, the equation is linearised into

∂ t (hu p ) i = -c f D + C b D,p h h ref |u n i | u p , p = 1, 2 (20)
with analytical solution

h n+1 i = h n+1/2 i (21a) (hu p ) n+1 i = exp -c f D + C b D,p h h ref |u n i | ∆t (hu p ) n+1/2 i (21b)
Remark on ux computation. The reader's attention is drawn to the computation of the uxes in (17a). The nite volume implementations presented in the literature [START_REF] Finaud-Guyot | An approximate-state Riemann solver for the two-dimensional shallow water equations with porosity[END_REF][START_REF] Cea | Unstructured nite volume discretization of two-dimensional depthaveraged shallow water equations with porosity[END_REF][START_REF] Guinot | Multiple porosity shallow water models for macroscopic modelling of urban oods[END_REF][START_REF] Guinot | Flux and source term discretisation in twodimensional shallow water models with porosity on unstructured grids[END_REF][START_REF] Sanders | Integral formulation of shallow water equations with anisotropic porosity for urban ood modelling[END_REF] use Riemann solvers to compute the ux tensor F. Noticing that the ux tensor in (??) can be rewritten as

F ≡ φ Γ   hU hV hU 2 + g 2 h 2 hU V hU V hV 2 + g 2 h 2   , U V ≡ β u v (22)
the structure of the ux tensor is identical to that of the classical two-dimensional shallow water equations. It stems as an obvious consequence that the uxes can be computed using any standard Riemann solver for the classical 2D shallow water equation, by dening the left and right states of the Riemann problem as

v L ≡ φ Ω,L   h L β L h L u L β L h L v L   , v R ≡ φ Ω,L   h R β R h R u R β R h R v R   (23)
where φ Ω,S (S = L, R) denotes φ Ω on the side S of the interface.

The discretization (17a) being explicit, it is subjected to a stability constraint based on the CFL of the fastest waves. The CFL condition for an unstructured grid is given in [START_REF] Guinot | An eigenvector-based linear reconstruction scheme for the shallow water equations on two-dimensional unstructured meshes[END_REF] and will not be recalled here. The formulas for the wave propagation speeds are derived in the next section.

Eigenvalue analysis

For the sake of clarity, the one-directional analysis of the equations is presented. The ux vector f is thus the rst column of the ux tensor F. Retaining only the conservation part of the equations, applying the integral form (5a, 5b, 11) to an innitesimal control volume, applying the divergence theorem under the assumption of continuous and dierentiable solutions leads to a system in the form

∂ t v + ∂ x (Mf ) = 0 (24a) v ≡ φ Ω   h hu hv   , f ≡ φ Ω   hu (1 -µ) βhu 2 + (1 -µ) g 2β h 2 (1 -µ) βhuv   = φ Ω   hu ηhu 2 + η β 2 g 2 h 2 ηhuv   (24b)
where η ≡ (1 -µ) β. It is interesting to notice that only the storage porosity φ Ω appears in the continuity equation (in contrast with the original IP model). This makes the new continuity equation consistent with that derived by Dena [START_REF]Two-dimensional shallow ow equations for partially dry areas[END_REF]. Table 1 summarizes the various available formulations for u and f , the corresponding Jacobian matrices A and their eigenvalues λ in the one-dimensional case.

Model v f A λ Shallow water   h hu hv     h hu 2 + g 2 h 2 huv     0 1 0 c 2 -u 2 2u 0 -uv v u   λ1 = u -c λ2 = u λ3 = u + c Single porosity φ   h hu hv   φ   h hu 2 + g 2 h 2 huv     0 1 0 c 2 -u 2 2u 0 -uv v u   λ1 = u -c λ2 = u λ3 = u + c Integral porosity φΩ   h hu hv   φΓ   h hu 2 + g 2 h 2 huv   1 β   0 1 0 c 2 -u 2 2u 0 -uv v u   λ1 = 1 β (u -c) λ2 = 1 β u λ3 = 1 β (u + c) Dual integral porosity φΩ   h hu hv   φΩ   hu ηhu 2 + η β 2 g 2 h 2 ηhuv     0 1 0 η β 2 c 2 -ηu 2 2ηu 0 -ηuv ηv ηu   λ1 = ηu -c λ2 = ηu λ3 = ηu + c
Table 1: Summary of the various formulations. For the sake of conciseness, only the one-dimensional version is presented.

Note: c ≡ (gh) 1/2 , c ≡ η 1 2 (η -1) u 2 + c 2 β 2 1 2 .
The present dual formulation implicitly rules out the conguration φ Γ > φ Ω . Assuming φ Γ > φ Ω induces β < 1, which leaves room for situations where η < 1. Then combinations (c, u) exist such that (η -1) u 2 + c 2 β 2 < 0 and the λ 1 and λ 3 are complex, inducing a loss of hyperbolicity and making the initial value problem ill-posed. This is a strong point in favour of the modied model in that meaningless porosity congurations are ruled out by assumption. Note that this is not the case with the IP model, that allows for β < 1, although such a situation is non-physical (see Section 2.1).

3 Verication Test Cases

Positive Wave

The purpose of this test case is to verify the DIP model formulation including the transient dissipation model proposed in the previous section. It in inspired from [START_REF] Guinot | Multiple porosity shallow water models for macroscopic modelling of urban oods[END_REF]. A periodic, orthogonal street network is generated (Figure 3, left). The spatial period in the longitudinal and transverse direction are respectively L and 2L, with L = 100m. The street width is W = 50m in both directions. The bottom is horizontal, friction is assumed negligible. The initial and boundary conditions are

v (x, 0) = [h 0 , 0, 0] T ∀x > 0 (25a) [hu, hv] T (0, t) = [q b , 0] ∀t > 0 (25b)
where h 0 and q b are respectively the initial water depth and boundary unit discharge. A positive wave is generated by setting q b = 1m 2 s -1 with h 0 = 1m. The two-dimensional frictionless shallow water equations are solved numerically over a 2.5m×2.5m computational grid (Figure 3, right). Symmetry considerations allow a complete solution to be obtained by meshing only half a period in the transverse direction, 0 ≤ y ≤ L. When solving the porosity-based shallow water equations, a much coarser computational grid is used, with cell widths (∆x, ∆y) = (20m, 100m). in Simulation (v) was calibrated to best reproduce both the average bore wave propagation speed and the water level at the model boundary. Note that, in the reference Simulation (i), the water level oscillates about an average value that is almost constant with x. The midpoint of the bore is located at x = 710m for t = 300s and x = 1400m for t = 600s. This conrms the self-similar character of the reference solution at the scale of the network period (L = 100m). As expected from the wave speed analysis, the IP model yields a signicantly underestimated wave propagation speed, even when the head loss is assumed zero (C b D = 0). As a result, the positive wave in Simulation (ii) is slower and the free surface elevation is higher than the reference case, Simulation (i). Increasing C b D would only slow down the wave propagation speed even more and is therefore not an option. The DIP model without transient dissipation (Simulation (iii)), in contrast, yields a slightly overestimated wave speed and a subsequently slightly underestimated free surface elevation.

In Simulation (iv), tuning the building drag coecient C b D improves prediction of the bore speed but destroys the self-similar character of the ow solution, with a sloping free surface and a time-dependent free surface elevation at the boundary. Only Simulation (v) predicts the bore speed and boundary water level accurately while preserving the self-similar character of the solution.

Negative Wave

The purpose of this test is to validate the absence of transient momentum dissipation when ∂ t h < 0, as formulated in the DIP model. The physical layout and initial/boundary conditions are identical to those in Section 3.2, with the exception that q b = -0.5m 2 s -1 in order to create a negative wave. The following simulations are carried out (Table 3): (i) the reference simulation whereby the classical 2D shallow water equations are solved on the 2.5m×2.5m rened grid, (ii) IP simulation without building As in the positive wave test case, the IP model (ii) yields signicantly underestimated wave propagation speeds. Consequently, the boundary value for the free surface elevation is underestimated. Simulation (iii) shows that when dissipation is added to the DIP model for ∂ t h ≤ 0, wave speeds and water levels are also underestimated in the region of constant state behind the wave. In Simulation (iv), the presence of a non-zero source term on the right-hand side of the momentum equation destroys the self-similar character of the solution. The water depth at the boundary decreases with time as a result of the increasing length over which friction occurs. The proposed DIP model formulation (Simulation (v)) is the only one where the self-similar character of the solution is preserved and the water levels and average location of the negative wave are accurately predicted. Comparing Simulations (iii) and (v) conrms that the proposed momentum dissipation model should be triggered only when the water level decreases.

Although Simulation (v) allows the average propagation speed of the negative wave to be estimated accurately, it fails to represent correctly the spreading of the wave. This failure is explained by inspecting the velocity eld in the reference simulation (i). Figure 6 shows the velocity distribution for 725m ≤ x ≤ 925m at t = 600s. The gure clearly indicates that the ow velocity is not uniform across the street. Along the building walls, a low velocity region can be observed, with a width approximately 20% of the width of the street. In this low velocity region, the velocity is approximately twice as small as the velocity along the centreline of the street. This low velocity region is due to the water being abstracted from the lateral streets into the main street. The nearly immobile water coming from the lateral street mixes with the fast owing water in the main street, causing a decrease in the ow speed near the interface between the two street networks. This non-uniform velocity distribution is not accounted for in the proposed closure model [START_REF] Cea | Unstructured nite volume discretization of two-dimensional depthaveraged shallow water equations with porosity[END_REF], that assumes a uniform ow velocity over the cross-sections between the buildings. Improving the closure model ( 8) would require that a correction coecient be applied to decrease φ Γ , which would result in increasing u Γ . Increasing u Γ would Figure 6: Propagation of a negative wave into an idealized urban network. Reference simulation (i).

Velocity at t = 600s next to the left-hand boundary. indeed lead to increase c (see formula in Table 1) and a subsequently increased spreading of the wave.

Dam-Break

The urban dam break problem was introduced in [START_REF] Guinot | Multiple porosity shallow water models for macroscopic modelling of urban oods[END_REF]. It is an initial value problem with a piecewise constant initial state solved over the periodic street network presented in Sections 3.1 and 3.2. It is a

Riemann problem for the porosity-based equations. The initial conditions are dened as

h (x, 0) = h L if x < 0 h R if x > 0 u (x, 0) = 0∀x (26)
In the present test the street network extends from x = -2km to x = +2km. The following simulations are carried out (Table 4): (i) classical 2D shallow water equations solved using the same rened grid as in Sections 3.1-2, (ii) single porosity (SP) model [START_REF] Guinot | Flux and source term discretisation in twodimensional shallow water models with porosity on unstructured grids[END_REF], (iii) IP model, (iv) DIP model with µ xx = 0.48, (v) DIP model with µ xx = 0.41. For µ xx = 0.48 the propagation speed of the shock wave in the DIP model is the propagation speed of the middle of the shock in Simulation (i). For µ xx = 0.41 the propagation speed of the shock wave in the porosity model is that of the head of the shock in Simulation (i). Figure 7 shows the water depth and unit discharge proles computed at t = 200s. In Figure 7a the computational results of Simulation (i) are plotted for all of the computational cells.

In Figures 8b-d, the results of Simulation (i) are averaged over the 100m spatial period of the street network to produce a solution at the resolution of the porosity model, which lters out high frequency uctuations associated with wave reections o building walls. Quite expectedly, the SP model yields strongly overestimated wave propagation speeds (Figure 8a) and unit discharges (Figure 7b). The IP model produces strongly underestimated wave propagation speeds (Figure 7a) and overestimated unit discharges (Figure 7b). The DIP model, in contrast, accurately computes all aspects of the water depth and unit discharge proles except the steepness of the forward moving shock wave (Figures 7b,d).

Modifying the momentum dissipation coecient µ xx leads to dierent water depths and unit discharges behind the shock wave and does not aect the upstream part of the rarefaction wave, where no momentum dissipation occurs. 4 Field-Scale Levee-Break Flood Modelling

(i) (ii) (iii) (iv) (v

Simulation objectives and setup

The purpose of this test case is validate eld-scale applicability of the DIP model and to analyse and discuss model performance. The test case chosen here involves a hypothetical scenario whereby an earthen levee protecting a 450m×650m neighbourhood of West Sacramento (CA) suddenly fails, as

shown in Figure 8. The Sacramento conurbation is one of many metropolitan areas around the world where vast areas of development are below sea level, protected by levees, and at risk of major damage and casualties in the unfortunate event of a levee failure e.g. [START_REF] Ludy | Flood risk perception in lands protected by 100-year levees[END_REF]. Consequently, models that eciently and accurately model the hazard are of great interest for emergency preparedness and response. For this hypothetical scenario, it is assumed that water in the Sacramento River Deep Water Ship Channel is at a stage of 8m compared to an average ground elevation of 4m East of the levee, and for simplicity, that the levee breaches instantaneously with a width of 100m and a bottom elevation of 5 m. Flooding dynamics are established by solving the classical shallow-water equations on a rened grid, and the IP and DIP models are run on specially designed coarse grids and compared to the rened grid solution to characterize porosity model errors [START_REF] Kim | Urban ood modeling with porous shallow-water equations: A case study of model errors in the presence of anisotropic porosity[END_REF]. Since the purpose is to analyse porosity model errors, as opposed to structural model errors associated with shallow-water models [START_REF] Kim | Urban ood modeling with porous shallow-water equations: A case study of model errors in the presence of anisotropic porosity[END_REF], a zero bottom friction is assumed in all models. Previous work has shown that in densely developed areas, the blockage eects of buildings are far more important than bottom friction relative to ooding predictions [START_REF] Schubert | Building treatments for urban ood inundation models and implications for predictive skill and modeling eciency[END_REF][START_REF] Kim | Urban ood modeling with porous shallow-water equations: A case study of model errors in the presence of anisotropic porosity[END_REF].

0 100 200 300 400 500m

Levee failure based surveys to properly document [START_REF] Gallien | Predicting tidal ooding of urbanized embayments: A modeling framework and data requirements[END_REF]. However, for the present study, errors in the approximation of building shape (as well as storage and conveyance) do not undermine the investigation because the rened solution by the classical shallow-water equations and the coarse grid solution by the porosity models reference the same building data.

Models

Fine grid model. The classical shallow-water equations are solved on a rened grid consisting of an unstructured mesh of triangular and quadrilateral cells that is constrained by building walls. This represents the so-called building-hole method of representing buildings in urban ood models [START_REF] Schubert | Building treatments for urban ood inundation models and implications for predictive skill and modeling eciency[END_REF].

The mesh consists of 78840 cells with an average cell area 3.6m 2 within the neighbourhood and 30m 2 outside. The average cell width within the built area is 2.5m, as shown in Figure 9. Three types of boundary conditions are used in the model (Figure 9). Boundary type 1 is a wall boundary condition which is prescribed along the levee and the canal bank, as well as all building walls. Boundary type 2 is a prescribed free surface elevation which is specied to be 8m across both the upstream and Coarse grid models. A coarse grid was created to support application of the IP and DIP models.

The grid consists of 2529 cells (Figure 10) with an average area of 345 m 2 within the neighbourhood (to be compared with the average 3.6m 2 in the rened model). A combination of triangular and quadrilateral cell types were used to most eciently span the domain and place vertices at the centroids of buildings, so model edges resolve the most narrow gaps between buildings to regulate conveyance [START_REF] Sanders | Integral formulation of shallow water equations with anisotropic porosity for urban ood modelling[END_REF][START_REF] Kim | Urban ood modeling with porous shallow-water equations: A case study of model errors in the presence of anisotropic porosity[END_REF] and to ensure that φ Γ < φ Ω . The mesh shown in Figure 10 is the result of manual processing.

The boundary conditions in this model are identical to those in the rened model. As in the ne grid model, a zero bottom friction is assumed, therefore the bottom drag coecient is set to zero and only the building drag coecient is to be calibrated.

Model parametrization requires that the porosities φ Ω and φ Γ be computed over each cell and each interface respectively. This is done by sampling the coarse model cells and interfaces. In each cell of the coarse mesh, a large number of points (N = 10, 000) is generated with a uniform density so as to cover the entire cell. For each of these points, the rened model mesh is scanned to determine whether the point belongs to a building (thus lying outside the rened model domain) or is within the rened model mesh. The domain porosity φ Ω is computed as the ratio of the number n of points lying within the rened model mesh to the total number N , φ Ω = n/N . A similar procedure is used for the cell interfaces. The order of magnitude of the imprecision in the areal and connectivity porosities is therefore 10 -4 .

The porosity models were implemented using spatially uniform drag parameters and transient momentum dissipation parameters. The former were calibrated under steady state conditions, by comparison to ne-grid model predictions, and the latter were calibrated under transient conditions.

Although the momentum dissipation model presented in this paper is not part of the original IP model, it was implemented here for the sake of model comparison. In this way, the IP and DIP models have the same number of calibration parameters. Benchmarking is thus free from model dimension-induced biases. Additional detail is provided in the following sections.

m

Model performance indicators

Several model performance indicators are used. The rst is the L-norm of the error between the porosity model and the ne grid simulation results. For a simulated variable U , the norm is dened as:

L p = i φ Ωi A i U i -U i p i φ Ωi A i 1/p , p > 0 (27)
where A i and φ Ωi are respectively the area and the areal porosity of the cell i of the coarse model mesh, U i and U i are respectively the values for U over the cell i computed by the porosity and the rened model. U i is computed by averaging U in all the cells of the rened mesh that belong to the cell i of the porosity model. Three dierent norms are tested: the L 1 , L 2 and L ∞ norms, but only L 1 norms are reported hereafter because the L 2 trends were similar and the L ∞ norm was found to be insensitive.

The second indicator is the Flood Extent Agreement F A , dened as

F A = i φ Ωi A (I) i i φ Ωi A (U ) i (28a) A (U ) i = A i if i ∈ A (U ) 0 if i / ∈ A (U ) , A (I) i = A i if i ∈ A (I) 0 if i / ∈ A (I) (28b)
where the superscripts (I) and (U ) denote cells that make up the intersection and union, respectively, of the ooded area predicted by the rened model and porosity model. A ood extent agreement F A = 1 indicates that both models give identical ooded areas (thus the best possible agreement), while F A = 0 indicates that the ooded areas computed by both models have no intersection (thus the worst possible t).

The third indicator, used for vector variables such as the ow velocity and unit discharge, is the polar histogram. The polar histogram for a vector eld w is computed as follows. The 360 degrees angular range is divided into N = 36 intervals. For each class i = 1, . . . , N the following contribution is computed

w i = j∈i A j w j , W i = w i N i=1 w i (29)
where j ∈ i is the set of cells in the model for which the direction of the ow velocity falls within the ith direction interval, A j and w j are respectively the area and the average value of the vector w in the computational cell j. The sum of the W i is therefore equal to unity and the polar histogram must be interpreted as a normalized contribution of the various ow directions.

Steady-state analysis

The classical shallow-water model and porosity model are rst run for a period of t = 2400s which establishes a steady state ow from the Deep Water Ship Channel across the urbanized development and over its outer perimeter (free overfall boundary condition). At a steady-state, the dissipation parameter µ = 0, so the role of drag coecients on model accuracy can be isolated as in previous applications of porosity models [START_REF] Schubert | Building treatments for urban ood inundation models and implications for predictive skill and modeling eciency[END_REF][START_REF] Sanders | Integral formulation of shallow water equations with anisotropic porosity for urban ood modelling[END_REF][START_REF] Kim | Urban ood modeling with porous shallow-water equations: A case study of model errors in the presence of anisotropic porosity[END_REF]. The building drag coecients in the xand y-directions for the IP and DIP models are assumed uniform over the neighbourhood and zero outside. Drag coecient values between 0 and 0.025 were tested to identify the optimal parameter value (i.e., calibration).

Figure 11 shows the maps of the L 1 norm of the water depth error as a function of the building drag coecients in the xand y-directions. Figure 11a and 11b show the L 1 norm of the water depth error for the IP and DIP models respectively. Errors for the IP model are slightly larger than for the DIP model, for example, IP errors range from ∼0.12-0.19 m while DIP errors range from ∼0.11-0.17 m. Figure 11 also shows that the sensitivity of each model to the drag coecient is relatively small, at least within the range of parameter values considered here, compared to the sensitivity of the depth predictions to the porosity model (IP vs DIP). Calibrating the drag coecient against the water depth leads to the optimal values shown in Table 5. Table 5 also shows the corresponding L 1 error norms for the ow velocity and the unit discharges. The optimal value of the drag coecient for IP and DIP yield water depth errors of 0.122 and 0.107 m respectively. From a practical perspective, depth errors of 0.10-0.15 m in a ood prediction would be considered acceptable given that aerial lidar data are typically characterized by errors of 0.10-0.015 m RMSE on hard surface such as concrete [START_REF] Sanders | Evaluation of on-line DEMs for ood inundation modeling[END_REF].

However, from a model development perspective, this error poses an opportunity for additional model improvements. The L 1 -norms of the ow velocities and unit discharge are seen to be larger for the DIP model than for the IP model. Figure 12 shows the ow elds computed by the rened, IP and DIP models at t = 2400s. Figure 13 shows the polar histogram for the steady state ow velocity and unit discharge. It is visible from

1 (h) (m) L 1 (V ) (m/s) L 1 (hV ) m 2 /s F A Optimal C b D,x Optimal C b D,

Transient analysis

The preceding modelling scenario is now analysed at t = 120s when the levee-break wave front has penetrated the neighbourhood but has yet to reach the model boundary. Under transient conditions, the DIP model requires a dissipation parameter µ > 0, and here a calibration procedure is used to estimate its value using the ne grid solution as a benchmark. It is assumed that µ is uniform within the neighbourhood and zero outside. Given building layout, the principal directions of the tensor are the xand ydirections. Therefore, µ xy = µ yx = 0 and only the coecients µ xx and µ yy are calibrated assuming a range of 0 to 1. The optimal parameter values that minimize the depth error (L 1 norm) are given in Table 6. The contrast between the IP and DIP models is striking. While the optimal xand ydissipation coecients for the IP model are small (0.0 to 0.2), the optimal ydissipation coecient is 0.7 for the DIP model, while no x-momentum dissipation is required to reach optimal calibration. This is somewhat surprising in that the geometry of the streets is roughly the same in the xand y-directions. The structure of the neighbourhood thus clearly involves two principal directions. A possible explanation is that the wave propagation properties play a key role in reproducing the ow features and that the momentum dissipation tensor is of secondary importance.

Another possible explanation is that the momentum dissipation model is incompletely parameterized (see the discussion in Section 5). This may also be seen as a conrmation of previous theoretical and experimental studies [START_REF] Guinot | Multiple porosity shallow water models for macroscopic modelling of urban oods[END_REF][START_REF] Kim | Urban ood modeling with porous shallow-water equations: A case study of model errors in the presence of anisotropic porosity[END_REF] where anisotropy was identied as as a key feature of ow dynamics.

Model Optimal parameter sets.

L 1 (h) (m) L 1 (V ) (m/s) L 1 (hV ) m 2 /
Table 3 shows that the L 1 norm for the water depth and discharge per unit width is smaller in the DIP model than in the IP model. On the other hand, the L 1 norm of the ow velocity error is signicantly larger in the DIP model than the IP model. This contrast in the performance with respect to V and hV is explained as follows: the DIP model makes a larger ow velocity error than the IP model in a limited number of cells. In these cells the water depth is very small, which contributes to a small error in terms of unit discharge. The ood extent agreement in the DIP model is also signicantly better than for the IP model. Figure 14 shows water depths at at t = 120s based on the rened model (Fig. 14a), the rened model coarsened to the scale of the porosity model (Fig. 14b), the IP model (Fig. 14c), and the DIP model with three combinations of dissipation parameters (Fig. 14d-e). This shows that an increase in µ yy reduces spreading of ood water in the ydirection and qualitatively improves the accuracy of the DIP prediction, compared with the coarsened ne grid prediction. Hence, inclusion of anisotropic dissipation parameters adds a second level of control over the spreading of ood water, in addition to the anisotropic porosity distribution resulting from the intersection of cell edges and building footprints.

Note that the IP model with the parameter set (µ xx , µ yy ) = (0, 0.7) gives L 1 (h) = 0.082m, L 1 (V ) = 0.69m s -1 and L 1 (hV ) = 0.32m 2 s -1 . The unlikely parameter set (µ xx = 0, µ yy = 0.7) for the DIP model is explained to some extent by the vector elds and the polar histograms in Figures 1516. The ow velocity eld simulated over the ne grid at t = 120s exhibits the same two preferential directions as those observed at t = 2400s. The y-direction, that is associated with small depths, contributes much less to the histogram of the unit discharge. This strongly anisotropic wave propagation behaviour is preserved after averaging the ne grid simulation results onto the coarse grid. This again is an indication that the ow propagation pattern exhibits two preferential directions on the scale of the building block. The velocity and unit discharge elds computed by the IP model give too much weight to the y-direction, while the ow velocity computed by the DIP model give too much weight to the x-direction. The polar histogram for the unit discharge in the DIP model is very similar to that of the rened model. This may explain the smaller L 1 error in the water depths. Figure 15 conrms the F A values in Table 3, in that the IP model slightly overestimates the extent of the ooded area compared to the DIP model. This may be attributed to several reasons: (i) the spatial scale of the perturbation (here the width of the dyke breach) is of the same order of magnitude as the size of the building block, which is practically relevant, but violates the scale separation principle (a basic requirement for upscaled, continuum-based model validity, [START_REF] Auriault | Eective macroscopic description for heat conduction in periodic composites[END_REF][START_REF] Auriault | Heterogeneous medium. Is an equivalent macroscopic description possible?[END_REF]), (ii) the impinging wave propagates along one of the principal directions of the drag and momentum dissipation tensors, thus resulting in a poorly conditioned inverse problem, (iii) the momentum dissipation tensor model may be incomplete. In the present model, the tensor µ is assumed independent from the ow variables (water depth and ow velocities). It is easy to imagine, however, that the dissipation rate of the hydraulic transients described in Section 2 depends on the ow conditions. Besides, it may well be that the 2 × 2 tensor formulation [START_REF] Gallegos | Two-dimensional, high-resolution modeling of urban dam-break ooding: a case study of Baldwin Hills, California[END_REF] does not suce to fully reproduce the ood dynamics. This point is detailed in the discussion section.

Discussion

There has been a debate as to whether the dierential form of the porous shallow water equations is meaningful, even though it is clearly valuable for demonstrating wave propagation properties [START_REF] Sanders | Integral formulation of shallow water equations with anisotropic porosity for urban ood modelling[END_REF][START_REF] Guinot | Multiple porosity shallow water models for macroscopic modelling of urban oods[END_REF].

Writing the equations in dierential form and inferring a continuous solution should not be considered a problem if the modeller does not try to interpret the solution at a scale smaller than the typical dimension of the heterogeneities (typically, the building size). Groundwater ow modellers deal with a similar problem [START_REF] Wen | Upscaling hydraulic conductivities in heterogeneous media: an overview[END_REF] in strongly heterogeneous aquifers where the REV does not exist, or, if it does, it is signicantly larger than the computational cell size [START_REF] Baghbanan | Hydraulic properties of fractured rock masses with correlated length and aperture[END_REF]. It is possible to rene shallow-water porosity model grids to minimize numerical truncation errors. However, mesh renement undermines the primary goal of porosity models to radically reduce the computational cost of urban ood modelling, and calls for new methods of estimating porosity parameter distributions based on the intersection of mesh cells and edges with building footprints. The authors acknowledge that the implications of mesh design in terms of numerical accuracy have not been fully explored for the IP and DIP model and that further research work is needed along this line.

The unusual transient momentum dissipation model presented in Section 2.3 is clearly superior to energy dissipation models usually found in the literature. In particular, it allows the large scale selfsimilar character of solutions to be preserved, while energy dissipation models do not. Nevertheless, the computational examples shown in Sections 3.1 and 3.3 indicate that further eort is needed to make the parameterization of this model fully explicit. Obviously, the dissipation tensor depends not only on the geometry, but also on the hydraulic eld. Indeed, while the computational examples in Sections 3.1 and 3.3 use the same geometry, the optimum value is µ xx = 0.156 in Section 3.1 and µ xx = 0.48 in Section 3.3. Likewise, a strongly anisotropic parameter set (µ xx = 0, µ yy = 0.7) is obtained in the eld scale test case (Section 4) over a model with almost identical xand y-geometries. More research will be needed to characterize the relationship between the dissipation tensor and the hydraulic variables.

Another path for further research should consist in exploring a wider range of urban geometries (in terms of building alignment, spacing, street orientation, etc.) than those explored in the present paper.

The DIP and IP models include parameters that cannot be measured and therefore must be inferred.

Whereas both models require a drag parametrization for buildings, the DIP also introduces momentum dissipation parameters that allow for additional control over ood spreading. If we assume that in any application of porosity models, a ne grid solution is needed to support calibration of drag coecients, then the key dierence between the DIP and original IP models is the number of parameters that need to be calibrated, and thus the number of porosity model simulations (trials or functional calls) needed to nd optimal parameter values. The two-dimensional application presented in Section 4.4 shows that neither the IP nor the DIP models achieve a satisfactory description of the anisotropic structure of the velocity distribution in the transient phase. The eld scale test presented in Section 4 shows that the directional behaviour of the street network acts on a scale at least as large as that of the block scale. Averaging the ow eld on the scale of a building block might have been expected to yield a rather isotropic velocity vector eld. This is not the case, as the polar histograms in Figure 16 indicate. At any rate, the unlikely optimal parameter sets obtained from the calibration of both models (zero y-drag coecient for the IP model, zero x-momentum dissipation for the DIP model) seem to point towards an inadequate model structure. This preliminary conclusion is supported by previously reported work on the consequences of attempting to calibrate overland ow models with inappropriate structures [START_REF] Guinot | Calibration of physically-based models: back to basics[END_REF]. Bearing in mind that the IP model was specically designed to address anisotropy issues [START_REF] Sanders | Integral formulation of shallow water equations with anisotropic porosity for urban ood modelling[END_REF], the reason for this limitation of the IP and DIP models is not clear at this stage.

Field studies with measurements of ood extent, ood depth and ood velocities are needed to deepen understanding of the performance of urban ood models in general, and porosity models in particular. Sensitivity analyses of the one-and two-dimensional shallow water equations [START_REF] Guinot | Sensitivity equations for the 1D shallow water equations. Practical application to model calibration[END_REF][START_REF] Guinot | Sensitivity analysis of 2D steady state shallow water ow. Application to free surface ow model calibration[END_REF] show that one-and two-dimensional shallow water models behave very dierently with respect to model geometry, parameters and boundary conditions. It should thus not come as a surprise that the DIP model exhibits excellent performance in one-dimensional test cases while showing moderate success in two-dimensional ow simulations. It has been shown at the laboratory scale that a combination of physical measurements, ne grid model predictions, and porosity model predictions can be used to distinguish between structural model errors in the shallow-water equations, scaling errors associated with use of relatively coarse grid, and porosity model errors [START_REF] Kim | Urban ood modeling with porous shallow-water equations: A case study of model errors in the presence of anisotropic porosity[END_REF]. However, eld-scale performance is the true test of these models and in this context there are many additional sources of uncertainty beyond what is encountered at the laboratory scale such as uncertainty in topography, ow volumes, ood control infrastructure such as storm drains and pumps, and the physical behaviour of buildings such as their capacity to store and convey ood water, e.g., through doors and windows. Previous work on shallow-water models suggests that topographic errors and ood volumes are the primary controls on urban ood model accuracy [START_REF] Bates | Integrating remote sensing data with ood inundation models: how far have we got?[END_REF][START_REF] Gallegos | Two-dimensional, high-resolution modeling of urban dam-break ooding: a case study of Baldwin Hills, California[END_REF], so a key remaining issue in the development of porous shallow-water models is whether and to what extent porosity models errors presently serve to limit model accuracy.

Conclusions

The The DIP model includes a new, transient momentum dissipation term that acts directly on uxes, in contrast to conventional momentum dissipation terms (i.e., bottom friction) that appear as a sink term. Verication test cases support the use of dissipation only for positive waves dened by ∂ t h > 0.

The transient dissipation parameter(s) can be identied through a calibration procedure using a ne scale shallow-water model as a reference. In two dimensions, this is introduced as a 2 × 2 tensor.

Field scale testing shows that the DIP model can achieve smaller errors in ood depth and discharge per unit with, and better ood extent accuracy, compared to the IP model. Additionally, eld-scale testing shows that transient momentum dissipation introduces an important control for shaping the directionality of ood spreading within arrays of buildings that obstruct ow. Hence, the DIP oers greater versatility for eld-scale modelling than the IP model in the sense that there is greater control over model solutions. Nevertheless, both the DIP and IP models fail to some extent to reproduce the preferential directions of the ow eld in the transient phase. The practice of shallow-water porosity model calibration and application is at its infancy, and thus more work is needed to provide guidance for accurate and ecient modelling and to deepen understanding of the trade-os between model structure (including the number of calibration parameters), mesh design and model performance at eld scale.

Figure 2 :

 2 Figure 2: Propagation of a positive wave into a built area. Top: plan view of wave propagation and bore reection. Bottom: side view along transect [AB]. The dashed line indicates the initial water level.

Table 2 :Figure 3 :

 23 Figure 3: Propagation of a positive wave into an idealized urban area. Left: denition sketch for the street network. Right: bird's eye view of the free surface elevation computed at t = 600s over a 2.5m×2.5m computational grid.

Figure 4

 4 Figure 4 shows side views of the free surface elevation proles computed at t = 300s and t = 600s. Several simulations are carried out (Table 2): (i) the classical shallow water model on the 2.5m square grid, (ii) the IP model [29] without friction (C b D = 0), (iii) the DIP model with µ xx , C b D = (0, 0), (iv) DIP with µ xx , C b D = (0, 0.01) and (v) DIP with µ xx , C b D = (0.156, 0). Simulation (i) is considered as the reference simulation. The value C b D = 0.01 used in Simulation (iv) is the numerical value for which the propagation speed of the centre of the bore is approximated best. The value µ xx = 0.156

  (iv) DIP model, mu=0, CDb=0.01 (v) DIP model, mu=0.156(iv) DIP model, mu=0, CDb=0.01 (v) DIP model, mu=0.156, CDb=0

Figure 4 :

 4 Figure 4: Propagation of a positive wave into an idealized urban network. Side views of the free surface proles computed at t = 300s and t = 600s.

3 :

 3 , (iii) DIP with µ xx , C b D = (0.156, 0), where the momentum dissipation tensor D is activated regardless of the sign of ∂ t h, (iv) DIP with a drag parametrization for building resistance but without transient dissipation µ xx , C b D = (0, 0.01), and (v) DIP with the momentum dissipation tensor active for only ∂ t h > 0 (i.e., proposed model formulation). Note that ∂ t h ≤ 0 in this test case, so momentum dissipation has no eect on the uxes. The values µ xx = 0.156 and C b D = 0.01 are used in this test case because they appeared as the best tting values in the previous test. Note that in Simulation (i), the locally dropping water level next to the left boundary is due to boundary eects and should be disregarded. Propagation of a negative wave.Overview of simulations. (C): parameter is calibrated; (F): parameter is xed; (NA): Not Applicable because the parameter does not exist in the model. (ii) IP model, CDb=0 (iii) DIP model, mu=0.156 active, (iv) DIP model, mu=0, CDb=0.01 (v) DIP model, mu=0.156(iv) DIP model, mu=0, CDb=0.01 (v) DIP model, mu=0.156, CDb=0

Figure 5 :

 5 Figure 5: Propagation of a negative wave into an idealized urban network. Side views of the free surface proles computed at t = 300s and t = 600s.

  Reference, block-averaged (ii) Single porosity model (iii) IP model, CDb=0 (iv) DIP model, mu=0.48Reference, block-averaged (iv) DIP model, mu=0.48, CDb=0 (v) DIP model, mu=0.41, CDb=0

Figure 7 :

 7 Figure 7: Urban dam break problem. Comparison between the rened 2D shallow water simulation and the various models. (a) Water depth proles, Simulation (i) displayed for all computational points. (b) Water depth proles, Simulation (i) results are averaged over each 100m × 100m building block. (c), (d): unit discharge proles, Simulation (i) results are averaged over each 100m-long building block.

Figure 8 :

 8 Figure 8: Real-World test case. Aerial view of the West Sacramento neighbourhood. (c) Google Earth. Neighbourhood centre coordinates: lat. 38°37'15N, long. 121°34'44W.

Figure 9 :

 9 Figure 9: Rened model. Model denition sketch and magnied view of the mesh. Boundary type 1: impervious. Boundary type 2: prescribed water level. Boundary type 3: free outow.

Figure 10 :

 10 Figure 10: Porosity model. View of the mesh.

Figure 11 :

 11 Figure 11: Calibration of the building drag coecient. L1-norm of the water depth error. Left: IP model. Right: DIP model.

L

Figures 12 and 13

 13 Figures 12 and 13 that the rened ow simulation exhibits two preferential directions along the main street directions. The polar histograms for the velocity and the unit discharge have very similar shapes. This indicates that the ow is aligned along the two preferential directions regardless of the water depth. Quite strikingly, averaging the rened simulation results onto the coarse grids used by the IP and DIP models does not eliminate this directional behaviour. This indicates that the ow eld is strongly directional on the scale of the building block. Both the IP and DIP models fail to identify these preferential directions to a large extent. In the IP model, the histogram spreads rather evenly over the fourth quadrant (spanning the South to East directions), while in the DIP model the contribution of the x-velocity is exaggerated. The cells with a large x-velocity have small depths, therefore their contribution to the unit discharge histogram is comparatively small. This explains that the unit discharge histogram (Figure13, lower right diagram) is more isotropic than the velocity histogram. Both the IP and DIP models overestimate the velocity and unit discharge in the diagonal direction.

Figure 12 :Figure 13 :

 1213 Figure 12: Calibration of the building drag coecient. Unit discharge (hu) and velocity vector elds (u) for the steady state ow velocity and unit discharge at t = 2400s. Top: ow velocity. Bottom: unit discharge.

Figure 14 :

 14 Figure 14: Real world test case. Water depths simulated at t = 120s. xand ytick mark spacing 100 m.

Figure 15 :Figure 16 :

 1516 Figure 15: Unit discharge (hu) and velocity vector elds (u) simulated for t = 120s. Top: ow velocity. Bottom: unit discharge.

  Dual Integral Porosity (DIP) model presented herein achieves closure of the previously introduced Integral Porosity (IP) model by dening an edge-based model solution (boundary of computational cell) in addition to a cell-based model solution (interior of a computational cell). The DIP model results in a modication of the mass and momentum uxes compared to the IP model, as well as the wave propagation speeds. A wave propagation speed analysis of the DIP model shows that problem well-posedness requires that the storage porosity be greater than or equal to the conveyance porosity. Several verication test cases including a positive wave, a negative wave and a dam-break wave clearly demonstrate a major accuracy advantage of the DIP model over the IP model when predicting the ow eld, including ow depth, ow velocity and the position of shock waves even without the proposed transient dissipation model. Furthermore, the transient momentum dissipation scheme proposed for the DIP model is shown to exert a signicant inuence on the propagation properties of positive waves. This model contributes to increase the response range of the DIP model compared to the original IP model. Furthermore, the transient dissipation model introduces another level of control for improving model accuracy at the expense of an additional parameter that must be calibrated.

Table 4 :

 4 Urban dam-break problem. Overview of simulations. C: parameter is calibrated; F: parameter is xed; NA: Not Applicable because the parameter does not exist in the model.

	)

Table 5 :

 5 Calibration of the building drag coecient against the steady state solution. Optimal para-

	y

Table 6 :

 6 Calibration of the momentum dissipation coecient against the transient solution at t = 120s.

	s	F A	Optimal µ xx Optimal µ yy

  One may question the value of a porosity model ood simulation if a ne grid model is nevertheless required; benets arise when large numbers of scenarios are run to characterize model uncertainties (i.e, Monte Carlo methods) and when test cases on a relatively small domain can be used to to develop functional forms for porosity model parameters that can be extrapolated over large domains. Additional research is needed to pursue these possibilities.When both the building drag and momentum dissipation models are incorporated in both the IP and DIP models, the DIP model gives more accurate results than the IP model. It is also noted that the DIP model is more accurate than the IP model in idealized 1D tests when the transient momentum dissipation is not activated, i.e., µ = 0. This indicates that the discharge closure model proposed in Subsection 2.2 is essential to an accuracte representation of wave propagation properties.

Acknowledgements

The authors with to thank the City of West Sacramento, California, USA for graciously providing aerial lidar data for use in this study. Participation by BFS and JES on this research project was made possible by a grant from the National Science Foundation (DMS-1331611), who support is gratefully acknowledged.