Developments of a flood inundation model based on the cellular automata approach: Testing different methods to improve model performance

Physics and Chemistry of the Earth 36 (2011) 266–280

Contents lists available at ScienceDirect

Physics and Chemistry of the Earth journal homepage: www.elsevier.com/locate/pce

Developments of a flood inundation model based on the cellular automata approach: Testing different methods to improve model performance F. Dottori ⇑, E. Todini Università di Bologna, Dipartimento di Scienze Della Terra e Geologico-Ambientali, Via Zamboni 67, 40126 Bologna, Italy

a r t i c l e

i n f o

Article history: Received 15 October 2010 Received in revised form 21 January 2011 Accepted 9 February 2011 Available online 15 February 2011 Keywords: Inundation modelling Flood propagation Cellular automata model Local time step Numerical analysis

a b s t r a c t Over the last decade, several flood inundation models based on a reduced complexity approach have been developed and successfully applied in a wide range of practical cases. In the present paper, a model based on the cellular automata approach is analyzed in detail and tested in several numerical cases, comparing the results both with analytical solutions and different hydraulic models. In order to improve the model’s performance, the original code based on the diffusive wave equations and a constant time step scheme is modified through the implementation of two techniques available in literature: an inertial formulation for the computation of discharges, originally developed for the LISFLOOD-FP model by Bates et al. (2010); and the incorporation of a local adaptive time step algorithm, based on a technique originally presented by Zhang et al. (1994). The analysis of the numerical cases showed that the proposed model can be a valuable tool for the simulation of flood inundation events. When applied to one-dimensional numerical cases, the model well reproduced the wave propagation, whereas it showed some limitations in reproducing two-dimensional flow dynamics in respect to a model based on the full shallow water equations. However, differences were found to be comparable with the uncertainty level related to available data for actual flood events. The use of the inertial formulation was very effective in all the cases, and reduced run time up to 97% as compared with the diffusive formulation, although it did not improve the overall accuracy of results. Finally, the incorporation of the local time step algorithm produced a speedup from 1.2x to 4x, depending on the simulation and the model version in use, with no loss of accuracy in the results. Ó 2011 Elsevier Ltd. All rights reserved.

1. Introduction Over the last decade, a great number of original flood inundation models have been developed, both for commercial and research purposes. The proliferation of modelling codes can be explained by considering the increasing availability of high quality topographic data (e.g. airborne laser altimetry (LIDAR) data) and the increasing computation power, which facilitates flood inundation analysis with greater degree of detail, and the simulation of flood events affecting wide areas, without excessive simulation time and data costs. In addition, the development of remote sensing imagery analysis (e.g. satellite Synthetic Aperture Radar (SAR) data) allows the measurement of the actual flood extent in a greater number of flood events, thus improving the assessment of model performances. Alongside the increasing ability to simulate flood events, there has also been considerable growth in the demand for flood inundation studies. For example, the Flood Directive 2007/60 of European Union (2007) prescribes the definition of flood hazard maps and

⇑ Corresponding author. E-mail address: [email protected] (F. Dottori). 1474-7065/$ - see front matter Ó 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.pce.2011.02.004

flood risk maps, and the development of flood management plans, in all the European river basins. Among the wide variety of modelling techniques, several approaches attempt to reproduce the dynamics of inundation events by reducing the complexity of the two-dimensional computation scheme. The solution of the full shallow water equations is therefore simplified, either through the omission of specific equation terms or through the approximation of the integration scheme, in order to capture only the essential dynamics while maintaining an adequate level of detail (Hunter et al., 2007). The main motivation for these approaches is that solving simpler equations should reduce the computational burden and therefore the simulation run times (Pender and Néelz, 2010). Actual results reported in the literature only partially agree with this hypothesis. Reduced complexity models are generally faster when applied on coarse grids (25–100 m) and wide areas (see, for example, Horritt and Bates, 2002), but run times can become prohibitive as the computation grid is refined, due to time step restrictions (Hunter et al., 2008). Recently Pender and Néelz (2010) have compared several 2D hydraulic models in different cases, and they did not observe a consistent saving in computational effort by applying simplified equations models as compared to full shallow water equations models.

267

F. Dottori, E. Todini / Physics and Chemistry of the Earth 36 (2011) 266–280

Nonetheless, their development and implementation in practical cases is relatively straightforward, and different applications to flood inundation events have shown that the use of simplified models does not necessarily imply a loss of accuracy and reliability in results (Horritt and Bates, 2002; Yu and Lane, 2006; Hunter et al., 2008; Pender and Néelz, 2010). Therefore, so far reduced complexity models have been extensively applied, and they will probably be an appropriate solution for many flood inundation analyses, especially for large scale events. Considering this framework, Dottori and Todini (2010) have presented a reduced complexity model based on the cellular automata (CA) approach and the diffusive wave equations, specifically designed to simulate flood inundation events involving wide areas. The model was tested in some numerical cases, obtaining good results both in terms of accuracy and computational speed. In the present paper, the CA model is analyzed in more detail. First, the theoretical background of the model is investigated and it is shown how the model structure can be derived from the finite volume approach under hypotheses which may be considered acceptable in most flood inundation events. Special attention is given to the discretization of head losses and different solutions, both original and taken from other hydraulic models, are presented and briefly analysed. Then the problem of the stability condition for diffusive models is discussed and two solutions for improving model performance are considered: the use of an alternative formulation for computation of fluxes, developed by Bates et al. (2010); and the incorporation of a local time step algorithm into the model, based on a technique originally presented by Zhang et al. (1994). Finally, the CA model including the proposed modifications is tested in different 1D and 2D numerical cases; the results are compared with analytical solutions and hydraulic models, either based on a reduced complexity approach or on the complete shallow water equations.

2. The original version of the CA model 2.1. The cellular automata approach and model description A cellular automaton (pl. cellular automata, abbrev. CA) is a discrete model for the description of physical dynamic systems regulated only by local laws. Typically it consists of a regular grid of cells, each in one of a finite number of states or properties, such as ‘‘on’’ and ‘‘off’’. For each cell, a set of cells called its neighbourhood (usually including the cell itself) is defined relative to the specified cell. An initial state is selected by assigning a state for each cell. The system evolves in time according to some fixed rule (generally, a mathematical function) that determines the new state of each cell in terms of the current state of the cell and the states of the cells in its neighbourhood. Typically, the rule for updating the state of cells is the same for each cell and does not change over time, and is applied to the whole grid simultaneously. Cellular automata are currently applied in computability theory, mathematics, physics, complexity science, theoretical biology and microstructure modelling, and also to fluid dynamics phenomena like debris flows (D’ambrosio et al., 2003). The model proposed in the present paper uses the macroscopic cellular automata approach for representing flood plain inundation events in which diffusion phenomena are dominant. The macroscopic CA approach partly differs from the classical CA structure so far described. Each single cell represents a volume of fluid to which the momentum and continuity equations may be applied. The cell lattice forms a grid of polygonal elements of either regular or variable shape and dimension, which represents the topography of the study area.

The adopted computational scheme can be described through the finite volume schematization. According to Garcia-Navarro and Murillo (2010), considering a single finite volume the following equation may be written:

@W þ rF ¼ S @t

ð1Þ

where W denotes the conservative variables, F the flux function, S the source term. Following a two-dimensional cell-centered finite volume scheme, Eq. (1) is integrated in a volume or grid cell X:

@ @t

Z

W dX þ

Z

X

rF dXþ ¼

X

Z

S dX

ð2Þ

X

In the CA model scheme originally proposed by Dottori and Todini (2010), W corresponds to the volume of water, V, stored in a cell, and F corresponds to the total discharge Q between cell i and the adjacent m cells, per unit width. Thus, Eq. (2) is discretized to give the solution at cell i at time level t + Dt:

V itþDt ¼ V ti þ Dt

m X

Q ti;j þ qt

ð3Þ

j¼1

where Dt is the time step in use and q is the total discharge entering or exiting the domain through cell i. The integration in time is obtained with the Euler explicit scheme, following the CA rule that the overall system state only depends on the previous time step. Considering two adjacent cells i and j, the computation of the discharge Q through the contact face is performed using the momentum equation written in the diffusive form, by neglecting the variation of velocity in time and space along the connecting direction: 5=3

