A simple 2-D inundation model for incorporating flood damage in urban drainage planning

An urban inundation model was developed and coupled with 1-D drainage network model (EPA-SWMM5). The objective was to achieve a 1-D/2-D coupled model that is simple and fast enough to be consistently used in planning stages of urban drainage projects. The 2-D inundation model is based on a non-standard simplification of the shallow water equation, lays between diffusion-wave and full dynamic models. Simplifications were made in the process representation and numerical solving mechanisms and a depth scaled Manning coefficient was introduced to achieve stability in the cell wetting-drying process. The 2-D model is coupled with SWMM for simulation of both network flow and surcharge induced inundation. The coupling is archived by mass transfer from the network system to the 2-D system. A damage calculation block is integrated within the model code for assessing flood damage costs in optimal planning of urban drainage networks. The model is stable in dealing with complex flow conditions, and cell wetting/drying processes, as demonstrated by a number of idealised experiments. The model application is demonstrated by applying to a case study in Brazil. Correspondence to: A. Pathirana (a.pathirana@unesco-ihe.org)


Introduction
Considering the large spatial variability and dynamic nature of urban floods due to sewer overflow, reliable assessment of their impacts demand spatially distributed inundation models.There are numerous modelling products that can be used for the purpose.However, most of the available models tend to be complex both in terms of setting up and data preparation and are computationally expensive.Probably due to these facts, employment of inundation models for optimal planning of urban drainage systems by balancing costs and benefits, is not widely practiced by the professional community.For this purpose, arguably, a simple and quick prediction of flood within an acceptable limit of accuracy is much more important than going for an accurate but complex and time consuming hydrodynamic simulation.
This work presents the development of a simple 2-D inundation model for quick prediction of flood plain hydraulics in urban area.The main focus in the model development was to archive simplicity and speed, within an acceptable level of accuracy.The inundation model was then coupled with a drainage network model, US EPA's Storm Water Management Model (SWMM 5.0) (Rossman, 2004).The coupled model result can be used for the assessment of flood damage costs in optimal planning of urban drainage networks.
The next section describes the 2-D model algorithm and its numerical implementation.Then the coupling of the 2-D model with SWMM is explained.The performance of the coupled model is examined for several hypothetic conditions and its application is demonstrated using a simplified case study in Brazil.

2-D inundation model code development
The inundation model is based on shallow water equations (SWE) that describe depth averaged 2-D flow.The assumption of a 2-D flow over inundated plain as slow, shallow phenomena can reduce the reproduction of flood plain hydraulics to the minimum necessary level to achieve acceptable predictions.This may lead to local inaccuracies but it has been shown that the uncertainties of data (e.g. the specification of topography and boundary roughness) are dominant and thus influence the model results to a greater extent than those incurred through simplified (depth-averaged) mathematics (Hunter et al., 2007).

Model formulation
Simplified inundation models relay on the kinematics wave approximation which neglect the local acceleration, convective acceleration and pressure term in the shallow water equation.The non-inertia or diffusion wave model neglects the local and convective acceleration but considers the back water effect.However in this paper, based on the assumption that the velocity of the water flow in urban flood plains is small compared the other terms, the convective acceleration term in the momentum equation is ignored from the full dynamic wave equation.So the flood flow representation is somewhat unique as it lays between diffusion wave and full dynamic equation.The simplified shallow water flow governing equations shown below are derived from the original equations (Chow, 1964).Continuity equation............ ∂h ∂t + d ∂u ∂x + ∂v ∂y = q (1) Momentum along x-dir........ ∂u ∂t + g ∂h ∂x + gS f x = 0 (2) Momentum along y-dir....... ∂v ∂t + g ∂h ∂y + gS f y = 0 (3) where h is water stage over the datum, d is the depth of flow, q is source or sink per unit area, t is flow time, u and v are velocities along x-and y-direction, respectively, S f x and S f y are friction slopes along xand y-direction, respectively.The friction terms are estimated by using the empirical resistance relationship in Manning's equation: where n is the Manning's roughness coefficient.

Numerical model
Finite difference (FD) method is employed for solving the flow equations.The following equation is obtained using temporal and spatial discretisation of simplified continuity equation (Eq. 1) for depth averaged shallow water.
where t is finite difference time step, x and y are regular grid spacings in the x-and y-directions, respectively, n is time step counter and (i,j ) is grid counter.Using Eqs. ( 2) and (3) rate of change of velocities in xand y-directions are written as follows: Further simplification of Eq. ( 5) and ( 6) is done by representing g tS f in term of the dependent variable-velocity component (let g tS f x = uS x ).The velocities at time step n + 1 in x and y directions are then written as: where . Equations ( 4), ( 7) and ( 8) can be combined to determine the water stage above the datum.
This equation leads to a set of tridiagonal system of linear algebraic equations which can be solved by forward sweep and a backward substitution.(See Sect.2.1.2)

