Detailed Numerical Modeling of Flood Flow in Floodplains with Complex Geometry

Numerical modcling of Jlood flou and the eaahtation of flood hautrds can be based on aarious numerical modzls and, rnodeling technigues. One-dirnensional ( 1 D), quasi two-dimen-sional ( I ,5 D), tuo-dimensional (2D) or three-ditnensional (iD) aaria,nts of nurneri.cal modcls can be used. While 3D modek are too demand:ing to be used for flood flou modeling on personal cornputers, quasi-2D and 2D modek can be more uidely used to solae eaen larger practical problems noua.days. Detailed two-dimensional numerical rnodeling of flood flow in flooded urbanized areas with cornplex geometry using the 2D dlpth aaeraged model is presmted in this paper. The gouerning equations of the model are expressed, with a set of depth aueraged Reynol.ds equations consisting of the continuity eEta,tion and two momentum equations for the horizontal aelocity components. The eddy uiscosity, which influences the horizontal turbulent momentun?, ercchange processes, 'is modeled with the aid of a fupth aaerage aersion of the two-equation h-e turbulence model. The partial dffirentinl equations are solved numeritallg with a control uolume method using f.ne non-orthogonal curuilinear grid"s and a non-staggered aaria,bk artangernent. The applicability and adaantages of this modcling approach for simulating, the Jlood flow in lloodplains with complex geometry, and in urban areas, are illustratedby the results of apilot study in Choceit and Usttnol Oilict, towns whichwere severely affected during the 1997 and 1998 floods in the Czech Republic.


I Introduction
Experience resulting from the extreme flood events in the Czech Republic in recent years have stressed the necessity of further development and practical evaluation of suitable methodologies enabling more detailed investigation of flood flow characteristics in domains with very complex geometry and in urbanised areas.Various numerical rnodels and mod- eling techniques [] can be used and are being developed and improved according to the growing performance of computers.The classical approach to the numerical modeling of flood events is based on the application ofvarious types of lD models.These models enable the longitudinal water elevation distribution to be estimated along the river axis for the flood peak discharge (steady case calculation), or analyses to be made of flood wave propagation (unsteady simulation).
Flooding lines were, and in most cases still are, estimated and offrcially delimited in the Czech Republic with the use of lD modeling technology.
Current trends in the field of numerical modeling of flood events are oriented towards the development and practical utility of multi-dimensional models.While 3D models are still too demanding to be used for flood flow modeling, quasi-2D and 2D models can nowadays be more widely used to solve larger practical problems, even on personal computers.These complex models can provide results of substantially higher quality than those of lD modeling tools.Togetherwith the water elevation distribution, 2D models provide a lot of other detailed information not available in the classical lD numerical model output (flow pattern evaluation, distribu- tion of velocity magnitudes and directions (velocity fields), streamlines, etc.), which is extremely important information for hydraulic and other follow-up (flood risk, flood damage) analyses of flood events.
The Faculry of Civil Engineering at the Czech Technical University is a co-bearer of the Czech-Swiss FI-AMOR (Flood   Analyses and Mitigation on the Orlice River) project, which was initiated by the flood events in the Czech Republic in 1997.The project is supported by Swiss Humanitarian Aid and Disaster Relief (SFIA+DR).Within the framework of this project a detailed case study dealingwith hydraulic and flood risk analyses in the urban areas of Chocei and Usti nad Orlici was elaborated.The modeling results for Chocei from this study are used as examples in this paper.

Numerical model
The applied numerical model is based on the FAST 2D model originally developed in Germany at the Institute for Hydrodynamics in Karlsruhe [2].The model enables the simulation of free surface steady water flow in domains with complex geometry Further development of the model, aimed at the design and programming of suitable preand post-processing tools, has been carried out by Hydroexpert Ltd., Prague, in cooperation with the Bundesanstalt ftir Wasserbau in Karlsruhe [3].