Q i;j ¼

bhm n

 1=2 H i  Hj Dx

ð4Þ

where Hi, Hj are the water stages in the two cells; Dx is the distance between centroids of the two cells; b is the width of the contact face; n is the Manning roughness coefficient; hm, is the arithmetic mean between water depths in the two cells. The reasons for choosing this specific water depth value are reported in detail in Section 2.3 It is important to note that the proposed computational scheme is equivalent to a pipe network, where the nodes represent the cells, and the links represent the contact faces; hence any generic 2D flow is decomposed in 1D flow components through the links, and the momentum equation is integrated along the link direction, Dx being a vector with x and y components. Each flux component is therefore decoupled from the others, so that Eq. (4) may be solved separately for each link. 2.2. Advantages of the proposed model Although the original CA model was developed starting from the cellular automata approach, the final model structure is similar to other modelling techniques, in particular the storage cell approach. Two well-known models based on this technique are LISFLOOD-FP (Hunter et al., 2005) and FLO2D (FLO-2D Software Inc., 2007), to which the proposed model resembles under some aspects. In particular, both models perform flux computation through the decoupling of the x and y flow components: the LISFLOOD-FP model uses Eq. (4) for discharge computation, although a different equation accounting for inertial effect was recently added (see Section. 2.6 for more details), while FLO2D uses the full dynamic momentum equation. However, the CA model structure presents a number of differences and advantages which are summarized here. First, the code structure allows for the use of both regular and irregular meshes within the model. By using a specific pre-proces-

268

F. Dottori, E. Todini / Physics and Chemistry of the Earth 36 (2011) 266–280

sor for input data preparation, the CA model can directly use regular grids derived from digital elevation models (DEM), as well as polygonal grids from TIN files. The code structure is also designed to take advantage of the analogy of the computational scheme with a pipe network; all the information and the computation cycles are partitioned between node elements and link elements, in order to optimize the code efficiency. In Section 4, the effectiveness of such approach to speed up the model is analysed. In the CA model, the same computation scheme is also applied both for the 1D and 2D flows. The ability of the model to simulate one-dimensional flows has been investigated in previous tests (Dottori and Todini, 2010); in the present paper two additional 1D cases, one of them with variable mesh, are analysed in Section 4. Focusing on the flux computation between cells, different formulations for computing the friction slope along the links were considered and tested for use in the CA model, and the results are briefly reported in Section 2.3. Finally, the technique for the reproduction of drying effects is described in Section 2.4.

2.3. Discretization of head losses The choice of a specific formulation for computing the friction slope J has a great impact on model performance. In order to find the best solution in terms of accuracy and stability, different formulations have been tested within the CA model, some of which are reported in Table 1. Please note that these formulations were obtained by considering a linear variation of the water depth along the link; other formulations developed in which the hypothesis of a logarithmic variation of water depth was supposed were also tested. The formulations were tested in different numerical cases, some of which are presented in this paper, namely cases 2 and 3; the simulations included 1D and 2D flow conditions over different bed slopes. The best result in terms of correct flow propagation was obtained using the arithmetic mean hm in Eq. (4), that is, considering a linear variation of the water depth along a generic link, and computing head losses in the middle section. The other formulations were unable to correctly simulate propagation in all the flow conditions and were therefore discarded. In the LISFLOOD-FP model a similar solution is implemented, but a different reference value for water depth instead of the arithmetic mean is used. The reference value hflow is defined as the difference between the highest water free surface in the two cells and the highest water elevation (Bates et al., 2010). Such formulation was chosen by the developers of the LISFLOOD-FP model because the use of hm produced highly unstable results for all except horizontal planes, while the use of hflow did stabilize the model also in regions with complex topography (Fewtrell, personal communication). Considering these previous findings, both these formulations are tested in Section 4. Table 1 List of formulations for the computation of the friction slope J tested within the CA model. Note that all these formulations were obtained under the hypothesis of a linear variation of water depth along a generic link. Computation of head loss in the middle section of a generic link

J¼

Average of head loss between upstream and downstream sections

J¼

Integral of the head loss along a generic link

J¼

210=3 n2i;j jQ i;j jQ i;j 2

2bi;j n2i;j jQ i;j jQ i;j 2

2bi;j

ðhi þ hj Þ10=3

10=3

ðhi

10=3

Þ

7=3

Þ

þ hj

7=3

2  hj 3 ni;j jQ i;j jQ i;j ðhi 2 7 ðhj  hi Þ b i;j

2.4. Simulation of wetting and drying effects Apart from fast propagating wetting fronts, like those produced by a dam break or in the proximity of a levee breach, during an inundation event wetting fronts are usually slow, wide extended and with reduced water depths, and are therefore similar to a rainfall induced overland flow. According to different authors (Horton, 1945; Guy et al., 1992; Bagarello and Ferro, 2006), such a condition is characterized by a mixture of laminar and turbolent flow in which infiltration, surface roughness and rainfall intensity play a major role (Bateman et al., 2010). Also, the gradually varied flow assumption can be not satisfied because of perturbation produced by microtopography, and the shallow water equations can lose validity (Tayfur et al., 1993). Drying processes are also strongly influenced by these factors, and their treatment is in some ways more problematic than of the wetting front, both in terms of developing analytical and finite element solutions (Horritt, 2002). Generally, inundation models handle wetting and drying effects by using simplified relationships. In fact, it should be considered that in actual flood events the quality of experimental data is rarely sufficient to assess model results, therefore a detailed reproduction of all the mentioned processes would be scarcely useful (Horritt, 2002; Hunter et al., 2007; Neal et al., 2009a). A current solution, adopted both by reduced complexity models and more complex models, is the use of a threshold value on water depth to distinguish dry and wet cells (Jia and Wang, 2001). This simulates the surface storage effect due to terrain roughness. In addition, in presence of reduced water depth the ordinary flow equations can be linearized or substituted with an empiric relationship (Hunter et al., 2005). In the CA model, wetting fronts are reproduced by the ordinary Eq. (4), while drying effects are handled by including a control on volumes at each discharge calculation step, to ensure mass conservation and prevent negative values of water depth. In addition, a surface storage value may also be used to distinguish between dry and wet cells. It must be noted that the volume control adopted within the CA model differs from ordinary flow limiting techniques. Generally, flow limiters are used to prevent oscillations of discharge values in deep water areas, and it was observed that such technique influences model results, that become highly sensitive to grid cell size and insensitive to surface roughness (Hunter et al., 2005). On the contrary, the volume control is applied in the model only in drying conditions, that is, where water depths are very reduced. This technique has been extensively tested in cases with both regular and complex topography and no significant instabilities have been observed. 2.5. Stability problems of diffusive models For an explicit time marching hydraulic model, the maximum allowable time step to preserve stability may be obtained from the Courant–Lewy–Friedrich (CFL) condition:

Dt ¼ C Dx=ðv þ cÞ

ð5Þ

where C is a coefficient which depends on the adopted explicit algorithm and c is the celerity. However Eq. (5) only provides the upper limit value, while a model may get unstable also with lower time steps, depending on the adopted computational scheme. In particular Hunter et al. (2005) stated that explicit models not based on the fully dynamic shallow water equations need more reduced time steps. For the LISFLOOD-FP model, based on the diffusive approximation, the Von Neumann condition was considered (Hunter et al., 2005):

F. Dottori, E. Todini / Physics and Chemistry of the Earth 36 (2011) 266–280

Dt ¼

 1=2 Dx2 2ni;j @Hi;j min 5=3 4 hmi;j @x

! ð6Þ