Alternating direction implicit (ADI) finite difference method
Implicit solving of shallow water equations leads to solving large sets of partial differential equations which is computationally expensive.In the case of fluvial flooding, the computational load of the implicit solutions is justified due to the possible offsetting by longer time steps.However, in case of small scale urban flooding (the subject of this paper), the modeller needs information at relatively finer time-resolution (due to the swift nature of urban floods), and therefore explicit approaches could often be more beneficial.In this paper, an alternating direction implicit finite difference (ADI) procedure is used to solve the governing equations.In spite of the postfix 'implicit' the ADI method behaves largely as an explicit method (due to implicit formulation only along one dimension at a time in the 2-D flood plain) and provide all the said advantages of the explicit approach.Of course this comes with the usual drawback of explicit approach as well, namely the limitation of time-step size imposed by a CFL condition (Caviglia and Dragani, 1996).According to Peaceman and Rachford (1955) ADI method provides high computational efficiency which requires less computing time because it involves a tridiagonal matrix.The method is based on splitting the differential equation in two parts and solves them sequentially in x-and y-directions within two half-time steps t n+1/2 and t n+1 , respectively.Simple diagrammatical representation of ADI routing procedure for one full time step is shown in Fig. 1.

Along the x-direction
From time step n to n + 1/2 the solution for water stage h n+1/2 along the x-axis is determined implicitly within a row while the value of other terms are expressed explicitly.Thus, the values for h n+1 i,j +1 and h n+1 i,j −1 are approximated from the previous time step values h n i,j +1 and h n i,j −1 respectively.Then Eq. ( 9) can be further simplified to a linear algebraic equation shown below: x-axis t-axis t n+1 t n+1/2

Y-direction routing
x-axis t-axis -direction routing X-direction Routing where +h n i,j + q t

Along the y-direction
From time step n+1/2 to n+1 the procedure is reversed so that the value of h along y-direction is determined implicitly while the other terms are expressed explicitly and then the water depths are updated from the calculated h values.
Analogues to the x-direction, in case of the y-direction the values for h n+1 i+1,j and h n+1 i−1,j are approximated from the previous half time step values h n+1/2 i+1,j and h n+1/2 i−1,j , respectively and Eq. ( 9) can be converted to linear algebraic equation shown below: Where Once the value of h is determined, the velocity components along the x-and y-directions will be updated by implementing Eqs. ( 7) and ( 8), respectively.

Courant condition
According to Caviglia and Dragani (1996) the ADI numerical scheme is not unconditionally stable when applied to the complete set of governing equations in an area with real topography and irregular boundaries.The instability arises because the routing method uses implicit values along one direction while the values in the orthogonal direction are determined explicitly from the previous time step.Thus courant limitation in time step t is implemented to ensure the necessary stability during simulation.
where t is finite difference time step, x and y are regular grid spacings in the x-and y-directions, respectively, and d is the water depth.In this inundation model the time step is defined as an input parameter and it has to be less than the courant t for each routing step.During each step of routing the minimum courant t of all cells in the 2-D domain, is calculated using Eq. ( 12) and used for verifying the input t value.

Process representation for wetting/drying cells
In this model the hydraulic parameters are represented for the mid-point of each square cell.However partially wetted cells i.e. cells without enough water to submerge all corners of the cell, the average depth are badly represented by depth at the centroid (Begnudelli and Sanders, 2006).Especially for a cell having very small water depth, the incidence of dry/wet condition will be very high and would frequently lead to model instability.In order to avoid this problem, we suggest a boundary depth value to classify wet and dry cells and depth factored Manning's coefficient for sudden drying wet cells.
Due to the fact that friction terms are factored by water depth, cells with very small depth value possibly create numerical instability of the model.Thus, threshold value of depth d min =0.0001 m is proposed for classifying submerged cells and updating hydraulic parameters.During each full time step the water depth is checked and the velocity and friction terms of those cells with depth less than d min are set to zero value to avoid instability.In addition to the above boundary condition, depth scaled factor for the Manning's coefficient, n modified =n(d min /d) is used to represent the reality of drying process in flood flow by slowing sudden drying of wet cells.The modification of Manning's coefficient is done for wet-cells with water depth below d min and suddenly dries within a time step.The proposed wetting/drying process representations avoid frequent instability in model simulation.

