Formulation of a fast 2 D urban pluvial fl ood model using a cellular automata approach

With the increase in frequency and severity of flash flood events in major cities around the world, the infrastructure and people living in those urban areas are exposed continuously to high risk levels of pluvial flooding. The situation is likely to be exacerbated by the potential impact of future climate change. A fast flood model could be very useful for flood risk analysis. One-dimensional (1D) models provide limited information about the flow dynamics whereas two-dimensional (2D) models require substantial computational time and cost, a factor that limits their use. This paper presents an alternative approach using cellular automata (CA) for 2D modelling. The model uses regular grid cells as a discrete space for the CA setup and applies generic rules to local neighbourhood cells to simulate the spatio-temporal evolution of pluvial flooding. The proposed CA model is applied to a hypothetical terrain and a real urban area. The synchronous state updating rule and inherent nature of the proposed model contributes to a great reduction in computational time. The results are compared with a hydraulic model and good agreement is found between the two models. doi: 10.2166/hydro.2012.245 om https://iwaponline.com/jh/article-pdf/15/3/676/387056/676.pdf er 2018 Bidur Ghimire Albert S. Chen (corresponding author) Michele Guidolin Edward C. Keedwell Slobodan Djordjević Dragan A. Savić Centre for Water Systems, College of Engineering Mathematics & Physical Sciences, University of Exeter, Exeter, EX4 4QF, UK E-mail: A.S.Chen@exeter.ac.uk

shows that the space between the water levels of cells is divided into four layers.L i is the free space between the water levels of cells ranked i and i þ 1 that can accommodate the water volume from cells with higher ranks.If the rank of the central cell is r c (e.g.rank 3 as shown in Figure 1 (b)), there will be r c -1 number of cells receiving water as flux through the NH cell boundaries.In the distribution among layers, there will be at most r c -1 layers to receive water, if enough water is available in the central cell.Thus outflow volume to the layer i can be given by the following formula that is applied locally for each cell considered: where V c is the water volume of the central cell in the previous time step; ΔV k is the volume distributed to layer k, P iÀ1 k¼1 ΔV k the total volume has been distributed to layers 1 to i À 1; V c À P iÀ1 k¼1 ΔV k represents the remaining volume available for distributing to layer i after filling i-1 layers.P i k¼1 A k is the total surface area of layer i; ΔWL i is the water level difference between cells ranked i and i þ 1; A k is the available space for storage in layer i.For the layer adjacent to the central cell, an additional term ΔV k ! is applied to limit ΔV i , which assumes that the water levels for all cells will reach an equivalent level.Thus, a cell with rank r receives water only from cells with higher ranks and the water received is added on top of its own water level.Thus the total outflow flux from the central cell to neighbouring cell with rank i is calculated as: For a regular grid, the areas of the central cell, A c and the neighbouring cells, A k , are constant over the domain.
However, the methodology is applicable to different grid settings.Therefore, a cell containing buildings that do not allow water to flow in can be described using a variable cell area to reflect the reduced space occupied by buildings.

Depth updating
A very important step in the CA approach is the execution of the state transition rule.In the present CA calculations, the global continuous state is the flow depth in a grid cell, which is updated for every new time step.This is done by algebraically summing the water depth from all its four neighbours.
The following transition rule is used to update the flow depth: where, θ is a non-dimensional flow relaxation parameter that can take values between 0 and 1, F is the total volume transferred to the cell under consideration as calculated from Equation (2) and A is the cell area.The purpose of the relaxation parameter is to damp oscillations that would appear otherwise.The effect of the relaxation parameter does not impart any effect on mass conservation rather it makes the flow smooth and gradual.The values of θ are determined by numerical experiments and calibration.

Time-step calculation
For most 2D hydraulic modelling, higher resolution DEM data are being used, the required time steps will be shorter to ensure the stability of model computations, which often leads to large computational burden, such that many studies have been focused on reducing the computational time of simulations.The time increment, determined as the largest that satisfies the stability criteria anywhere in the whole domain, implies that for most of the cells only a fraction of the locally allowable time steps is used to integrate the solution in time.This represents a waste of where, d* is the water depth of flow available at the interface, which is the difference between higher water level and higher ground elevation of the central cell and its neighbour cell to the interface where, WL and z are the water level and ground elevation respectively and the subscripts C and N represent central and neighbouring cells respectively.
To prevent the velocity from overshooting, a cap on the local allowable velocity is applied as given by Equation ( 6) based on the Manning's formula and critical flow condition as: where, the hydraulic radius R is taken to be equal to the water depth d and S is the slope of water surface elevation and is always positive for outflow calculation.If v is less than v*, the interfacial flux F is recalculated by replacing v* with v in Equation ( 4).
The global time step is then calculated based on the global maximum velocity to satisfy the conventional CFL criteria.Therefore, each time the state transition rule is applied, the global time step is updated using maximum velocity calculated from all cell interfaces, as given by where v j is the velocity calculated for jth cell interface for the entire domain.