that is, for every cell i, the term between parenthesis is computed in all the j connections with adjacent cells, and the minimum time step found is chosen to proceed with the simulation. The use of Eq. (6) in practical applications is not straightforward, since when particular flow conditions occur it may produce too low or high time step values. Hunter et al. (2005) observed that time step values tend to become very small in presence of still water surface. Therefore they introduced in the LISFLOOD-FP model a linearization of the flow equation, that is applied to areas where the water surface slope is below a specified threshold. Fewtrell et al. (2007) performed a statistical analysis of the time step distribution in an experimental application of the LISFLOODFP model, and evaluated the effect of relaxing the time step criterion given by Eq. (6). The relaxation produced some degradation in the model’s results as compared with the full adaptive time step solution. However, the difference was of the order of the error of elevation data. Previous applications of Eq. (6) within the CA model as a time step criterion gave partially different results. Dottori and Todini (2010) applied the CA model in a 2D numerical case, and observed that local time step values could be up to 100 times smaller than the mean value in computation grid. However, it was experimentally observed that the model could use much larger time steps (up to 10 times) without a degradation in the accuracy of results. In other applications of the CA model, numerical problems were observed also in presence of very low water depths, like in a wetting front, where very high values of time step were computed. Considering these previous results, a more accurate analysis of spatial distribution of time step values computed by Eq. (6) was performed and it is described in Section 4. Finally, it must be noted that, despite the possibility of a partial relaxation of the time step criterion, the CA model is still subject to the Von Neumann condition, which limits the model performance in terms of computational speed. 2.6. An alternative formulation for flux computation by Bates et al. (2010) A first solution to improve the performance of diffusive models like CA is to modify the model equations, in order to increase the maximum allowable time step up to the limit given by the CFL condition (Eq. (5)). Since its only variables are the grid size and the characteristic velocity, the CFL condition is not affected by the problems discussed in Section 2.5. Heading in this direction, Bates et al. (2010) developed an equation for discharge estimation which partly accounts for inertial effects. The derivation is herein briefly summarized: the starting point is the momentum equation from the one-dimensional De Saint Venant equations:

" # @Q @ Q2 gA@ ðh þ zÞ gn2 Q 2 þ þ þ 4=3 ¼ 0 @t @x A @x R A

ð7Þ

It should be noted that Eq. (7) was multiplied by A, which has been considered constant in time and space to be extracted from derivation signs. In so doing, Eq. (7) becomes non-conservative. Neglect2 @ Q ing the advection term @x ½ A , dividing through a constant flow width and approximating the hydraulic radius, R, with the flow depth, h, Eq. (7) can be discretized with respect to time step, Dt, to give:

q  gh @ðht þ zÞ gn2 q2t tþDt  qt  7=3 ¼ t @x Dt ht

269

ð8Þ

By rearranging Eq. (8) an explicit equation for q at time t + Dt is obtained:

qtþDt ¼ qt  ght Dt

" # @ðht þ zÞ n2 q2t þ 10=3 @x ht

ð9Þ

The same authors developed a semi-implicit version of Eq. (9) by replacing a qt in the friction term by a qt+Dt. After rearranging they obtained the following explicit form for discharge computation:

qtþDt ¼

qt  ghflow Dt @ðhþzÞ @x 10=3

1 þ ghflow Dtn2 jqt j=hflow

ð10Þ

where hflow is the difference between the highest water free surface in the two cells and the highest water elevation, as in Section 2.3. Bates et al. (2010) reported that the use of Eq. (10) in the LISFLOOD-FP model dramatically increased the model stability, with respect to the previous version based on a diffusive equation. However, the use of Eq. (10) was not always favourable in terms of accuracy with respect to the diffusive version. It must also be noted that, although Bates et al. (2010) defined Eq. (10) as an inertial formulation, such definition is not completely correct, since only the local acceleration term is included, but not the advection term. However, in order to be consistent with the previous paper by Bates et al., in the present paper the Eq. (10) is called as inertial method, following the original denomination. The structure of Eq. (10) is particularly suitable for the reduced complexity models; therefore it was introduced into the CA model, and tested in the numerical cases presented in Section 4. 3. The introduction of a local time step algorithm 3.1. Why a local time step algorithm? A further method to improve the performance of a hydraulic model is discussed in this section. The great majority of explicit time marching hydraulic models include an algorithm for modifying the current time step according to flow conditions, in order to use, during the simulation, the greatest value which guarantees stability over the whole computation grid. The time step is usually adapted considering both theoretical criteria, based on the stability analysis (e.g. the CFL condition), and empirical controls (e.g. a maximum value of water depth change in a cell during a single time step). The use of a global, adaptive time step (GTS) has a major drawback, especially in the presence of concentrated flow areas and complex topography, or when using unstructured meshes with large variations on cell size: the time step can locally assume values that are significantly smaller than values computed in the rest of computational grid, thus resulting in a considerable slowing of the simulation. A viable approach to overcome this problem is through the use of local adaptive time step (LTS) for the time marching of simulation, whereby each grid cell is integrated according to the local flow physics and numerical stability (Kleb and Batina, 1992). Early results on local time stepping have been published by Osher and Sanders (1983) for one-dimensional scalar conservation laws and, since then, many LTS methods have been proposed for numerical models solving the Euler and Navier–Stokes equations. In the field of hydraulic engineering, LTS has demonstrated utility both in 1D channel flow modelling (Crossley et al., 2003) and in 2D

270

F. Dottori, E. Todini / Physics and Chemistry of the Earth 36 (2011) 266–280

applications involving irregular topography and wetting and drying processes (Sanders, 2008). In both works, the local time step algorithm increased the computational efficiency from 25% to 70% with respect to GTS methods, with no loss in accuracy. LTS techniques were first developed and implemented in models with unstructured meshes and local grid refinement schemes. However, according to Sanders (2008) and Crossley et al. (2003), models that use a uniform grid may also benefit from LTS schemes in applications involving a variable wave speed. The implementation of a LTS algorithm within the CA model could be particularly favourable, since the Von Neumann stability condition (Eq. (6)) implies the use of locally small time steps. Hence, a local time step algorithm based on an existing method has been implemented within the CA model and is described here. 3.2. The proposed LTS formulation The implementation of a local time step is obviously not straightforward: in order to maintain a time-accurate solution, a method of accurately exchanging information between cells being integrated at disparate time steps needs to be developed (Kleb and Batina, 1992). Several methods for managing interface between regions at different time steps have been proposed in literature. Kleb and Batina (1992) presented a LTS technique based upon advancing individual cells to the level allowed by the CFL condition, then using linear interpolation at the interface to extract the information at the correct point in time. A similar approach was used by Sanders (2008), whereas Müller and Stiriba (2007) used a predictor–corrector approach to compute fluxes at the interface. All these techniques require to identify regions where cells are integrated with the same specific time step, and to define interfaces between regions where time integration has to change gradually. In addition, in order to avoid numerical errors the ratio between minimum and maximum time step needs to be limited (Crossley et al., 2003; Sanders, 2008). This results in an increase of computation cost and complication of code, although all the mentioned authors reported considerable increases in overall code speed. Considering the existing techniques and the characteristics of the CA model, the developed algorithm is based on a flux-updating procedure, similar to that proposed by Zhang et al. (1994). Time steps are computed on links between nodes, assigning to each link a time step value Dtj, which is a multiple of the smallest time step throughout the domain DT:

kj DT 6 Dtj 6 ðkj þ 1ÞDT

ð11Þ

The discharge on each link (Eq. (4) for the original CA model) is updated considering the local time step value: a computed value at a specific simulation time t is used until the simulation has progressed to t þ Dt j , then new values of time step and discharge are calculated. Temporal update of volumes in each cell (Eq. (3)) is performed considering the smallest time step DT, using the current discharge values stored for the corresponding links; thus, Eq. (3) is modified and becomes:

V i;tþDt ¼ V i;t þ

m X Dt j Q kj i;j j¼1

ð12Þ