Process representation for flow barrier
Features like roads, buildings and dykes have great effect on flow dynamics and flood propagation.These barriers cause numerical instability in most of the hydrodynamic models.In this model algorithm, the flow interaction near topographic obstacles (barriers) is treated in such a way that an obstacle cell will not contribute to the velocity calculation of the nearest lower cell till it acquires a threshold depth of water d=0.01 m.
In addition, a realistic representation of flow near high barriers is achieved by modifying coefficients of the linear algebraic equations (Eqs.10 and 11) for the cells under consideration.The water stage of high barrier neighbouring is replaced with water stage value of itself ( cell under consideration) to avoid the numerical instability of the model, and to avoid spurious flows from barriers cells to any neighbouring cells.

Flood damage analysis block
Inundation models could help decision makers to have a good perception about the flood extent and to identify the effective flood mitigation measures through flood risk assessment.Providing an elementary means of estimating the loss due to flood is important for this purpose.A simple damage calculation block is integrated within the inundation model code.The maximum flood depth at each grid cell is used for the calculation of the damage in urban areas.The velocity of the water flow in urban flood plains is very small and its contribution is neglected in the calculation of flood damages.A regression equation for the possible flood damage curves is implemented in the model code.This equation is shown below.
where d is the flood depth and A to E are consecutive coefficients for the possible flood damage curves.Coupling which supports bi-directional interactions between sewer flow and 2-D inundation models is optimal to simulate sewer flooding.However, we assumed that the implementation of back flow into the sewer system is not essential, particularly for urban flood damage estimation work, where the prime concern is the maximum flood depth (not the duration as in the case, for example agricultural land).In swift urban floods, the maximum flood height is not heavily influenced by the ability of sewer network to receive return flow.Thus, coupling is achieved by mass transfer from the drainage network system to 2-D system.
There are different frameworks that can be used for model coupling.According to Bulatewicz Jr. ( 2006), they can be organised into four different categories of coupling approaches: monolithic, scheduled, communication, and component.In this paper, a unidirectional communication is used for coupling the 2-D-Inundation model with 1-D-SWMM.The unique features of this approach are: -Model codes remain independently executing programs that interact only by exchanging data via message passing during execution.
-It allows communication libraries that support direct model-to-model communication as well as model-tocoupler communication.
-It supports existing models to be coupled with minimal changes to the model source codes.

Creating a coupled model
The sewer overflow is the main interacting physical process between the two selected models.This physical process usually occurs at the point where outfall and surcharged manholes are located.The features which are integrated in the interacting physical process between 1-D and 2-D system are: -  -One of the coupling functions receives overflow information and transforms it into a flow rate per unit grid cell area within the computational domain.
-Because of the coarser time step in SWMM, the inundation model steps a number of times for a single SWMM report time step.

Coupled model algorithm
The model algorithm, which is used for surface flow routing and flood damage calculation, is written in C/C++ programming language, the choice of language largely influenced by the fact that SWMM code is available in C language.The simplified algorithm is shown in Fig. 3.

Initial and boundary conditions
The 2-D model allows setting up of several different boundary conditions (e.g.Free flow, wall, specified velocities) and initial condition (initially dry, specified flood level).However, most common initial condition for sewer overflow situation is that of initially dry condition.

Idealised experiments
The model was subjected to a number of tests using idealized situations.The flow process representations and simplifications used in this model are checked against mass conservation, stability and stationary of the result.
Several idealised test problems were utilised to examine the hydrodynamic behaviour, such as boundary conditions, 2-D flood wave diffusion, cell drying and wetting, and flow interaction with topographic obstacles.Hypothetical examples are used for quick visual inspection of the model's sensitivity and its performance against the simplifications made in the model code.
Several hypothetical tests were conducted for checking the boundary condition and the model result confirms that the numerical representation works properly by preserving the continuity and stability at boundaries.

Model performance against mass balance error
The presence of mass-balance errors in flood flow is a known problem in non-conservative numerical schemes (as opposed to conservative methods like Finite Volume approach) and improving the accuracy and mass-balance properties of fluid flow in general is an active area of research (Kippe et al., 2007).In this model the percentage of mass balance error is checked at different instants of simulation.massbal= T wdepth −T source /T source × 100 where massbal is the mass balance error in percentage , T source is the total input volume and T wdepth is the total volume on 2-D grids at same time of simulation.
The mass-balance check for different hypothetical terrains is done at different flood routing times and model outputs for flat, sloping and complex topographies are illustrated in Fig. 4, 5 and 6.The result shows the error increase with increasing complexity of the terrain (range up to 3.35 % for a very irregular topography).Inaccuracies associated with mass balance error are inevitable and the model accuracy is expected to decrease with increasing simplification in model formulation and numerical solving mechanism, with increasing time step of model integration and with increasing grid size.