Goaeming eqntions
The model is based on the set of depth averaged Reynolds equations often called'shallowwater equations', which can be derived by integration of the Reynolds equations for three dimensional flow over the depth of the water layer [4].The resulting set of governing equations consists of the continuiry equation and two momentum equations for the horizontal velocity components.It can be written in the following forrn: 0u; 0u; A ,, -:J+ u i-= -t -lh+26)+ At tAx, "dr' 1: .**(o4r)" s, ; i, j =r,2 (r) on ox; .J dh.dl hu;l :---r _-\--i-lz =n At 0*.; where I is the time [s], A is the water depth [m], z, denotes the vertical coordinate of the bottom level [m], p is the rvater density [kg.--'], g is the acceleration due to gravity [*.r-'].
The velocities ulandu2tm's-'] are the depth averaged veloc- ities in thexl andx2 directions, the source term S; [rn.s-'] in the momentum equations includes the influence of external forces and stresses, such as bed shear, stress due to wind on the water level and the Coriolis acceleration.The depth-aver- aged components of the effective stress tensor are defined as the sum of laminar stresses r/, turbulent stresses t, and stresses t2 resulting from the non-uniform velocity distribution in the vertical direction: The individual components are defined as Solution ttrocedrne 'fhe partial differential equations are solved numerical- ly with a finite volume procedure based on a method for pure two-dimensional flow, which is described in detail in Majumdar [6] and Rodi et al [7].The method employs non--orthogonal curvilinear grids and a non-staggered variable arrangement, i.e., all variables are stored at the center of the control volumes.This technique requires the use of a special momentum interpolation scheme in order to avoid unnatural pressure oscillations (for details, see Majumdar [6]).The scherne further uses Cartesian velocity components and hybrid central/uprvind diflerencing for the discretization of convective tenns in the momentum equations.The pressure-velocity coupling is achieved with the SIMPLEC iteration algorithm, in which the momentum equations are solved first with the guessed pressure field, and an improved pressure distribution is then calculated via a pressure correc- tion equation.The resulting system of linear equations at each iteration is solved with the Thomas tridiagonal matrix algorithm (TDMA).The implicit coupling between the pressure and water elevation (or water depth) in the depth average model is solved in an outer iteration loop.

Input data
The main input data for the model includes: o geometrical data describing the computational grid, o geometrical data describing the terrain topography, o a definition of the fixed non-overflorved obstacles (build-  ings, dikes, roads, etc.), r distribution of bottom shear stress parameters (Manning roughness coe{ficients), o boundary conditions.
3 Application to case study Stud.yarea Chocei lies on the Tichii Orlice river at 290 m.a.s.l.It is an industrial center with textile and machinery industry and an irnportant railway junction.The town has a population of 9 200 inhabitants (in 1530 houses).The urbanised area (residential and public buildings) covers a surface area of 2.1 square km.Fig.I shows an aerial photograph of the flooded town area during the flood in 1997.

Input dala and, model construction
The model created for Chocei covers the flooded area about 3.8 km in length and I km in width.The left and right side boundaries were chosen in such a way that the whole width of the flood plain inundated during the extreme flow events fits into the model area.
'Ihe curvilinear computational grid was generated using the model boundaries as input data.The number of control volumes in the longitudinal direction was 1470, while in the transverse direction 500 control volumes were used.The total number of computational celis was 735 000.The large number of control volumes enabled sufliciently detailed modeling ofstructures in the urban areas (the average size ofa control volume was about 2.5x2 m).In the next step in constructing the model, all the control volumes iocated inside buildings z6+h 1 | , , .)' ,l / z\/ z\ ri =-p(ui -ui )\uj -uj), where v in the laminar stress term is the kinematic viscosity [t rt.r-t].Turbulent Reynolds stresses are expressed using the turbulent fluctuations of the velocity components zli, and the bar denotes an averaging operation.Due to these terms, the system of governing equations is not closed and has to be supplemented with some suitable model of turbulence.The effective stresses T;.i are modeled by the eddy viscosity ap- proach in the model: ( ^". i,, ) rii =,,p1?.? l, " \oflj oxi ) where v, is the eddy viscosity [-2.r-'].The eddy viscosity which models the horizontal turbulent momentum exchange processes is calculated with the aid ofa depth average version of the A-e turbulence model developed by Rastogi and Rodi   [51.In this model the eddy viscosity is related to the turbulent kinetic energy fr [m'-s-l of the turbulent motion and to its dissipation *t. e 1*2.r-31 by -9 vt =cp e. (5) The horizontal distribution of A and t over the modeled domain is calculated by solving partial differential transport equations for these wo variables.Due to this modeling ap- proach the final set of governing equations solved in the model consists of a total of five partial differential equationstlvo momentum equations with a continuity equation and two transport equations for the turbulence parameters.

#i1
Fig. 2: Detail of the computational grid and other obstacles were blocked out, and the boundalies of the blocked sub-domains were then treated as vertical walls during the simulations.Frg. 2 shorvs one detail of a small com- putational grid part of the Chocei 2D model.The crossed rectangles represent the blocked out computational celis.
The terrain approximation in the final numerical model was created rvith help of a digital terrain model clcated on the basis of an evaluation of the aerial snaps.The raw digital model data had to be complemented with the topography of the river channel and some other geometrical entities rrefining the terrain reptcsentation in the model.The aerial photos were taken in the sumtner of 1998 by fugus Geosyst6m Ltd.The evaluation of the photos was done by Topograf Ltd., Prague, in order to create a Digital Terrain Model (DTM).
The raw data includes the elevations of all breaklines and also terrain elevations at a spacing of l0 metres.The declared data accuracy is 0.5 m in the horizontal dircction and 0.25 m in the vertical dircction The aerial photo DTM data n'as complemented with sur- vey information: eightl' cross-sec[iorls at intetlals of 50 to 100 m.The cross-sections at hydraulic structulcs as rvell as the geolnetln o[ the sturctuucs themselves \vel'e specificalh' sur- veyed.Missing topological information in small parts of the terrain outside the evaluated aelial snaps also had to be added to the terrain data, using maps on a scale of l:5 000.
All these changes and tcfinements were perfonned inter- actively, rvith the aid of the model Preprocessor.Finally, the terrain data rvas interpolated to the nodes of the com- putational grid.The final assembled model is graphically documented in Fig. 3, which contains a 3D vierv of the model geometry (with the geometry of the buildings and roads sym- bolically depicted).The boundary condition prescribing the total inflow into the modelwas specified at the inlet boundaries alongdith the prescribed flow direction.Non-uniform distribution of the discharge along the inlet boundary was generated automatically inside the model algorithm according to the actual water depth at each boundary node.The fully developed turbulent flow with equilibrium between the production of turbulent kinetic energy and its dissipation was used as a boundary condition for the turbulence parameters.At the downstream boundary the water elevation was prescribed as a boundary condition, For velocities and turbulent parameters a Neu- mann type boundary condition was used, assuming null derivatives of variables in the normal direction.

Modcl calihrotinn
The model was calibrated using the available calibration data.During the calibration a discharge equal to the peak discharge of the 1997 flood was simulated, which is almost equal to the 1OO-year flood discharge.Sub-domains with different characteristics of hydraulic roughness were identified using aerial photographs and other available information.
Frfteen different zone types were used.The corresponding computational cells for each zone were interactively identified in the computational grid and were given appropriate Manning coeflicient values.The values of the roughness coel ficients were adjusted by trial and error during the calibration simulation runs, with the aim of getting the best approximation of the observed water levels.The spatial distribution of $9"-. the calibration points is shown in l1g.4, together with a graphical evaluation of the final calibration results.The calibration results are expressed as the differences betlveen the calculated and observed water levels.In the same figure, there is a cornparison of the observed (dashed line) and computed (solid line) extent of the flood.

Evaluation and utilization of results
In order to make practical application of the numeri- cal model easier and more eflicient, a system of preand post-processing tools fbr the FAST 2D model has been developed.The system, called PREFAST [3], is oriented towards the use of personal computers, and was programmed as an application based on the ADS (AutoCAD development system) for AutoCAD graphical software.The system inherits interactive user-friendly tools with a graphical interface for all steps of the model designgrid generation and modification, creating, editing and exploiting the digital terrain model, specifying obstacles, specifying the distribution of bottom roughness coeffrcientsand for a graphical evaluation ofthe numerical simulation results.
The complete set of standard graphical interpretations of the simulation results includes: o streamlines.This standard evaluation set (flood hazard maps) plovides detailed information about the spatial distribution of flood flow characteristics, and enables the analyses of flood flow conditions in the whole flooded area.The results enable complex flood evaluation, including for example the identification of main flow paths, specification of active and passive flood flow zones, evaluation of flood intensities, identillcation ofweak points in the system, and recognition of main obsta- cles influencing the flood flow.Moreove4 the primary results can be further used as input data for more complex flood analyzing methods, i.e., flood risk analyses and flood damage predictions.
Together with standard flood hazard mapping, some other special graphical representations of the results can also be applied.These evaluations provide a more detailed analysis of the results and include for example a local evaluation of the velocity fields (see example in Fig. 6), or they can be used for visual presentation to the public (see the 3D presentation of the extent of the flood in Chocei in lig.7, as an example).
The results presented in this study document the current potential of very detailed and non-simplified numerical models of flood flow in domains with very complex geometry and in urban areas.Analog studies can be directly applied to practical problems encountered in water management and land-use planning, e.g., in formulation of guidelines for land use in flooded areas, categorization of flood planes, preparation of municipal plans, evaluation of impacts of human activities in flood plains on flood flow conditions, elaboration and evaluation offlood control projects, etc.

Fig. 3 :
Fig. 3: A 3D view of the numerical model