With such procedure, the number of temporal updates is the same of a global time step method, but the number of flux evaluation, and hence computation cost and time, is considerably reduced. The proposed local time step algorithm presents differences with respect to the original version of the approach (Zhang et al.,

1994) and successive applications (Crossley et al., 2003). In particular, there is no need to divide cells in groups with different time steps, since the maximum allowable time step is directly assigned to each link and stored after each computation step. Therefore, there is no control on the difference in time step values between adjacent nodes, as it is supposed that flow conditions always generate a gradual variation in time step. Although such hypothesis may be not always true, e.g. in areas with a transition from subcritical to supercritical flow, it is reasonable since flow dynamics in most inundation events is normally slow and the Von Neumann condition imposes small time steps. The algorithm should be suitable also for the CFL condition, if adequate values of the coefficient C are chosen (Eq. (5)). The validity of such an hypothesis was tested in the numerical cases presented in Section 4. Furthermore, there is no need to limit the ratio between maximum and minimum time step during a generic simulation step. Hence, the algorithm is particularly suitable for the spatial distribution of time steps given by Eq. (6). For example, considering a levee breach simulation, the minimum time step may be less than 0.1 s near the inflow area, due to high water depths and rapidly varying flow, while the maximum value may be 10 s in areas with low water depths or still water. The only limitation introduced in the algorithm regards the minimum and maximum allowable time step (normally 0.01 s and 10–20 s), which are necessary considering the numerical problems given by the Von Neumann condition (see Section 2.5). 4. Numerical cases and discussion In the present paragraph, the original version and the proposed modifications of the CA model are tested in several numerical cases. The results are compared with analytical solutions whenever possible, or alternatively with hydraulic models based on the full shallow water equations. No experimental cases are considered since this would introduce further uncertainty sources, thus complicating the assessment of the model. Numerical cases allow to isolate single flow dynamic effects, and to assess the computation scheme under simplified and well controlled conditions, reducing the number of variables to be considered and the uncertainty concerning their influence on model performance (Hunter et al., 2005). 4.1. Case 1: horizontal channel The first benchmark for the different model versions is the numerical case proposed by Hunter et al. (2005) and applied also by Bates et al. (2010). The authors developed a one-dimensional analytical solution for the propagation of a non-breaking wave over an horizontal plane where the full De Saint–Venant equations can be simplified to yield an ordinary non-linear differential equation. The equation for the water depth h, at any point in space x and in time t is therefore:

hx;t ¼

  3=7 7 C  n2 u3 ðx  ut Þ 3

ð13Þ

where u is the component of depth-averaged velocity in the x direction and C is a constant of integration, which can be fixed by referring to the initial conditions of the problem. Here the solution was implemented using parameter values of u = 1 m s1 and n = 0.03 m1/3 s, for a simulation of duration 3600 s. The computation grid has a length of 6 km in the flow direction, and a width of 0.8 km; the grid dimensions are therefore comparable to an ordinary 2D simulation, although the flow propagation is one-dimensional due to the constant discharge along the width.

271

F. Dottori, E. Todini / Physics and Chemistry of the Earth 36 (2011) 266–280

2.5

Δx 10m. Δx 100m analytical sol.

water depth [m]

2

1.5

1

20 min

40 min

60 min

0.5

0 0

1000

2000

3000

4000

5000

x [m] Fig. 1. Case 1: comparison between the analytical solution and flow profiles computed by the diffusive version of the CA model with hm for different grid resolutions.

2.5

Δx 10m. Δx 100m analytical sol.

water depth [m]

2

1.5

1

20 min

40 min

60 min

0.5

0 0

1000

2000

3000

4000

5000

x [m] Fig. 2. Case 1: comparison between the analytical solution and flow profiles computed by the diffusive version of the CA model with hflow for different grid resolutions.

2.5

Δx 10m. Δx 100m analytical sol. 2

water depth [m]

The CA model has been tested using different grid resolutions, from 10 m to 100 m; for each grid the LTS and GTS algorithms were used in combination with the diffusive and the inertial versions. In addition, the diffusive formulation was tested using the two reference values of water depth hm and hflow, respectively adopted by CA and LISFLOOD-FP (see Section 2.3 and the variables of Eqs. (4) and (10)). An attempt to use the arithmetic mean hm within the inertial formulation was also made, but this caused strong instability and accuracy problems, so the obtained results are not represented. The results compared with the analytical solution are shown graphically in Figs. 1–3. The results of the LISFLOOD-FP inertial and diffusive versions obtained by Bates et al. (2010) are not represented in Figs. 2 and 3, since they are very similar to the analogous CA model versions using hflow. Table 2 summarizes the root mean square error (RMSE) and the minimum time steps (DT) computed for the diffusive CA model with hm, applied both with the LTS and GTS algorithms. Finally, Tables 3 and 4 summarizes the RMSE and DT values for the different versions of the CA and LISFLOOD-FP models; results for the CA model are referred to the LTS version since, as can be seen from Table 2, the LTS and GTS algorithms have the same level of accuracy and stability. As can be seen from Figs. 1–3 and the RMSE values in Tables 3 and 4, the three CA model versions produce quite different results. The original version, based on the diffusive approximation and the arithmetic mean of water depth hm, seems to be substantially insensitive to grid resolution and reproduces well the analytical solution. Differences are visible on the wave front when using the 100 m resolution, where the grid is too coarse to reproduce well the flow profile (Fig. 1). On the contrary, the water profile computed with versions using hflow is significantly influenced by the grid resolution, and the results are less accurate, especially near the wave front (Figs. 2 and 3). The inertial version produces better results than the diffusive version with hflow in coarser grids, while the opposite is true for the 10 m and 25 m grids. A similar result is observed with the analogous versions of LISFLOOD-FP, as can be seen from RMSE values in Table 4. Hence, in this particular case, the assumption of a linear variation of water surface along the link produces more accurate results, probably because of the marked differences in water depth between adjacent cells; on the contrary, the use of hflow seems to introduce some errors in the computation of head losses. A further explanation for the errors produced by the inertial formulation can be made from the hypotheses assumed for the

1.5

1

20 min

40 min

60 min

0.5

0 0

1000

2000

3000

4000

5000

x [m] Fig. 3. Case 1: comparison between the analytical solution and flow profiles computed by the inertial version of the CA model for different grid resolutions.

derivation of Eq. (10): the area A is considered as a constant in time and space between two cells, and this means that the equation becomes non-conservative, thus introducing a further error in addition to the use of hflow. It must be noted that in the numerical cases performed by Bates et al. (2010), the use of Eq. (10) was not always favourable in terms of accuracy with respect to the diffusive version, although in cases with complex topography the differences were limited. As can be seen in Tables 3 and 4, the maximum allowable time steps values are coincident for the diffusive versions of the CA and LISFLOOD-FP models, and are directly derived from the Von Neumann condition (Eq. (6)). Minor differences between the two inertial versions were observed, but this depends on the values of the coefficient C used within the CLF stability condition (Eq. (5)). For the CA model, values of C between 0.4 (for 10 m mesh resolution) and 0.5 (100 m mesh resolution) were found appropriate, while the in the LISFLOOD-FP model a value of 0.7 was adopted. Hence, the introduction of the inertial term within the CA model produces a remarkable improvement in model stability. Finally, Table 5 report run times for the different versions of the two models. Simulation times for LISFLOOD-FP model are taken from Bates et al. (2010). The simulations with the CA model were

272

F. Dottori, E. Todini / Physics and Chemistry of the Earth 36 (2011) 266–280

Table 2 case 1: minimum time step during the simulation (DT) and RMSE values for simulation steps 20, 40, 60 min for the GTS and LTS versions of CA model. Global time step algorithm

Dx (m) DT RMSE 20 min) RMSE (40 min) RMSE (60 min)

10 0.006 0.011 0.006 0.004

Local time step algorithm

25 0.037 0.012 0.006 0.005

50 0.149 0.014 0.012 0.011

100 0.602 0.052 0.04 0.034

10 0.006 0.011 0.006 0.005

25 0.037 0.013 0.007 0.006

50 0.149 0.019 0.015 0.013

100 0.602 0.051 0.037 0.031

Table 3 Case 1: minimum time step during the simulation (DT), and RMSE values at 60 min for the diffusive and inertial versions of CA model. CA diff. version with hm x (m) DT RMSE (60 min)

10 0.006 0.004

25 0.037 0.005

CA diff. version with hflow 50 0.149 0.011

100 0.602 0.034

10 0.006 0.012

25 0.037 0.033

CA inertial version

50 0.149 0.057

100 0.602 0.10

10 1.0 0.064

25 2.9 0.048

50 5.1 0.032

100 10.2 0.055

Table 4 Case 1: minimum time step during the simulation (DT), and RMSE values at 60 min for the diffusive and inertial versions of the LISFLOOD-FP model (from Bates et al., 2010). LISFLOOD-FP diffusive version

Dx (m) DT RMSE (60 min)

10 0.006 0.013

25 0.04 0.03

LISFLOOD-FP inertial version 50 0.15 0.06

100 0.6 0.09

10 1.450 0.065

25 3.62 0.05

50 7.25 0.03

100 14.51 0.05

Table 5 Total simulation times, in seconds, for CA and LISFLOOD-FP models versions with different grid resolutions (LISFLOOD-FP values are taken from Bates et al., 2010).