Model performance in wetting and drying problems
The model was examined over hypothetical topography with the various flow conditions that may occur in actual floodplains and shown to be capable to simulate wetting and drying processes that will occur as the flood flows over an urban area.This is achieved by using a threshold value for water depth and depth scaled Manning's coefficient for wetting/drying process representation (see Sect. 2.3 ).This representation has the following advantages: -Avoids possible numerical instabilities due to very small water depth values.
-An implicit representation of the reality of high friction at very small water depths (boundary layer, laminar flow).
-Avoid sudden drying of wet cells.

Model performance in uneven topography including flow barriers
Flow barriers are often challenges for numerical stability.
The inundation maps for a hypothetical terrain involving flow barriers are presented in Fig. 8 and 9.The model performance is checked on hypothetical irregular topographies.For problems involving large flow barriers, the flood wave diffusion is visually inspected for synthetic terrains (e.g. the one shown in Fig. 7) and the probable mass balanced errors were examined at cell level.Close attention is paid for the behaviour of flood wave near barriers and other sudden changes of topography (e.g.falls, walls) and hydraulics (e.g.wetting/drying).The model works properly with global mass balance error for irregular topography (Fig. 6) of less than 3.4 percent.The result shows that the model performs well on complex topographies and frequent wetting and drying.The main drainage network runs below the street grid, except for a small part that intersects a housing block and passes under the foundations.The drainage system in the basin of Arroio da Areia can be divided into two distinct systems: one drained by gravity and the other by the pumping station Silvio Brum.The areas with a level above 8,13 m are drained by closed conduits, while the pumping station drains an area of 139,2 ha that is below the level of 8,13 m.In the end the drainage from the upstream basin flows inside a pressured pipe.
The actual capacity of the urban drainage in some parts of the Areia basin is not enough to discharge the upstream increase in flood peak and volume as a result of the urbanisation process.The heaviest inundations happen on the intersection of the roads "Nilo Peanha" and "Texeira".On this point that is the lowest of the region the drainage system, which transports the water to the "Arroio the Areia", overflows and inundation levels can reach one metre.Previous floods had resulted in damage to property and even in the loss of life.It is this general area that was used for the demonstration here.(Fig. 10)

Input data for the coupled model
Delineation of the model domain for 2-D computation can be challenging because the modeller doesn't know the region of inundation before execution (Begnudelli and Sanders, 2007).The approach adopted here was to delineate a boundary significantly far from the flooded manholes.Radar Topographic Mission (SRTM) elevation data (Farr and Kobrick, 2000) is used as the primary digital elevation model (DEM) for the area.The google-earth aerial photography was used together with site-visit experience to digitize the buildings and roads.These were superimposed on the DEM (buildings as flow barriers, roads as having a slightly lower elevation than the surrounding topography).Therefore, the buildings and roads were explicitly considered (Fig. 11) in the 2-D flow space to represent urban topography appropriately (An important requirement for small-scale urban flood modeling).Addition of buildings and roads to the elevation model introduced flow barriers and sharply bending flow paths in the topography which could lead to instabilities in the model.Thirteen nodes (physically these are man-holes) in the drainage network were connected with the inundation model.The locations based on grid counters (i,j ) along x-and ydirection are shown in Table 1.(See also Fig. 14) The model was simulated with 50 yr-2 h design rainfall event (Fig. 13).

Flood damage data
Damage data, which relate flood damage to flood inundation parameters, are different for different classes of land use and properties.This is an essential component is flood damage estimation.In most cities the central agencies have data and experience in making damage estimates but often no comprehensive guides are available at the local level.In this study residential class damage curve developed by Nascimento et al. (2006) using a survey data in the city of Itajubá, a town located in the South-eastern region of Brazil, during the year of 2002 was used as an input for completing the damage analysis (see Fig. 12).D = 130.9 + 56.3ln(d)where D is damage cost in Brazilian Real currency per sq.m and d is depth of inundation in m.