MODEL TEST CASES
The model we developed has been applied to two case studies: a hypothetical terrain case and a real-world case study in the area of Keighley, UK.The hypothetical terrain consists of 30 by 20 cells with a cell resolution of 50 m, whereas the Keighley area consists of 377 by 269 cells with a cell resolution of 2 m.The excess rainfall was directly applied over the surface of the catchment uniformly in both case studies.

Hypothetical terrain
The model was tested on a hypothetical terrain shown in  were reproduced in a manner one would expect.It can also be observed that the water exchange between the two pond areas took place through the narrow strip on the lower side of the area.
The water depth hydrographs computed by both UIM and CA models at various check points are compared in

Keighley area
The overall results obtained from the UIM and the CA model for maximum inundation depth show very good agreement (Figure 6).Higher inundation depths can be observed along the low lying areas, which was mainly due to the accumulation of water from higher ground.The  the UIM required 98 min of processing time on a desktop computer (Intel Pentium D CPU 3 GHz, 1 GB of RAM).
The root mean square error (RMSE) calculated for maximum depths over the domain, assuming the UIM results as a reference, was as small as 0.015 m.

CONCLUSION
The proposed 2D CA model has been applied to a hypothe-

θ
scheme is ensured by checking the Courant condition in every time step and choosing the computational time step accordingly.

Figure 1
Figure 1 | (a) Cells ordered in NH according to their ranks, L1-L4 are layers of free spaces between the water levels of two cells that are available within NH for water distribution, the numbers shown are cell ranks.In this diagram the ground level for each cell is shown in dark grey and the water level light grey, (b) an example of outflow fluxes (shown by arrows) from the central cell having rank 3 to its neighbouring cells.
computational effort (Zhang et al. ) and limits the use of the method.A spatially varying time step can increase solution accuracy and reduce computer run time (Wright et al. ; Crossley et al. ).In this implementation we use maximum permissible velocity which ensures the minimum time steps required to distribute the applied flux.The interfacial velocity v* is determined based on the flux transferred through a cell boundary given by:

Figure 2 .
Figure 2. The terrain consists of both forward and reverse slopes of 0.2%.It also has a lateral slope of 0.1% toward

Figure 5 .
Figure 5. Figure 5(a) shows the hydrographs at the point labelled as 'Pond' in Figure 2. The depth increased rapidly at the beginning (after the onset of rainfall), but when water started transferring to the downstream sub-catchment through the crest, the water depth became lower and then remained almost constant.Nonetheless, some discrepancies occurred at the later time stage due to the tendency of the CA model to release water quickly from higher elevation areas.The variation of flow depths at other points are shown in Figures 5(b)-5(d).The plots show the depths at the crest cell and its left and right cells.The CA model results at these points show a smoother variation than UIM results.This is attributed to the use of a flow relaxation parameter, which was taken as 0.7 for this implementation.The overall CA model results are in good agreement with Figure 3) were carefully selected so as to represent the different flow dynamics with respect to both space and time.The temporal variation of inundation depths at the selected points are shown in Figure 7.For Points 1 and 3 (Figures 7 (a) and 7(c)), the CA results were marginally higher than the UIM results after the peak depth occurs.The later horizontal portion of the plot indicated that the chosen sample points are local pond cells.For Points 2 and 5 (Figures 7 (b) and 7(e)), the temporal variations of depths are in good

Figure 7 |
Figure 7 | Temporal variation of depths at various points from UIM (solid line) and CA (dotted lane).
tical terrain and an urban area.Numerical results obtained were compared with those of a physically based model UIM that employs Saint Venant's equations.The depth hydrographs obtained at various check points for the hypothetical terrain show consistent results to the UIM, which demonstrates the CA model's capability to simulate the spatio-temporal evolution of runoff.It captured the dynamics of wetting and drying over the terrain well.In the case of Keighley, the results for maximum inundation depths show a good agreement between the CA and the UIM results.The comparisons of depth hydrographs at various sampling points also show a consistent spatio-temporal evolution of the inundation process.The state updating algorithm strongly contributes to the reduction of computational time due to the straightforward algorithm to distribute water in the NH, thereby avoiding the computation of a time consuming iterative solution which would have been incurred in hydraulic models based on partial differential equations.In both the cases mentioned above, the CA model run times were over 30 times shorter than that for the UIM for no significant discrepancies in the results.One can expect greater discrepancies, which requires further investigation, in high velocity flows where the momentum of flow speed governs the process.It should also be noted that the proposed CA algorithm is particularly suitable for parallelization which would increase the computational power of the model and therefore the range of applications.The CA approach will be further developed to include both overland flow and the sewer network flow.The 2D CA model will be enhanced by including hydrological losses, e.g.infiltration.The methodology will be further improved to minimize the discrepancy and an investigation will be conducted into the relaxation parameters in detail so that they can be determined automatically.The model presented here might also be extended to incorporate the influence of buildings based on the concepts of BCR and CRF introduced in Chen et al. (a).