10 m 25 m 50 m 100 m

CA-GTS diffusive (s)

CA-LTS diffusive (s)

CA-GTS inertial (s)

CA-LTS inertial (s)

LISFLOOD diffusive (s)

LISFLOOD inertial (s)

6711 135 8 1

5275 110 7 1

177 6 1 1

135 6 1 1

1921.2 156 19.8 3

21 7.2 1.2 1.2

performed on a 1.86 Ghz Intel Core2 processor with 0.98 Gb of RAM, while simulations with the LISFLOOD-FP model were run on a single 2.66 Ghz node of a dual-core Intel Core2 Duo processor with 3 Gb of RAM. Run times for the diffusive versions of the CA model were approximately the same. As can be seen, when using a grid resolution of 25 m or coarser, the CA model versions resulted faster than the corresponding LISFLOOD-FP versions. On the contrary, the LISFLOOD-FP versions were considerably faster on the 10 m grid. The difference between the run times for this resolution is partly due to the data preparation time. For application of the inertial version to the 10 m grid, the effective run time is less than 40 s, comparable to LISFLOODFP overall run time. It must be also considered that LISFLOOD-FP simulations were run on a more powerful processor, so the comparison of run times is just qualitative. As for the LISFLOOD-FP model, the implementation of the inertial formulation produces a dramatic decrease in CA model run times, with an overall speedup of 39x for the 10 m grid. Finally, in this case the LTS technique produces a limited advantage compared to GTS version. This can be explained considering the small number of cells and the one-dimensional flow dynamics. 4.2. Case 2: channel with variable grid resolution This numerical case was presented by Dottori and Todini (2010), and simulates the one-dimensional wave propagation over a regular channel with a mild bed slope; the geometry is summarized in Table 6. The inflow is given by the following upstream hydrogram: the simulation starts from an initial condition of uniform flow along all the channel, with a constant discharge of

Table 6 Geometry and simulation settings for case 2. Bed slope

10–4

Section shape

Rectangular

Channel width Channel length

250 m 50 km

Manning roughness simulation time

0.05 m1/3 s 144 h

10 m3 s1; at t = 0, the inflow discharge from upstream is linearly increased from 10 to 100 m3 s1 in 5 h; the discharge remains constant for the next 30 h and then decrease linearly in 5 h from 100 to 10 m3 s1. The downstream boundary condition is the uniform flow. The results produced were compared with the well known model HEC-RAS (HEC, 2001). For the present simulation a computation grid made by a single row of cells and with a variable longitudinal resolution was chosen: starting from upstream, there are 50 cells of 250  250 m size, then 50 cells of 500  250 m and 50 cells of 250  250 m size. Fig. 4 shows the results in terms of flow profile, compared with HEC-RAS results. The profiles computed by the different versions of the CA model are practically coincident, and the maximum RMSE from the profiles given by the HEC-RAS model is below 0.5 cm. Such results also coincide with those obtained by Dottori and Todini (2010) using a regular mesh, thus showing the reliability of CA model when using variable resolution grids. In Table 7 a comparison of overall simulation times is presented. The simulations were performed on a 1.86 Ghz Intel Core2 processor with 0.98 Gb of RAM. Like in case 1, for the diffusive version the time steps computed by the Von Neumann condition correspond to experimental maximum values preserving simulation

water depth [m]

273

water depth [m]

water depth [m]

water depth [m]

F. Dottori, E. Todini / Physics and Chemistry of the Earth 36 (2011) 266–280

Fig. 4. Case 2: comparison between flow profiles computed by CA and HEC-RAS models.

Table 7 Total simulation times for case 2, in seconds. CA-GTS diffusive (s)

CA- LTS diffusive (s)

CA-GTS inertial (s)

CA-LTS inertial (s)

14

9

3

1

stability. For the inertial version the CFL condition with a 0.6 correction coefficient could be used. Unlike case 1, no relevant differences between the results of the model versions were observed. Focusing on computation times, the inertial formulation produce a large speed up of model. On the contrary, the advantage given by the LTS algorithm is significant but not as important. 4.3. Two-dimensional numerical cases In literature, the ability of models based on the diffusive equations to correctly simulate flooding phenomena is still under discussion (Hunter et al., 2007). To date, these models have been tested successfully against measurements from actual flood events, often performing as well as models based on complete shallow water equations (Horritt and Bates, 2002; Yu and Lane, 2006). As such, it is hypothesized that uncertainties over the data set (especially topography and roughness), influence model results to a greater extent than errors and approximations due to a simplified mathematical description (Yu and Lane, 2006; Xanthopoulos and Koutitas, 1976). Models based on diffusive equations also provided good results when tested against both analytical solutions and results from physical or alternative numerical solutions (Di Giammarco et al., 1996; Hunter et al., 2005, 2007; Xanthopoulos and Koutitas, 1976), but it must be noted that most of these works present one-dimensional analyses. In fact, 2D numerical solutions normally chosen to test hydraulic models are related to dam-break events, thus requiring models based at least on the fully dynamic shallow water equations. An alternative option to establish accurate and not-trivial benchmark cases for reduced complexity models is to test them

against models of greater complexity, in particular models based on the full shallow water equations. For example, Hunter et al. (2008) compared six models of varying complexity in a single urban flood scenario, while Pender and Néelz (2010) examined the results of several models in different benchmark test cases, including both numerical cases and actual scenarios like river floodplains and urban areas. However, according to us, more analyses are needed in order to investigate the behaviour of reduced complexity models when simulating two-dimensional flow dynamics typical of inundation events. In particular, the actual effect of the diffusive approximation on reproducing the water surface profile and the velocity field are generally not considered (the recent work by Pender and Néelz (2010) is an exception). Considering this framework, in the present paper the CA model is tested in a number of numerical 2D cases, and compared with a 2D model based on the complete shallow water equations. The chosen model for the comparison is CCHE2D, which is available free of charge on the NCCHE website. CCHE2D is a two-dimensional hydrodynamic and sediment transport model for unsteady open channel flows over loose bed, based on depth-averaged Navier–Stokes equations. The turbulent shear stress are approximated using Boussineq’s approximation and the turbulent eddy viscosity, with the use of three different closure schemes. The set of equations is solved implicitly using control volume approach and efficient element method. A detailed description can be found in Jia and Wang (2001). The comparison between the two models is carried out in three numerical cases. Case 3 is a horizontal plane made of 20  20 cells, with a 10  10 m size; case 4 is a horizontal plane as well, but the flow domain is much larger than case 3, being made of 50  50 cells with a 100  100 m size. Lastly, case 5 has the same grid size and dimension of case 4, but the plane has a 104 slope in x-direction. In all the three cases the domain is closed, except for the outlet zone, while the inflow is given by a constant discharge entering from a limited number of boundary cells (see Figs. 5, 8 and 11). The geometry and simulations settings are resumed in Table 8. It must be noted that it was not possible to set exactly the same boundary condition for the outlet within the two models.

274

F. Dottori, E. Todini / Physics and Chemistry of the Earth 36 (2011) 266–280 0

0

1

4

1.9

1

7

50

4

1.85

7

50

1.8

100

2

5

Y [m]

Y [m]

1.75

8

100

2

5

1.7

8

1.65 150

150

3 200

0

6 50

100

1.6

3

9 150

200

200

0

6 50

X [m]

100

1.55

9 150

200

X [m]

1.5

H [m]

Fig. 5. Case 3: water levels computed by CA model (left) and CCHE2D model (right) after 30 min from simulation start. Inflow and outflow areas are indicated by arrows.

point-3 0.9

0.8

0.8

water depth [m]

water depth [m]

point-1 0.9

0.7 0.6 0.5 0.4

0.6 0.5 0.4 0.3

0.3 0.2

0.7

0

10

20

0.2

30

0

10

time [min] point-5

20

30

point-9

water depth [m]

0.8

water depth [m]

30

0.7

0.9

0.7 0.6 0.5 0.4

0.6 0.5 0.4 0.3

0.3 0.2

20

time [min]

0

10

20

0.2

30

0

10

time [min]

time [min]

Fig. 6. Case 3: water depths computed in control points 1, 3, 5, 9, by the CCHE2D model (grey line) and the CA model (black line).

0

0