Simulation results
The 1-D pipe flow simulation results of the existing urban drainage network, show that flooding occurs at a number of nodes (shown in Fig. 14).The flood hydrographs for some of the manholes is shown in Fig. 13.There was no data for proper calibration and validation of the inundation mode, so this part was not done.Reliable inundation data are exceptionally scarce for most cities in the world, in terms of both events and spatial and temporal coverage during an event.
The coupled model simulation is performed to see the spatial and temporal variation of flood flows.Two hours of 1-D/2-D coupled simulation was done at 2-D model time step of 1 s and grid size of 20 m.The inundation model generates flood map based on maximum water depth (Fig. 15).
In addition to the flood depth, the model estimates monetary value of the damage associated with flooding, which is one of the fundamental pieces of information upon which expenditure decisions are based.For this case study the stage-damage curve developed in 2002 was used and the  model output showed the extent of flooding is 19 ha and the associated damage is about 297 000 BRL (Brazilian Real currency).

Conclusions
A new 2-D inundation model was developed and was then coupled with 1-D SWMM model to simulate surcharge induced inundation in urban areas.Due to the high computational efficiency of ADI numerical scheme, the model is fast, while maintaining stability and an acceptable level of numerical accuracy.
The 1-D/2-D coupled model was examined with a number of numerical experiments involving idealised topography, including large flow barriers, sloping terrains and irregular topography.It was found that the model performs well against stability.It has also been shown to be capable of dealing with the various flow conditions that may occur in actual urban floodplains (e.g.flowing over barriers, sharp turns), as well as being able to simulate wetting and drying processes that will occur as the flood flows over an urban area.
The model was applied for a case study to determine the detail inundation zones, depths and velocities due to surcharged water.It was also examined for some of the important hydrodynamic features, such as handling of various boundary conditions, 2-D flood wave propagation, cell drying and wetting, and flow interaction with topographic obstacles.Since the inundation model is based on detailed spatial information i.e. land use and topography, its output will be more realistic to determine flood damage.A flood damage calculation block is integrated within the inundation model code.This tool can help the planner to make quick, informed decisions on flood damages control measures on cost-benefit basis.
In planning of urban drainage a simple and fast tool for prediction of flood within an acceptable level of accuracy is often useful for the practitioner as an gentle introduction to the (currently uncommon) use of inundation models in optimal planning of drainage infrastructure.In the current proach adopted here was to delineate a boundary significan        state-of-the-art of urban storm drainage planning, there is a large polarisation of models in terms of complexity and precision.On one extreme there are simple 1-D models that do not provide any inundation information and on the other there are sophisticated 1-D/2-D (often commercial) models that sacrifice speed/simplicity to precision.There is not much middle ground!Our objective was to contribute filling this void.Sacrificing some of the precision for accessibility for a wider audience (in terms of simplicity, speed) can indeed improve current situation of under-employing inundation modelling in urban drainage planning among practitioners, particularly in the developing world.
Cost benefit analysis and optimisation of drainage infrastructure design could make significant savings and increase the reliability of drainage projects.However, optimisation problems demand large number (thousands) of simulations runs and many existing inundation models are too numerically heavy to employ for that purpose.Therefore, often these studies are done based on the volumetric flooding (as overflow from sewers), without considering the inundation effects explicitly.In the process the large errors are introduced in damage estimation.However, the current model is a good candidate for optimisation problems due to its numerically light implementation.
Edited by: N. van de Giesen

Fig. 1 .
Fig. 1.An alternating direction flood routing during each half time step.
Figure 3:1 Interacting physical systems which can be simulated by 1D-2D coup model The following functions are devised for developing a relationship between v incorporated in the main simulating engine of the 2D model.
Algorithm showing computation sequences of the 1-D-2-D coupled model (Wher ation model time step, t r and t s are the SWMM reporting and routing step).

Fig. 3 .Fig. 4 .
Fig. 3. Algorithm showing computation sequences of the 1-D/2-D coupled model (where t 2D is inundation model time step, t r and t s are the SWMM reporting and routing step, respectively).

Fig. 7 .
Fig. 7.A The synthetic terrain with a cylindrical notched flow barrier used to test the mode.
Figure 4:4 A 3D image of the hypothetical terrain data
Figure 5:1 selected flood plain area

Fig. 15 .
Fig. 15.Maximum flood depth at each grid during the 2h simulation.
Porto Alegre is the capital city of Brazil's southernmost state, Rio Grande do Sul.Porto Alegre itself has a population of 1.3 million inhabitants.It covers an area of 470 square kilometres, 40 per cent of which is urban and 60 per cent rural.The Areia basin is located in the north of the city.It covers an area of 21 km 2 , of which approximately half corresponds to the basin of Arroio da Areia and the rest belongs to the Airport polder.

Table 1 .
Location of nodes in 2-D grid domain.

Table 5 :
1. Location of nodes in 2D grid domain j