Developing an effective 2-D urban flood inundation model for city emergency management based on cellular automata

Flash floods have occurred frequently in the urban areas of southern China. An effective process-oriented urban flood inundation model is urgently needed for urban storm-water and emergency management. This study develops an efficient and flexible cellular automaton (CA) model to simulate storm-water runoff and the flood inundation process during extreme storm events. The process of infiltration, inlets discharge and flow dynamics can be simulated with little preprocessing on commonly available basic urban geographic data. In this model, a set of gravitational diverging rules are implemented to govern the water flow in a rectangular template of three cells by three cells of a raster layer. The model is calibrated by one storm event and validated by another in a small urban catchment in Guangzhou of southern China. The depth of accumulated water at the catchment outlet is interpreted from street-monitoring closed-circuit television (CCTV) videos and verified by on-site survey. A good level of agreement between the simulated process and the reality is reached for both storm events. The model reproduces the changing extent and depth of flooded areas at the catchment outlet with an accuracy of 4 cm in water depth. Comparisons with a physically based 2-D model (FloodMap) show that the model is capable of effectively simulating flow dynamics. The high computational efficiency of the CA model can meet the needs of city emergency management.


Introduction
Storm-water runoff simulation models are of fundamental importance for assessing urban flooding dynamics (Chen and Adams, 2007), risk assessment (Dunn et al., 2014;Mahmoudi et al., 2013;Kubal et al., 2009;Taubenböck et al., 2011;Tsakiris, 2014) and the evaluation of runoff quantity and quality-control performance (Gilbert and Clausen, 2006;He et al., 2011).In the past decades, flash floods have occurred frequently due to land cover and climate change (Chang et al., 2010), especially in the urban areas of southern China, and have caused serious damage to urban properties.These problems present many challenges to urban management, such as storm-water pipeline maintenance, non-point source contamination, massive traffic jams.Therefore, a simple, efficient and process-considered urban flood inundation model is urgently needed for urban storm-water and emergency management.However, due to the complexity of the urban infrastructure and built environment, urban flood modeling that is sufficiently efficient for emergency management has always been a challenging task (Bidur Ghimire, 2013).
Compared to lumped hydrological models, twodimensional (2-D) hydraulic models with high spatial details based on hydraulic equations gain better accuracy and credence in simulating the spread of water onto a wide surface area that contains complex urban feature and can provide detailed information about the flooding dynamics.However, these models require a large amount of data, which are often difficult to obtain (Chen et al., 2012;Hunter et al., 2005a;Neal et al., 2009;Schumann et al., 2011;Zhang and Pan, 2014).Indeed, some cities do not have detailed information on their sewer system, including exact pipe junction locations and depths, pipe diameters and the intricate morphology and boundary conditions of channel and sewer networks (Zhang and Pan, 2014), especially in the context of ancient and emergent urban areas as a result of rapid urbanization of China.Thus, establishing a physics-based model is a difficult task for most urban management.In addition, urban flooding is a fast process at a time scale as small as tens of minutes.The calculation of the differential equations in these physics-based models is complex and requires an immense computation time.In contrast, floodplain inundation is most often characterized by a much slower phenomenon (Leandro et al., 2014); thus, the physics-based hydrologic model via calculation of the differential equations is tolerable in a floodplain inundation simulation but impractical in the urban flooding inundation modeling and emergency management at present (Versini et al., 2010).
A variety of measures have been developed for improving the model performance, including grid coarsening (Yu and Lane, 2006a), parallelization (Lamb et al., 2009;Sanders et al., 2010;Yu, 2010), adaptive grid-based methods (Hunter et al., 2005a;Wang and Liang, 2011), complexity-reduced models (Liu and Pender, 2010) and simplified governing equations (Bates et al., 2010).However, the loss of information with low resolution in grid coarsening method often leads to less accurate results, thus a subgrid treatment (Yu, 2010;Yu and Lane, 2006b) is developed to regain necessary information.Although an unstructured mesh method can effectively reduce the computing burden among these methods and has the potential to describe surface features, the preprocessing of model setup data is complex, especially in an urban area (Chen et al., 2012).Therefore, the requirements of these methods do not fully meet the need of city emergency management (Zhang and Pan, 2014).Accordingly, rapid and simplified flood inundation models (RFIM) have appeared to be common and straightforward for simulating flood inundation recently, such as RFIM (Krupka et al., 2007), rapid flood spreading model (Lhomme et al., 2008), urban inundation models based on geographic information system (GIS) (Chen et al., 2009) and flood-connected domain calculation (Zhang et al., 2014).Although these models have the advantages of computing efficiency and robustness, they only simulate the maximum inundation extent and fail to provide the flow and inundation processes.
Cellular automaton (CA) provides an alternative solution by simplifying the physics-based 2-D models, with gridbased simulation frame and allowing both time varying and geo-algebraic solution of spatial dynamics (Parsons and Fonstad, 2007).Instead of solving the differential equations iteratively within the simulation domain, the CA model describes a system on the basis of the local interactions of their constituent parts through incorporating gravitational rules to represent surface-flow runoff (Burks, 1970;Di Gregorio and Serra, 1999;Wolfram, 1984).Furthermore, CA is naturally implementable through parallel computing, allowing for efficient simulations.All these characters of CA models have the tacit accordance with the requirements of city storm-water and emergency management.While it is a sensible choice to only use the gravitational term of the 2-D shallow water equations for urban emergency management, it should be noted that the algorithms in CA models do not solve all terms in the equation necessary to replicate the full dynamics.
Since the storage cell approach was proposed by Zanobetti et al. (1970) to predict floodplain inundation, many CA models have become available in recent years for simulating braided rivers (Parsons and Fonstad, 2007;Thomas and Nicholas, 2002), river channels (Coulthard et al., 2002), floodplain rivers (LISFLOOD-FP) (Bates and DeRoo, 2000) and urban environments (Fewtrell et al., 2008(Fewtrell et al., , 2011)).These methods discretized flood plains into fine regular grids, among which the fluxes of water are calculated according to many uniform flow formulae, mostly simplified from the Saint-Venant or Manning equations (Bates et al., 2010).For instance, the LISFLOOD-FP is one of the most popular models and has been tested and improved by various scholars, including comparing with other approaches (Horritt andBates, 2001, 2002), calibrating using different data (Hunter et al., 2005a), proposing the adaptive time step (Hunter et al., 2005b) and evaluating in an urban area (Fewtrell et al., 2008(Fewtrell et al., , 2011;;Sampson et al., 2012).Bates (2010) presented a reasonable approach which considered only the inertial term for low computational cost.Thus, in some circumstances such as urban emergency management, it is pragmatic to simplify the flow-governing terms.
Therefore, this study aims to (i) develop an effective CA model for efficient simulation of storm-water runoff for city storm-water and emergency management by considering detailed urban features with little preprocessing on commonly available urban geographic data; and (ii) assess the feasibility and effectiveness of the model by using in situ observed street water depth.The study area is a small independent catchment in the downtown of Panyu district, Guangzhou, China.Development of the model is described in Sect. 2. Model calibration and validation are presented in Sect.3. Results and discussion are elaborated in Sect.4, followed by conclusions in Sect. 5.

Model description
According to the underlying theory (Wolfram, 1985a(Wolfram, , b, 1986a, b) , b) of CA models, this study develops an urban surface runoff simulation model.The model considers four parameters: CA =< L, N, P , T > . (1) L is the set of regular cells used to discretize the study area.
The model employs a 2-D lattice of regular cells to maintain consistency with digital elevation models (DEM) in space and to integrate with GIS software.The Von Neumann neighborhood (Von Neumann, 1966) N is employed, allowing water to be distributed to any of the four cells orthogonally surrounding a central cell at each iteration.At a specific time, each cell can exchange flow with all the four surrounding cells, depending on the water surface elevation.Exchanges between the central cell and the diagonal cells require two consecutive iterations of the Von Neumann neighborhood exchanges.P is the finite set of CA model parameters (Table 1) for describing the transition functions and the cell states.While cell size, elevation and infiltration threshold do not change, other state variables of the cell may be updated throughout the simulation.For instance, water is transmitted from the central cell to its neighbors at each iteration, resulting in a corresponding change of its water depth.T represents the transition rules that determine the flow routing from one cell to another.During each iteration, precipitation, infiltration and discharge to storm sewer inlets are considered within each cell.The amount and timing of flow routing from one cell to its neighbors are elaborated in Sect.2.4.
During an urban storm event, surface water flow dynamics consist of three major processes: infiltration to pervious surface, discharge to storm sewer system and surface runoff of excess water.Evapotranspiration can be neglected due to the short duration and atmospheric conditions of a storm event.Sections below present the individual transition functions associated with each process.

Infiltration
The process of infiltration is represented by a spatially varying infiltration rate for each pixel, limited by an infiltration threshold beyond which soil saturation is assumed to be reached and infiltration is no longer allowed.The infiltration threshold describes the maximum amount of surface water loss to infiltration.Different land-cover types have different infiltration rates and thresholds.Because the percentage of pervious land is very small (only 16.6 % of the study area) and the duration of the storm is relatively short, the following simple function is used to present the infiltration in the model for more efficient calibration.Once the cumulative infiltration amount is equal to or greater than the infiltration threshold, no further infiltration will occur: where R a is the infiltration rate of cells of land cover a and t is the CA time step for iteration.

Mass loss to storm sewer system
In addition to infiltration, discharge to a storm sewer system needs to be considered during an urban rainstorm event.
Inlets are typically presented along the road edges in an urban environment, through which water can discharge to urban storm sewerage systems.In this CA model, storm sewer inlets are explicitly represented by individual pixels.Instead of infiltration, the pixel that contains an inlet has a sink term associated with the mass loss to the underlying storm sewer system.Storm sewer surcharge is not considered in the model and it is assumed that the storm sewer system fully functions during a storm event.The rate of mass loss at inlets depends on a number of factors, including the local topography, storm sewer capacity and inlet configuration.In this model, the rate of discharge at the inlets to storm sewer system is determined by a method by Sun (1999) (Eqs.3 and 4): where Q is the inlet discharge flow (m 3 s −1 ), w is the area of inlet (m 2 ), C is orifice coefficient, g is gravitational acceleration, h is allowed storage of water head of an inlet, k is orifice obstruction coefficient and Ar is the area of a cell; and where D is the discharge rate of an inlet and t is the CA time step.

Surface-flow routing
After considering infiltration and inlet discharge, the rest of the storm water can form surface flow, often termed as surface runoff.An empirically derived hydrological rule known as the minimization algorithm is used in the model to represent surface runoff.This empirical method is proposed by Di Gregorio and Serra (1999) in the SCIARA model to simulate lava flow.D'Ambrosio et al. ( 2001) employed a similar method in the SCAVATU model and specified three internal transformations and two external interactions to simulate the behavior of soil erosion by water.Avolio et al. (2003) extended this method and proposed a new solution to overcome the problem of constant velocity and lack of precision.However, the specific time step is not given in the original algorithms of the three aforementioned publications.In this paper, the Manning equation is combined with the minimization algorithm to calculate the amount of water routed from a central cell to its orthogonal neighboring cells in the CA model, which is an improvement of the original methods.
In the CA model, the distribution of water from a central cell to its neighboring cells is based on an algorithm that aims to minimize the water surface elevation difference between cells.The minimization algorithm is given below.a. wel i denotes the water surface elevation of the central cell i = 0 and its neighboring cells (i = 1, 2, 3, 4).The average water surface elevation of the cells is computed by Eq. ( 6): where m is the total cell count of the central cell plus the neighboring cells involved in water redistribution -5 in step (a) -and term A refers to the set of neighboring cells that receive water from the central cell.
b. Eliminate cells with wel i > av 1 .
c. Calculate the average water elevation of the remaining cells using Eq. ( 6) and further eliminate cells where wel i > av 1 .
d. Repeat step (c) until no further cells can be eliminated.
e. Partition the outflow from the central cell such that the remaining neighboring cells have the same elevation as the central cell.
This is illustrated in Fig. 1 using a set of cells with hypothetical elevation values.
To determine the time step of the model iteration, Manning's equation is used to calculate the velocity for each pixel, defined as a function of water depth (d), surface slope (S) and roughness coefficient (n) (Eq.7).The time step is calculated for each pixel as the ratio between cell size and velocity.The minimum time step among the calculated values is used as a global time step so that stability can be maintained in the simulation.
To summarize, the combination of the minimization algorithm and Manning's equation allows for the calculation of the water elevation, the amount and direction of water exchange and time taken for the exchange.At each time step, the model applies the CA rules to each cell in the following order: (i) precipitation is added to cells according to the length of the time step and precipitation rate; (ii) infiltration or inlets discharge is subtracted to obtain the current water depth in the cell; (iii) the velocity of each cell is calculated according to Manning's equation, the traverse time of each cell is obtained and the minimum transverse time is chosen as the CA iterative time step to prevent the water from crossing a pixel in less than one time step; and (iv) the movable water is distributed to its neighboring cells following the minimization algorithm.

Case study and data availability
The model is tested in a 0.2 km 2 urban catchment (Fig. 2), the downtown of Panyu district, Guangzhou, China.While common flooding events in a natural watershed inundate large ar- eas and last for hours and days, urban flooding events usually occur in a small or even tiny subcatchment and on discrete sites with lower elevation due to natural topology or civil engineering projects.Most urban flooding events normally last for tens of minutes.The test site is one of the frequently flooded spots.It is located at a low section of a main avenue in the downtown area, and the flooding often causes serious traffic jam.There is a closed-circuit television (CCTV) camera monitoring the site and a rain gauge nearby, thus providing good rainfall data and water inundation (depth and volume) data on the street.data are critical for model calibration.
The land-cover map of the catchment is obtained from the local government, and the primary land-cover types are grassland, building, road and concrete surfaces.Around 83 % of the catchment surface is impervious.There are 96 individual storm sewer inlets in the catchment.According to Eqs. ( 3) and ( 4), the drainage rate of an inlet is calculated as 4.6 mm s −1 .All inlets have the same dimension (rectangle of 0.21 m 2 ), giving an orifice coefficient of 0.6 and obstruction coefficient of 0.67.Mark et al. (2004) suggested values between 1 and 5 m as optimal grid size to capture the main urban topographic features.The CA model is established using 5 m DEM grids.The elevation of buildings is universally raised by 10 m above the ground for the simulation purpose.
Since the resolution and accuracy of the DEM are sufficiently high, water runoff over the roads is explicitly considered in the simulation, including eventual lateral water inflows coming into the study area through the streets.As the study area is an independent catchment area, the only way for water to flow out of the system is through the 96 individual storm sewer inlets.
The model is calibrated by a storm event occurred on 19 April 2012 and validated by another storm on 4 September 2012.The first storm event lasted for 1.2 h and the peak rainfall intensity reached 91 mm h −1 .Precipitation hyeto-graph of a 1 min interval is presented in Fig. 3. Street monitoring CCTV videos are available at the outlet of the catchment for both events.Time series of water depth is estimated with uncertainties of ±4 cm based on the traffic separation metal fences and other street marks visible in the videos reasonably well.Peak depths are validated on-site after the events by measuring the watermarks left on the metal fences.Time series of the observed water depth are presented in Fig. 3a and b, depicting the inundating process and providing critical data for model calibration and validation.No surcharge was observed during the flood.

Model calibration and sensitivity analysis
Model calibration and sensitivity analysis are undertaken for the 19 April 2012 event to obtain the optimal set of parameter values for the study site.Roughness value and infiltration rate are identified as the key parameters in the CA model.Given that this is an urban site, four relatively low values of Manning's n are tested (0.01, 0.02, 0.03 and 0.04) for each of the four land-cover types (grassland, building, road and concrete surfaces), and two infiltration rates (5.8 and 11.7 mm h −1 ) are tested for grassland (Rossman, L., 2004).A total of 512 (4 × 4 × 4 × 4 × 4 × 2 = 512) simulation runs are generated in the calibration process.
Nine representative simulations are presented below to elaborate the sensitivity analysis (Table 2).They are classified into four groups.The impermeable land and paved roads account for 51.61 % of the total area.Simulations 1-4 are used to analyze the impact of the roughness of concrete surface on model predictions.Around 31.82 % of the study site is occupied by buildings.Roughness is varied in simulations 1, 5, 6 and 7 to test its impact on flow routing.Simulations 1 and 8 are used to test the impact of rough-   ness specification for grassland on flow prediction.Moreover, two infiltration values are used for the grassland surface (11.7 and 5.8 mm h −1 ) to test the model response to the infiltration parameter.Results of the sensitivity analysis are shown in Fig. 4. The model's responses to various roughness values and infiltration rates are illustrated in Fig. 4. With the increase of Manning's n for concrete surface (Fig. 4a), the simulation curve becomes smoother, and the peak water depth decreases proportionately (39.55, 38.50, 37.32 and 36.17 cm).As expected, the time of peak depth occurs later into the event (11:22,11:24,11:27 and 11:28 LT), due to the slower water movement associated with higher surface roughness.The model demonstrates a similar trend in the sensitivity analysis when the roughness for the building roof is varied (Fig. 4b).However, the magnitude is slightly smaller due to the coverage of buildings.Due to the small size of the grassland (17 %), the model is relatively insensitive to the variation of the roughness parameter (Fig. 4c).When the infiltration rate of grassland is halved from 11.7 to 5.8 mm h −1 , the peak water depth only increases by 1.6 cm and the peak depth timing is delayed by 1 min, indicating the model is relatively insensitive to the change of infiltration rate as well for the small grassland area.
Root mean square error (RMSE) is calculated by comparing the model-predicted depth with the observed depth at a 1 min interval using Eq. ( 7): where, for each simulation, q is number of the paired values and d pred i and d ref i are the predicted and reference water depths, respectively.The optimal set of parameters with a Manning n value of 0.04 (grassland), 0.01 (impermeable land and road) and 0.02 (building) and an infiltration rate of 11.7 mm h −1 is found to produce the lowest RMSE value (3.30 cm).
A series of maps (Fig. 5) is created to demonstrate the change of flooded area during the inundation process.At the beginning of the event (Fig. 5a, 10:50 LT), shallow water appeared but was constrained to the road surface.As the storm continued, water started to inundate the lower parts of the catchment (Fig. 5b and c, 10:58-11:06 LT).Thereafter, water started to accumulate around the outlet and expanded to the largest areas with maximum depth in the catchment (Fig. 5d  and e ,.At the end of the storm event, the inundated water began to decline (Fig. 5f,11:30 LT).The overall simulated pattern confirms the pattern observed from the CCTV video.

Model validation
To avoid the potential problem of model overfitting, the optimal set of parameters identified in the model calibration and sensitivity analysis for the storm event on 19 April 2012 (Fig. 6a) is used to simulate the inundation processes of another storm event on 4 September 2012 (Fig. 6b).
The hydrograph of the inundation processes simulated by the established CA model for both storm events are illustrated in Fig. 6.The CCTV videos obtained from the local government only cover 35 and 13 min for the two events.The overall patterns of water depth predicted in the two simulations agree well with the observations.For the 19 April event, the model prediction for the rising limb of the event agrees well with the observation.The model-predicted peak water depth is about 5 cm higher than the observed value from the first storm on 19 April and is in agreement with the observed values from the second storm on 4 September.It is likely that part of the difference could be attributed to the intrinsic error of the DEM model.In terms of timing, the predicted maximum water depth is reached 4 min later than the observation for the first storm, and the falling limb is much slower than the observed pattern.For the second storm, the timing of the peak value is almost identical between the observed and simulated results.
Overall, the Pearson correlation coefficient between model prediction and observation of the first storm is 0.98, with the two-tailed p value less than 0.01.Moreover, the Nash-Sutcliffe model efficiency (E) coefficient (Nash and Sutcliffe, 1970) is calculated to assess the predictability of the CA model and evaluate how consistent the predicted values match the observed.The Nash-Sutcliffe efficiency is calculated to be 0.88, indicting a good level of temporal agreement.For the second storm event (Fig. 6b), the Pearson correlation coefficient is 0.95 (p = 0.01) and the E is 0.82, also underscoring a good agreement between the predicted and observed depth values in the time series.
The results above show that the model works well for both storm events.Therefore, the model is validated and can be applied to other storm events.

Model comparison
To test the effectiveness of this CA model, a 2-D distributed model, the FloodMap (Yu, 2010;Yu and Lane, 2006a, b), is used to simulate the same storm event on 19 April 2012 in the study catchment.Three CA grids (Fig. 7a) were carefully selected to compare the inundating process simulated by both models.The three grids are located at the discharge point of this catchment (P1), on the flow routing paths along the street (P2) and on the concrete surface between build-ings (P3).The overall results simulated by both models show very good agreement (RMSEs = 0.43, 0.1, 0.02 cm, respectively).Since the two models have different discharge mechanisms, the simulated results in Fig. 7 do not consider the discharge of sewer network in order to eliminate its influence on the inundation process.
For point 1, the inundated water depths by CA are marginally higher than those of FloodMap before the peak depth of CA is reached but lower toward the end of the simulation.The inundated water depths at point 2 and 3 are much lower than point 1 due to the difference of their location in the catchment, while both models give consistent simulations.According to the above comparison, the CA model is capable of simulating the urban inundation process with similar performance to a common physically based 2-D distributed model.

Discussion
Although the simulated water depth on 19 April 2012 was in agreement with the observations in the rising limb, it fell slower than the real situation in the falling limb.This is likely related to the constant inlet discharge rate used in the CA model.Thus, a dynamic discharge rate varied with water depth could be a good choice to improve the agreement in the falling limb.The relationship between discharge rate and water depth is nonlinear, and we are collecting more observations to derive a reliable curve that can be used in the future.
In addition, DEM accuracy and grid size may have impacts on inundation simulation, especially on the catchment boundary, routing path and inundated water volume/depth.The optimal resolution of DEM for explicitly expressed topography in the CA model can be tested with various resolutions to improve the accuracy of the model.Currently, the finest available DEM grid size is 5 m, and we are implementing an airborne LiDAR (Light Detection And Ranging) data collection that can obtain finer DEM grid size of 0.5 m with higher accuracy of elevation in a much larger area of 20 km 2 .Thus, we can test the effects of different DEM grid size on the performance of the CA model and its computing efficiency.

Conclusions
This study develops an effective and efficient CA model to simulate storm-water runoff and the flood inundation process in an urban setting for city emergency management.This CA model is established primarily by using a DEM data set that considers buildings being elevated, plus the detailed urban land cover and inlet distribution along the sewerage pipe lines, while assuming that the pipeline can drain all the inflow from the inlets and ignore their possible surcharge.These two assumptions are confirmed as reasonable in the field observations.The model is calibrated and validated by the two storm events on 19 April and 4 September 2012, and simulated water depths are comparable to those observed by CCTV videos and a physically based 2-D model (FloodMap).These results suggest that the CA model is capable of generating reliable predictions to the urban storm runoff and the inundation processes, which can provide useful information for city emergency management.
This CA model is tested at a grid size of 5 m in a small urban catchment for the purpose of quick model setup and validation.At such a small catchment, the computing time on Java platform for simulating the 1.2 h storm event on 19 April 2012 is 35 s (Intel i5 CPU 3.1 GHz, 4 GB of RAM, NO parallelization).This CA model is so speed efficient precisely because it ignores important flow governing terms; nonetheless, it is pragmatic to only use the gravitational term of the 2-D shallow water equations for urban emergency management.
Future studies will examine the use of higher-resolution DEM, such as LiDAR-derived elevation data, for more accurate simulation of urban flooding.Other forms of CA models will be explored for possible improvement in efficiency.Furthermore, the model will be tested in a parallel computing environment for large metropolitan areas for the purpose of citywide storm water and emergency management.

Figure 2 .
Figure 2. (a) Location of the catchment in Panyu, Guangzhou; (b) the surface elevation DEM; (c) the components of this catchment.

Figure 3 .
Figure 3.Time series of precipitation (mm h −1 ) and observed water depth (cm) measured from the lowest point of the catchment during the (a) 19 April and (b) 4 September 2012 events.

Figure 4 .
Figure 4. Results of the sensitivity analysis: (a) various roughness values for the concrete surface; (b) various roughness values for the buildings; (c) various roughness values for the grassland; and (d) various infiltration rates for the grassland.

Figure 5 .
Figure 5. Process of flood inundation during the 19 April 2012 event.

Figure 6 .
Figure 6.Time series of observed and simulated water depth (cm), measured from the lowest point of the catchment, during the (a) 19 April and (b) 4 September 2012 events.

Figure 7 .
Figure 7. (a) Locations of selected grids and (b, c, d) temporal variation of water depth simulated by FloodMap and CA model at three selected grids.

Table 1 .
Parameters of the CA model.

Table 2 .
Representative parameter combinations used in the sensitivity analysis.