50

50

1.5

Y [m]

Y [m]

1 100

100 0.5

150

200

150

0

50

100

150

200

X [m]

200

0

50

100

X [m]

150

200

0

V [m/s]

Fig. 7. Case 3: velocity fields computed by CA model (left) and CCCHE2D model (right) after 30 min from simulation start.

Therefore in all the simulations local differences are visible in this area. However, tests with different boundary conditions showed that the influence of the outlet area is limited to the nearest cells and does not condition the overall propagation of flow.

As can be seen from the description, the geometry and the boundary conditions were chosen to reproduce a slow propagation inundation event. Moreover the simple geometry allows for a better comprehension of flow dynamic, and for an easier comparison

275

F. Dottori, E. Todini / Physics and Chemistry of the Earth 36 (2011) 266–280 0

0

1

4

1

1000

4

7

1000

2000

2

5

Y [m]

Y [m]

2

7

8

1.8

2000

1.6

2

5

8

3000

3000

1.4

4000

4000

1.2

5000

3 0

1000

2000

6 3000

9 4000

5000

5000

3 0

1000

2000

X [m]

6 3000

9 4000

5000

X [m]

1

H [m]

Fig. 8. Case 4: water levels computed by CA model (left) and CCCHE2D model (right) after10 h from simulation start. Inflow and outflow areas are indicated by arrows.

point-2

point-1

0.6

0.4

point-3

0.55

0.5

0.5

0.45

water depth [m]

0.8

water depth [m]

water depth [m]

1

0.45 0.4 0.35 0.3

0

2

4

6

8

0.2

10

0.35 0.3 0.25

0.25 0.2

0.4

0

2

time [h]

4

6

8

0.2

10

time [h]

0

2

4

6

8

10

time [h]

point-5 0.5

0.4 0.35 0.3

water depth [m]

water depth [m]

water depth [m]

0.45

0.25 0.2

0

2

4

6

8

10

time [h]

time [h]

time [h]

Fig. 9. Case 4: water depths computed in control points 1, 2, 3, 5, 6, 9 by the CCHE2D model (grey line) and the CA model (black line).

of models results. Since the CCHE2D model can manage only a limited number of dry cells, these simulations don’t include wetting– drying effect. Like in cases 1 and 2, the LTS and GTS algorithms were used in combination with the diffusive and the inertial versions; the diffusive formulation was also tested using both the two reference values of hm and hflow. The results are compared in terms of water surface, water depth and velocity. In order to highlight the differences between the two models, for each case the water depth evolution was plotted in different control points (Figs. 6, 9, 12 and 13). Please note that, for the CA model, velocity vectors are obtained simply by computing the vectorial sum over all the links connected to a specific cell.

4.3.1. Case 3 results The graphs in Figs. 5–7 refer to case 3. The results of the diffusive CA model version with hm and the inertial version are very similar: the maximum difference is below 1 cm for water depth values, and 0.1 m s1 for velocity values. Moreover, the inertial version and the diffusive version with hflow provide almost the same results.

As expected, the results produced by CA and CCHE2D show some differences. As can be seen from graphs in Fig. 6, the difference between the water surface computed by the two models is below 10 cm throughout the grid, going from 8–9% to 12–13% of the CCHE2D water depth. Focusing on the velocity field (Fig. 7), it can be seen in the CCHE2D results that the incoming flow produces a jet effect on the velocity distribution near the upper left boundary, because of the momentum conservation. On the other hand, the velocity distribution of CA model does not show such effect, and the flow has a clear diffusive behaviour.

4.3.2. Case 4 results As for case 3, the differences between the diffusive CA model version with hm and the inertial version in terms of water surface are quite limited: the average difference is below 0.5 cm, with local values up to 1.8 cm. Also the velocity fields are similar, with differences below 0.1 m s1 in all the grid. Again, the inertial version and the diffusive version with hflow provide almost the same results. Since the overall results of the original, diffusive version of the

276

F. Dottori, E. Todini / Physics and Chemistry of the Earth 36 (2011) 266–280 0

0

1000

1000

2000

2000

0.8

Y [m]

Y [m]

0.7 0.6 0.5 0.4

3000

3000

4000

4000

5000

5000

0.3 0.2 0.1

0

1000

2000

3000

4000

5000

0

1000

X [m]

2000

3000

4000

5000

X [m]

0

V [m/s]

Fig. 10. Case 4: velocity fields computed by CA model (left) and CCCHE2D model (right) after10 h from simulation start.

0

0

4

1000

1000

2000

2000

2

5

8

3000

3000

4000

4000

5000

3 0

1000

2000

2.5

7

Y [m]

Y [m]

1

6 3000

9 4000

5000

X [m]

5000

1

4

7

2

5

8

3 0

1000

2000

6 3000

X [m]

2

9 4000

5000

1.5

H [m]

Fig. 11. Case 5: water levels computed by CA model (left) and CCCHE2D model (right) after10 h from simulation start. Inflow and outflow areas are indicated by arrows.

Table 8 Geometry and simulation settings for cases 3–5.

Grid dimensions Cell dimension Slope Manning roughness (m1/3 s) Incoming discharge (m3 s1) Initial water depth (m) Simulation time

Case 3

Case 4

Case 5

100  100 m 10  10 m None 0.05

5000  5000 m 100  100 m None 0.05

5000  5000 m 100  100 m 10–4 0.05

20

200

200

0.2 30 min

0.2 10 h

0.2 10 h

CA model are slightly better, in the next graphs only the results of this version are displayed. As can be seen from graphs in Figs. 8–10, in CCHE2D simulation higher water depth values are computed with respect to CA results; however, the difference between the water surface computed by the two models is everywhere below 10% of the CCHE2D water depth. The velocity fields are similar (Fig. 10), although differences are visible near the inlet and outlet area, since it was not possible to set exactly the same boundary condition for the two models. 4.3.3. Case 5 results The presence of a mild terrain slope in case 5 modifies the results computed by the CA and CCHE2D models in respect to case 4. As can be seen from Figs. 12 and 13, the differences between the results can be locally up to 15 cm (points 1 and 9), therefore greater than in cases 3 and 4. Focusing on the overall distribution, CCHE2D computes lower water depths where the water surface slope is steeper, while in

the rest of the grid it computes higher water depths, accordingly with the results obtained in case 4 (Fig. 11). On the contrary, the differences between velocity fields are somewhat limited, except for the inlet and outlet areas (Fig. 14). Again, no significant difference between the different versions of the CA model was observed, both in terms of water surface and velocity field.

4.3.4. Discussion Some considerations can be made about the observed results. Since the inertial and the diffusive versions of the CA model provided similar results, the differences between the results of CA and CCHE2D models should be explained as an effect of the advection term of the momentum equation. Previous research works have analysed the relative importance of this term in one dimensional channel flow, and it has been observed that advection term can be relevant in channels with a variable cross section (Dottori et al., 2009), and when representing small scale features (Hunter et al., 2007). The results observed in the present paper suggest that the importance of the advection term is not negligible in 2D flow, even though further numerical analyses are needed to investigate its magnitude in a wider range of flow conditions. It is worth noting that the CA model maintains a diffusive behaviour also when using the inertial formulation, as can be inferred from the observation of results in case 3 (Figs. 4 and 6). This can be explained considering that in the model structure no information about the flow direction is transferred between adjacent cells, since the hydraulic head is stored in cells simply as water level. Therefore, the introduction of the inertial term within the model is very effective for stabilizing the flow along the connections, but it seems to add little further description of the overall flow dynamics, at least in the presented cases. Similar results were

277

F. Dottori, E. Todini / Physics and Chemistry of the Earth 36 (2011) 266–280 point-2 0.35

0.6

0.3

water depth [m]

water depth [m]

point-1 0.7

0.5 0.4 0.3 0.2

0.25 0.2 0.15

0

2

4

6

8

0.1

10

0

2

4

time [h]

6

8

10

8

10

time [h]

point-3

point-4

0.2

0.5

water depth [m]

water depth [m]

0.45 0.15

0.1

0.05

0.4 0.35 0.3 0.25

0

0

2

4

6

8

0.2

10

0

2

4

time [h]

6

time [h]

Fig. 12. Case 5: water depths computed in control points 1, 2, 3, 4 by the CCHE2D model (grey line) and the CA model (black line).

point−5

point−6

0.5

0.4

water depth [m]

water depth [m]

0.45 0.4 0.35 0.3 0.25

0.35 0.3 0.25 0.2

0.2 0

2

4

6

8

10

0

2

4

time [h]

8

10

8

10

point−9 0.7

0.6

0.6

water depth [m]

water depth [m]

point−7 0.7

0.5 0.4 0.3 0.2

6

time [h]

0.5 0.4 0.3

0

2

4

6

8

10

time [h]

0.2

0

2

4

6

time [h]

Fig. 13. Case 5: water depths computed in control points 5, 6, 7, 9 by the CCHE2D model (grey line) and the CA model (black line).

observed also by Bates et al. (2010) in applying the inertial and diffusive versions of the LISFLOOD-FP model to a urban flood scenario. Further comparisons of the two formulations in actual flood events are therefore needed to assess their performances. These observations can be used to draw some considerations about the effective influence of the reduced complexity of the CA model in the simulation of actual flood events. Results obtained in case 3 indicate that the CA model can reproduce locally concentrated flows only with some approximation. It must be also remembered that the CA model does not in-

clude a shock capturing numerical scheme; therefore the model would not probably be able to reproduce particular flow conditions, like supercritical to subcritical flow transitions in the proximity of a levee break, or the wake effect of buildings (Pender and Néelz, 2010). Nevertheless, the overall results in cases 3–5 suggest that these localized flow conditions may have a minor influence on overall 2D flood propagation, and previous research studies on an urban flood scenario by other authors came to a similar conclusion (Hunter et al., 2008; Pender and Néelz, 2010).

278

F. Dottori, E. Todini / Physics and Chemistry of the Earth 36 (2011) 266–280 0

0

1000

1000

2000

2000

0.8

Y [m]

Y [m]

0.7 0.6 0.5 0.4

3000

3000

4000

4000

5000

5000

0.3 0.2 0.1

0

1000

2000

3000

4000

5000

X [m]

0

1000

2000

3000

X [m]

4000

0 5000 V [m/s]

Fig. 14. Case 5: velocity fields computed by CA model (left) and CCCHE2D model (right) after10 h from simulation start.

In particular, the observed differences between the CA and the CCHE2D models are comparable with the uncertainty normally associated to available ‘‘real world’’ data for flood events. For example, according to Fewtrell et al. (2007) and Hunter et al. (2008), LiDAR data have a typical vertical error of ±15 cm, comparable with the maximum difference in water depths observed in case 5. High level of uncertainty are also associated with calibration and evaluation data, like inundation extent maps from SAR imagery, or measurements of water level from water marks (Horritt et al., 2007; Neal et al., 2009a). Therefore, if detailed localized analyses are not required, the CA model should be able to reproduce with an adequate detail 2D flow over wide areas, both in terms of water surface and velocity. A further consideration is needed about the use of hm and hflow within the equations for discharge computation (Eqs. (4) and (10)). In all the two-dimensional cases, the differences the CA model versions were somewhat lower as compared to the one dimensional case described in Section 4.1. However, it has been reported by the developers of the LISFLOOD-FP model that the use of hm results in highly unstable results in case with complex topography, while the use of hflow did stabilize the model (Fewtrell, personal communication). On the contrary, previous applications of the CA model in cases with complex topography did not record such problems: minor instabilities in water surface were observed only in areas with very low water depths and steep slopes (above 1%), but they did not affect the overall simulation results. Therefore, further comparisons will be necessary to establish the advantages and drawbacks of the two formulations. 4.3.5. Simulation times Table 9 reports the computation time for the simulations performed with CCHE2D and the CA model versions. All the simulations were run on a 1.86 Ghz Intel processor with 0.98 Gb of RAM. The CCHE2D model uses an implicit time marching scheme with a constant time step, so the maximum allowable values reported in Table 8 were chosen with a trial-and-error procedure.

Table 9 Total simulation times, in seconds, for two-dimensional cases (3–5). The results obtained by setting a minimum time step value are identified by ‘‘(lim.)’’.

CA diffusive v. GTS (s) CA diffusive v. GTS (lim.) (s) CA inertial v. GTS (s) CA diffusive v. LTS (s) CA diffusive v. LTS (lim.) (s) CA inertial v. LTS (s) CCHE2D (s)

Case 3

Case 4

Case 5

170 42 1 57 33 1 30

221 57 12 68 39 9 112

450 55 12 105 28 9 113

For the diffusive versions of the CA model, simulation times were first computed considering the Von Neumann condition alone. However, during the simulation analysis it was observed that the Von Neumann condition locally computed very low time step values, thus producing a considerable slowing of simulation (see discussion in Section 2.4). Such low value were found to be not necessary for stability: in case 3 for example, according to the Von Neumann condition minimum time step values around 5  103 s are needed; however, it was experimentally observed that a minimum value of 0.02 s was sufficient to preserve stability, and producing the same results both in terms of water surface and velocity. Therefore the CA simulations were run again introducing a minimum time step value, which obviously reduced computation times (Table 9). It must be pointed out that this finding does not allow for a generalized use of larger time steps, since the Von Neumann condition is still valid in most of the computation grids. In particular, it was observed that the velocity field is more subjected to instability than the water surface; for example, in case 3 the water surface computed with a time step limit of 0.05 s does not show oscillations, whereas relevant oscillations are visible in velocities, both in module and direction. A further research could be necessary to establish a general rule for the application of the lower limit to time step, as a function of flow condition and cell size. However, such limitations can be completely avoided if the inertial formulation is used, since the maximum allowable time step can be obtained directly from the CFL stability condition (Eq. (5)), using an appropriate value for the coefficient C. For cases 3–5 the adopted values were between 0.4 and 0.5. As for the onedimensional cases, the inertial formulation produced a dramatic speed up of the model especially in case 3 (33x) but the improvement was remarkable also for cases 4 (4.3x) and 5 (3.1x). Focusing on the comparison between the LTS and GTS algorithms, it can be seen that the diffusive LTS version is significantly faster than the GTS version, with a speed up from 2.5x to 4x. When the inertial formulation is used the run time reduction is less important although it should be noted that the numerical cases included only regular grids with a limited number of cells. A further considerations can be made about the comparison between CCHE2D and diffusive CA run times. In case 3, the CCHE2D model was significantly faster than the diffusive CA model, mainly because the computation grid includes few cells and the CA model needs to use very low time step values to preserve stability. On the contrary the diffusive CA model with a lower limit to time step value was faster in cases 4 and 5, in which the grid extension is larger. However, it must be remembered that CCHE2D is an implicit time-marching model, therefore the procedure for finding the optimum value of time step in a practical case could require running different simulations, thus being time-consuming. This is avoided

F. Dottori, E. Todini / Physics and Chemistry of the Earth 36 (2011) 266–280

with an explicit model like CA, since the maximum available time step is directly computed by the model. 5. Conclusions The analysis of the numerical cases presented in this paper have allowed us to investigate the effective ability of the CA model to reproduce two-dimensional flow dynamics. Some of the considerations herein drawn could also be applied to similar reduced complexity models like LISFLOOD-FP and FLO2D. The CA model could well reproduce the wave propagation in one-dimensional cases, and the run times were comparable to the similar and established LISFLOOD-FP model. When applied to two-dimensional numerical cases, the model showed some limitations in reproducing flow dynamics with respect to a 2D model based on full shallow water equations. Such results were partly expected since the original formulation of the model does not account for acceleration and advection effects. However, such limitations were found to be acceptable for inundation modelling, considering the uncertainty level related to available data for actual events, e.g. topography and inundation extent maps. In order to increase model performance and to reduce run times, two existing techniques were implemented and tested, both with overall good results. A local adaptive time step algorithm based on the frozen flux method by Zhang et al. (1994) was implemented within the CA model. The incorporation of the algorithm reduced the computation time in all the simulations, with no loss of accuracy in the results. The advantage against the version incorporating a global adaptive time step was variable, according to the flow dynamic and the model version: the local time step technique was more effective with the diffusive version of the model, since in twodimensional cases the speed up went from 2.5x to 4x. The reduction of run time was less relevant when the inertial formulation was used, although it should be noted that the numerical cases included only regular grids with a limited number of cells. Moreover, an inertial formulation proposed by Bates et al. (2010) was implemented in place of the diffusive formulation originally used in the CA model, and tested; the new formulation was found to be very effective in all the presented cases, reducing computation times up to 97%. Such reductions of run times were comparable with those of the LISFLOOD-FP model. However, compared to the original diffusive formulation, the inertial formulation did not improve the accuracy of results. In particular, the analysis of two-dimensional cases suggested that the inclusion of the advection terms within the model equations will be necessary for a more accurate simulation of flow dynamics. In conclusion, the results suggest that the CA model can be a valuable tool for the simulation of flood inundation events, especially with the implementation of the local time step algorithm and the inertial formulation. Although only numerical cases have been presented and discussed, the CA model has already been applied to a wide variety of practical applications that will be presented in the future. At the present time, the model is under further development and will be used by the Civil Protection of Emilia-Romagna region for the simulation of inundation scenarios and to prepare flood management plans. Also, it must be noted that the CA model computation scheme is particularly suited for parallelization, which could increase dramatically the computational, allowing for model application over wide areas (Neal et al., 2009b). A further development could be the use of CA model for the reproduction of overland flow and associated processes, like soil erosion and sediment transport, both in plain and hill slope areas. In fact, the model structure is potentially able to manage wetting and drying effects over wide areas, although the effectiveness of the model in practical applications is currently under evaluation.

279

Acknowledgements The original version of the CA model was developed in collaboration with the Civil Protection Agency of Emilia Romagna. The authors wish to thank Timothy Fewtrell for his valuable suggestions and comments, and for providing additional helpful references. An anonymous reviewer is also acknowledged. References Bagarello, V., Ferro, V., 2006. Erosione e conservazione del suolo. McGraw-Hill, Milano (in Italian). Bateman, A., Medina, V., Velasco, V., 2010. Soil infiltration effect in flat areas floods simulation. In: Barcelona, J. Carrera, (Ed.), XVIII International Conference on Water Resources, CMWR 2010. Bates, P.D., Horritt, M.S., Fewtrell, T.J., 2010. A simple inertial formulation of the shallow water equations for efficient two-dimensional flood inundation modelling. J. Hydrol. 387, 33–45. Crossley, A.J., Wright, N.G., Chris, D., Whitlow, C.D., 2003. Local time stepping for modeling open channel flows. J. Hydraul. Eng. 129 (6), 455–462. D’Ambrosio, D., Di Gregorio, S., Iovine, G., 2003. Simulating debris flows through a hexagonal cellular automata model: SCIDDICA S3hex. Nat. Hazards Earth Syst. Sci. 3, 545–559. Di Giammarco, P., Todini, E., Lamberti, P., 1996. A conservative finite elements approach to overland flow: the control volume finite element formulation. J. Hydrol. 175 (1–4), 267–291. Directive 2007/60/EC on the assessment and management of flood risks, 2007. Official J. Eur. Union, L., 288, 27–34. Dottori, F., Martina, M.L.V., Todini, E., 2009. A dynamic rating curve approach to indirect discharge measurement. Hydrol. Earth Syst. Sci. 13, 847–863. Dottori, F., Todini, E., 2010. A 2D flood inundation model based on cellular atomata approach. In: Barcelona, J. Carrera (Ed.), XVIII International Conference on Water Resources CMWR 2010. FLO-2D Software Inc., 2007. FLO-2D Users Manual (2007). Fewtrell, T.J., Bates, P.D., Horritt, M.S., Trigg, M.A., 2007. The effect of temporal and spatial coarsening on storage cell predictions of urban flood inundation. In: 32nd Congress of the International Association for Hydraulic Engineering and Research, Venice, Italy, vol. 1, pp. 37–47, ISBN: 8889405066. Garcia-Navarro, P., Murillo, J., 2010. Finite volumes for the simulation of unsteady shallow water flows. In: Barcelona, J. Carrera (Ed.), XVIII International Conference on Water Resources CMWR 2010. Guy, B.T., Rudra, R.P., Dickinson, W.T., 1992. Process-oriented research on soil erosion and overland flow. In: Parsons, J., Abrahams, A.D., (Eds.), Overland Flow, pp. 225–242. Horritt, M.S., 2002. Evaluating wetting and drying algorithms for finite element models of shallow water flow. Int. J. Numer. Meth. Eng. 55, 835–851. Horritt, M.S., Bates, P.D., 2002. Evaluation of 1-D and 2-D numerical models for predicting river flood inundation. J. Hydrol. 268 (1–4), 87–99. Horritt, M.S., Di Baldassarre, G., Bates, P.D., Brath, A., 2007. Comparing the performance of a 2-D finite element and a 2-D finite volume model of floodplain inundation using airborne SAR imagery. Hydrol. Process. 21, 2745–2759. Horton, R.E., 1945. Erosion development in streams and their drainage basins; hydrophysical approach to quantitative geomorphology. Bull Geol. Soc. Am. 3 (3), 340–349. Hunter, N.M., Horritt, M.S., Bates, P.D., Wilson, M.D., Werner, M.G.F., 2005. An adaptive time step solution for raster-based storage cell modelling of floodplain inundation. Adv. Water Res. 28 (9), 975–991. Hunter, N.M., Bates, P.D., Horritt, M.S., Wilson, M.D., 2007. Simple spatiallydistributed models for predicting flood inundation: a review. Geomorphology 90, 208–225. Hunter, N.M., Bates, P.D., Neelz, S., Pender, G., Villanueva, I., Wright, N.G., Liang, D., Falconer, R.A., Lin, B., Waller, S., Crossley, A.J., Mason, D.C., 2008. Benchmarking 2D hydraulic models for urban flooding. Water Manage. 161, 13–30. Hydraulic Engineering Centre (HEC), 2001. HEC RAS Hydraulic Reference Manual, US Army Corps of Engineers, Davis, California, USA, 377 pp. Jia, Y.F., Wang, S.S.Y., 2001. CCHE2D: Two-dimensional Hydrodynamic and Sediment Transport Model for Unsteady Open Channel Flow Over Loose Bed. NCCHE Technical Report. NCCHE-TR-2001-01. Kleb, W.L., Batina, J.T., 1992. Temporal adaptive Euler/Navier–Stokes algorithm involving unstructured dynamic meshes. AIAA J. 30 (8), 1980–1985. Müller, S., Stiriba, Y., 2007. Fully adaptive multiscale schemes for conservation laws employing locally varying time stepping. J. Sci. Comput. 30 (3), 493–531. Neal, J.C., Bates, P.D., Fewtrell, T.J., Hunter, N.M., Wilson, M.D., Horritt, M.S., 2009a. Distributed whole city water level measurements from the Carlisle 2005 urban flood event and comparison with hydraulic model simulations. J. Hydrol. 368, 42–55. Neal, J.C., Fewtrell, T.J., Trigg, M.A., 2009b. Parallelisation of storage flood models using OpenMP. Environ. Modell. Softw. 24 (7), 872–877. Osher, S., Sanders, R., 1983. Numerical approximations to nonlinear conservation laws with local varying time and space grids. Math. Comput. 41 (164), 321–336. Pender, G., Néelz, S., 2010. Benchmarking of 2D hydraulic modelling packages. SC080035/R2 Environment Agency Bristol, 169 pp. Sanders, B.F., 2008. Integration of a shallow water model with a local time step. J. Hydraul. Res. 46 (4), 466–475.

280

F. Dottori, E. Todini / Physics and Chemistry of the Earth 36 (2011) 266–280

Tayfur, G., Levent Kavvas, M., Govindaraju, R.S., 1993. Applicability of St. Venant equations for two-dimensional overland flows over rough infiltrating surfaces. J. Hydraul. Eng. 1, 51–63. Xanthopoulos, T., Koutitas, C., 1976. Numerical simulation of a two-dimensional flood propagation due to dam-failure. J. Hydraul. Res. 14 (4), 321–331.

Yu, D., Lane, S.N., 2006. Urban fluvial flood modelling using a two dimensional diffusion-wave treatment. Hydrol. Process. 20 (7), 1541–1565. Zhang, X.D., Trepanier, J.Y., Reggio, M., Camarero, R., 1994. Time-accurate local time stepping method based on flux updating. AIAA J. 32 (9), 1926–1929.