Optimisation of LiDAR derived terrain models for river flow modelling

Airborne LiDAR (Light Detection And Ranging) combines cost efficiency, high degree of automation, high point density of typically 1–10 points per m 2 and height accuracy of better than ±15 cm. For all these reasons LiDAR is particularly suitable for deriving precise Digital Terrain Models (DTM) as geometric basis for hydrodynamicnumerical (HN) simulations. The application of LiDAR for river flow modelling requires a series of preprocessing steps. Terrain points have to be filtered and merged with river bed data, e.g. from echo sounding. Then, a smooth Digital Terrain Model of the Watercourse (DTM-W) needs to be derived, preferably considering the random measurement error during surface interpolation. In a subsequent step, a hydraulic computation mesh has to be constructed. Hydraulic simulation software is often restricted to a limited number of nodes and elements, thus, data reduction and data conditioning of the high resolution LiDAR DTM-W becomes necessary. We will present a DTM thinning approach based on adaptive TIN refinement which allows a very effective compression of the point data (more than 95% in flood plains and up to 90% in steep areas) while preserving the most relevant topographic features (height tolerance ±20 cm). Traditional hydraulic mesh generators focus primarily on physical aspects of the computation grid like aspect ratio, expansion ratio and angle criterion. They often neglect the detailed shape of the topography as provided by LiDAR data. In contrast, our approach considers both the high geometric resolution of Correspondence to: G. Mandlburger (gm@ipf.tuwien.ac.at) the LiDAR data and additional mesh quality parameters. It will be shown that the modelling results (flood extents, flow velocities, etc.) can vary remarkably by the availability of surface details. Thus, the inclusion of such geometric details in the hydraulic computation meshes is gaining importance in river flow modelling.


Introduction
Hydrodynamic-numerical (HN) modelling -also referred to as Computational Fluid Dynamic (CFD) modelling -is widespread in water sciences.Simulation tools are used for ecological studies (Hauer et al., 2008a,b) as well as for hydraulic engineering purposes such as river restoration and the like (Denoth, 2003;Gouda et al., 2002;Khan and Barkdoll, 2000;Fischer-Antze et al., 2001).Moreover, the precise delineation of endangered or vulnerable areas in flood hazard analysis is also based on hydrodynamic-numerical models (Habersack and Gaul, 2008).In many countries, land use regulations and building codes are directly linked to flood hazard maps, which show areas affected by events with a certain return period.
In the UK, France, United States, Canada, and New Zealand, the area affected by a 100-year flood (HQ100) plays an essential role for flood mitigation (Marco, 1994;Watt, 2000).In Spain spatial planning of flood areas is included in the Water Act and institutional regulations (Menendez, 2000).Switzerland heavily engages in mapping activities to identify zones that are prone to natural hazards (e.g., Petraschek, 2002).In Norway maps are created for six different G. Mandlburger et al.: Optimisation of LiDAR DTMs for river flow modelling flood levels (HQ10/20/50/100/200/500) following uniform technical guidelines (Hoydal et al., 2002).In Austria hazard zone maps are developed at a scale of 1 : 2000, in which three zones (red, yellow and blue) are distinguished (RIWA-T, 2006).Red zones are determined by a velocity-depth criterion, whereas yellow zones are bounded by the extent of a hundred-years flood.Blue zones exhibit areas of specific flood management.Additionally, rest risk areas (HQ300) and the spatial extent of the thirty year event (HQ30) are presented in Austrian flood hazard mapping (RIWA-T, 2006).Given the importance of event simulations for planning, their quality is crucial.
The most influential input parameters for HN models are the topography provided as a Digital Terrain Model of the Watercourse (DTM-W) and the flow resistances usually provided as roughness coefficients.According to technical advances in computer aided hydrodynamic-numerical modelling and the related multidimensional (2-D/3-D) applications, the need for high quality 3-D terrain data is increasing (Sinha et al., 1998).Airborne LiDAR (Light Detection and Ranging) provides high point density and remarkable height precision (Csanyi and Toth, 2007).It is especially well suited for capturing flood plains as well as river banks under low flow conditions and is therefore widely used as data basis for HN models.
However, the application of LiDAR poses problems as well.Contrary to traditional manual data acquisition techniques like stereo-photogrammetry or tachymetry, LiDAR data includes both terrain and off-terrain points on buildings, vegetation, power lines, etc.The quality of the derived DTM and, consequently, the quality of subsequent HN modelling depends crucially on how well off-terrain points have been eliminated within the filtering process.Another issue besides filtering is the fusion of LiDAR and river bed data.This involves the derivation of the water-land-boundary and the interpolation of river bed cross sections.
Nowadays, most HN models are solved using a Finite Element (FE) or Finite Volume (FV) approach on the basis of unstructured or hybrid geometries, i.e. a computation grid based on irregularly distributed points.Many hydraulic mesh generators are available, which build up a net of polygonal surface elements covering the entire project area.Mainly triangles and arbitrary quadrilaterals are in use.Hereby, the vertices of the quadrilaterals don't need to be co-planar.These mesh generators consider hydraulic parameters like angle criterion and aspect ratio, but normally disregard topographic details as contained in data from modern LiDAR sensors.The heights are mapped to the hydraulic grid a posteriori.Many modelling tools can only handle a limited number of points and elements (Nujic, 1999).Thus, the direct use of the high resolution LiDAR DTM as computation grid is impossible due to the enormous amount of data.
This article presents a new method for constructing hydraulic computation grids by preserving geometric surface details as well as considering physical mesh quality param- eters.An adaptive TIN refinement method is used to thin out the high resolution LiDAR DTM, which is professionally conditioned in a subsequent postprocessing step to meet the required mesh quality parameters.The method, thus, overcomes the limitations of mesh generation approaches considering physical parameters only.
In the following Sect. 2 the study area and data sets are introduced.Thereafter, the major steps of the proposed workflow are described in detail, comprising LiDAR DTM processing (Sect.3), DTM data reduction and conditioning (Sect.4) and hydrodynamic-numerical modelling (Sect.5).The results and application examples are presented and discussed in the subsequent Sects.6 and 7.The final section concludes the major findings and points out the benefits of our approach.

Study area and data sets
Special meteorological conditions in Lower Austria and the south of the Czech Republic were responsible for an extreme flood in summer 2002 at the river Kamp (Gutknecht et al., 2002).A precipitation event of 150-370 mm was documented in the first 12 days of August 2002, 3-4 times higher than the average monthly precipitation in the catchment area of the river Kamp (Holzmann, 2002).These massive rainfalls caused a flood with two distinct peaks.The first flood wave occurred on 8 August (c.f.Fig. 1).Within 24 h the discharge increased from 100 m 3 /s to 800 m 3 /s (gauging station Stiefern, HQ100 at Stiefern: 490 m 3 /s).The recorded discharges were nearly twice as high as the 100-year event for gauging station Stiefern (estimated return period 500-2000 years).On 14 August the second wave occurred with a maximum discharge of about 500 m 3 /s.During the flood The study reach is a 2 km section midstream of the river Kamp situated in the northern part of Lower Austria at the village Gars am Kamp (Fig. 2).The origin of the river Kamp is near the village Karlstift and after a 160 km long course it discharges into the Danube River (182 m a.m.s.l.) at Altenwörth.The total catchment area measures 1753 km 2 and the hydrologic regime is described as pluvio-nival (Mader et al., 1996) with an average annual precipitation of 900 mm.The average bed slope of the selected study reach is approx.2.5‰.The investigated river section of Gars am Kamp was heavily affected by the catastrophic flooding in August 2002.A high percentage of buildings (20.1%, 293 of 1460) were damaged or even destroyed.As entire losses a sum of more than 13 million C was reported by the Lower Austrian Federal Government (Habersack et al., 2004).
As a result of the extreme 2002 flood event a LiDAR campaign was carried out for the entire Kamp valley.The data were captured during 25 November 2002 and 2 March 2003 by the company TopoSys (TopoSys, 2009) using a Falcon II sensor at an average flying altitude of 850 m above ground.The reported accuracy of the captured 3-D-points was ±0.5 m in planimetry and better than ±0.15 m in height.The Falcon II scanner measures at a pulse repetition rate of 83 kHz.The scanner is mounted on an oscillating platform, which produces a meandering scan pattern on the ground to improve the inhomogeneous point distribution of the scanner (more points in flight direction than in scan direction).First echoes and last echoes together with their signal amplitude were recorded at an average point density of 1 point/m 2 .A Digital Surface Model (DSM, c.f. Fig. 2c) was derived from the original point recordings consisting of approx.1.16 million grid points for the entire study reach, and 700 000 points were classified as terrain points.Additionally to the LiDAR points, cross sections at an average distance of 60 m were captured by terrestrial survey to describe the topography of the river bed.

LiDAR DTM processing
In this section the state of the art in DTM acquisition in general and obtaining terrain models from airborne LiDAR in particular will be reviewed.
For the specification of terrain surface products the Digital Terrain Elevation Data (DTED) levels from 1 to 5 have been suggested.For HN simulations especially the resolution and accuracy levels are of interest.Level 2 is fulfilled by DTMs with a horizontal resolution (grid edge length) of 1 (approx.20-30 m) and an absolute vertical accuracy of 18 m (90% error bound).It may be reached by spaceborne Interferometric Synthetic Aperture Radar (InSAR) using the X-band (e.g.SRTM, Rabus et al., 2003), which is, however, scattered back at the tree crowns.Level 3, also referred to as High Resolution Terrain Information (HRTI-3), with 12 m resolution and 10 m accuracy, is aimed at by new satellite InSAR missions (Krieger et al., 2007).HRTI-4 (6 m resolution, 5 m accuracy) can be achieved by airborne InSAR and stereo photogrammetry (Kraus, 2007).At this level it makes a notable difference if the ground surface is recorded or the first visible surface from the bird's perspective (tree crowns, roofs, etc.).Short wavelength radar and passive optical imaging can provide high resolution, but cannot take measurements from the ground surface below the forest cover (Baltsavias et al., 2008).HRTI-5 can currently only be provided by airborne LiDAR for larger areas in a cost efficient manner and prescribes a resolution of 1 m, an absolute vertical accuracy of 5 m, and a relative (point-to-point) precision of 25 cm.Especially the accuracy requirements are the fact advocating for airborne LiDAR.
In airborne LiDAR the travel time of a short pulse of laser energy from the sensor to the ground surface and back is measured.The position and orientation of the sensor platform is determined with GNSS and IMU (Global Navigation Satellite System, Inertial Measurement Unit).The laser pulse can "see" through small gaps in the foliage and travel to the ground below the tree canopy.The result is a cloud of points originating from the ground surface and the canopy in a forest environment.In open areas the points are only lying on the ground surface.The task of separating the points on the ground surface from the other, off-terrain points is called filtering (Kraus and Pfeifer, 1998).Different approaches have been suggested (Pfeifer and Mandlburger, 2008) making different assumptions on the terrain surface and either proceed by gradually adding or by gradually removing unclassified points.Finally, advanced algorithms also consider the random measurement noise.
The methods of mathematical morphology (Kilian et al., 1996;Vosselman, 2000;Zhang et al., 2003) classify points depending on vertical and horizontal distances between points.The principle is that a large height difference at small horizontal distance can only occur, if the higher point is off the ground.In progressive densification (von Hansen and Vögtle, 1999;Axelsson, 2000;Sohn and Dowman, 2002) the terrain surface is described by a triangulation, that is gradually built up.Starting from a few known terrain points, new points are added, if they are within certain distances from the triangulation surface.With these algorithms also a notion of surface is introduced.In robust interpolation (Kraus and Pfeifer, 1998) a surface model is interpolated starting from all points, more precisely approximating all points, e.g. by linear prediction or kriging (Kraus, 2000;Journel and Huijbregts, 1978).This surface model runs in an averaging way between terrain and off-terrain points.The residuals, i.e. the vertical distances from the measurements to the surface model, are computed and become input for a weight function.This way, high weights are associated with points below the current surface, and small weights with points above it.These weights are considered in the next iteration step of surface computation.Points with larger weights attract the surface, while points with lower weights have less and less influence in subsequent iterations.The approach has been embedded in a hierarchic setup, for increasing robustness and efficiency (Pfeifer et al., 2001).A comparable method has been suggested by Elmqvist et al. (2001).The advantage of this approach is that the surface properties, expressed e.g. by the variogram of kriging, are clearly separated from the properties of LiDAR data (terrain and off-terrain points), which are expressed in the weight function.Random measurement errors are considered, if a suitable surface interpolation technique like linear prediction is chosen.Finally, segmentation approaches were suggested (Sithole and Vosselman, 2005;Nardinocchi et al., 2003;Tóvari and Pfeifer, 2005) for separating terrain and off-terrain points, which first group the points to larger entities (e.g.planar areas) and then classify those.They are mainly advantageous in city areas, where a large number of well-defined, planar surfaces do exist.
As reported in a comparison of filter algorithms (Sithole and Vosselman, 2004) all filters have their problems at step edges.As shown by Doneus and Briese (2006) also low vegetation can be very problematic for correctly identifying the terrain points automatically.With the advent of commercial full-waveform LiDAR systems (Wagner et al., 2004) a step forward in filtering became possible.In full-waveform Li-DAR not only the travel time of the pulse is recorded, but the entire shape of the backscattered echo is digitised.In open areas the ground surface is well-defined and the emitted shot of laser energy is scattered back undistorted.For vegetation, however, the situation is different.The "surface" of bushes or tree crowns is vertically much more extended with respect to the laser footprint (50 cm-1 m), backscattered echoes are therefore widened.This can be used as a pre-classifier to detect points lying on the vegetation, and thus, ground detection becomes more robust and accurate (Doneus and Briese, 2006).

DTM data reduction
Modern LiDAR sensors allow mapping of topographic details at the price of a highly increased data volume.In order to obtain a manageable amount of data for FE meshes used in HN modelling data reduction and surface simplification is inevitable.The main goal is to reduce the vertices of the DTM (grid points, breaklines, etc.) without losing relevant geometric details.
For polygonal simplification of 2.5-D surfaces mainly regular grid algorithms, decimation and refinement techniques are in use (Heckbert and Garland, 1997).The regular grid approaches are widespread, simple and fast, but they are not adaptive and produce poor approximation results.Better approximation quality can be achieved by applying decimation and refinement methods based on general triangulation algorithms such as Delaunay triangulation.Decimation methods work from fine-to-coarse and are not suited for processing large high resolution LiDAR DTMs since they require a triangulation of the entire point set.Refinement methods represent a coarse-to-fine approach starting with a minimal initial approximation.In each subsequent pass one or more points are added as vertices to the triangulation until the desired approximation tolerance is met or the desired number of vertices is used.
The performance of DTM data reduction is highly influenced by the existence of systematic and random measurement errors.For deriving thinned DTMs from original Li-DAR point clouds, refinement or decimation approaches as described above are not the first choice.Systematic errors have to be removed first by rigorous sensor calibration and fine adjustment of the LiDAR strip data (Schenk, 2001;Filin, 2003;Kager, 2004).The random measurement errors of the LiDAR points should rather be eliminated by applying a DTM interpolation strategy with measurement noise filtering capabilities such as linear prediction (Kraus, 2000) or kriging (Journel and Huijbregts, 1978).High reduction ratios can only be achieved for DTMs, where systematic errors have been strictly eliminated and random errors have been reduced during surface interpolation.
In contrast to approaches which rely on the original point cloud, our adaptive TIN refinement approach uses the filtered hybrid DTM (regular grid and breaklines, structure lines, and spot heights) as input (Mandlburger, 2006).The basic parameters for the DTM simplification are a maximum height error z max (L ∞ -norm) and a maximum planimetric point distance xy max .The latter avoids triangles with too long edges and narrow angles.The algorithm starts with an initial approximation of the DTM comprising all structure lines and a coarse regular grid (cellsize = xy max = 0 ), which are triangulated using a Constrained Delaunay Triangulation.Each 0 -cell is subsequently refined by iteratively insert-ing additional grid points until the height tolerance z max is reached.
The additional vertices can either be inserted hierarchically or irregularly.In case of a hierarchical breakdown, a coarse rectangular grid is used as initial surface approximation.Each coarse grid cell is divided into four parts in each pass, if a single grid point within the regarded area exceeds the maximum tolerance.This leads to a quad-tree like point distribution in the resulting surface approximation.By contrast, for irregular division a regular grid of (approximately) equilateral triangles is used as initial approximation.This is especially suited for FE-meshes since the connecting line of adjacent cell centers (i.e.direction of calculation) always intersects the joint cell boundary at a right angle.In a single pass the grid points with the maximum positive and maximum negative deviation are inserted (parallel greedy insertion).Furthermore, the process of point insertion and recalculation of the error measures is optimised by considering the locality of the mesh changes.Higher compression ratios (up to 99% in flat areas) can be achieved with irregular point insertion, whereas the hierarchical mode is characterised by a more homogeneous vertex distribution.The hierarchic division, thus, produces an adaptive cartesian grid whereas irregular division yields an unstructured grid.

Data conditioning
To get a high quality computation grid for HN modelling, the simplified polygonal surface mesh has to be conditioned considering the physical phenomena being simulated (flow direction, bottom shear stress, . . .).According to Ferziger and Peric (2002) the main quality parameters for a good FE grid are: (i) angle criterion, (ii) aspect ratio and (iii) expansion ratio.The aspect ratio describes the proportion of the longest and shortest edge within a mesh face, whereas the expansion ratio is a measure describing the area ratio between adjacent faces.Requirements concerning angles in a computational mesh appear in two different ways.First, the angles between the mesh edges and, second, the alignment of the cells with respect to the direction of force have to be considered.Ferziger and Peric (2002) report that small angles of less than 10 • should be avoided, the cells should be aligned to the principal flow direction, the aspect ratio should not exceed 10 (optimum<3) and the expansion ratio must not be greater than 3 (optimum<1.2).
Thus, to achieve physically reliable results the cells of the computation mesh have to be aligned alongside the current for the entire wetted perimeter of a river under high flow conditions.This applies to the river bed and the river bank.The river axis can serve as a sufficient approximation of the principle flow direction.Quadrilateral cells (forming a hyperbolic paraboloid) with the longer sides in and the shorter ones perpendicular to the flow direction with edge lengths proportional 3:1 (optimum aspect ratio) have turned out to produce good simulations results (Mandlburger, 2006).Beyond the embankment, the flow directions are no longer strictly parallel, thus, irregular data distribution as described in the previous subsection is appropriate.
However, especially for flood mapping different zones can be designated within the inundation area depending on the respective vulnerability.Residential areas in the vicinity of a river, for instance, should be mapped in more detail than remote or elevated areas.The concept of varying mapping precisions can be applied to the DTM simplification approach by introducing spatially varying maximum height tolerances.Mathematically the tolerances can be expressed as a bivariate function ( z max = f (x, y)).Distances from the river bank or relative height differences, respectively, can be used to control the tolerances, but a simple zonal model has turned out to be suited best.The delineation of the different accuracy zones can be derived from a pilot survey (e.g. a preliminary 1-D-CFD simulation), from DTM visualisations such as hill shadings, from existing maps or the like.
Finally, the resulting TIN has to be analysed with respect to the adherence of angle criterion, aspect and expansion ratio.Therefore, the entire triangulation is traversed and, if necessary, additional vertices are inserted or hydraulically unessential vertices are removed.Figure 3 shows the original high resolution LiDAR DTM-W whereas a simplified grid, professionally conditioned to preserve geometric details and fulfill physical quality criterions, is illustrated in Fig. 4.

Hydrodynamic-numerical modelling
To evaluate the impact concerning the availability of geometric details on the results of river flow simulations, the two-dimensional depth-averaged hydrodynamic-numerical model Hydro AS-2-D (Nujic, 1999) was applied in the study reach.Basic principles of two-dimensional mathematical flow modelling are the 2-D depth-averaged equations (simplified Navier-Stokes-equations), also referred to as shallowwater equations.Within the Hydro AS-2-D model the shallow water equations are spatially discretised based on the finite volume approach.Simplified in vector form the shallow water equation can be described as follows: where: H is the water surface elevation [m], the sum of the water depth (h) and the terrain height (z).u and v are the velocities in x-and y-direction [m/s], g the acceleration due to gravity [m/s 2 ] and ν the kinematic viscosity term [m 2 /s].
Hydro AS-2-D uses the Surface water Modelling System (SMS) as a pre-and postprocessing tool.The convective flow of the two-dimensional model is based on the Upwindscheme by Pironneau (1989) and the discretisation of time is done by an explicit Runge Kutta method in second order.Several methods for implementing viscosity in numerical modelling exist.In the applied Hydro AS-2-D model the viscosity is calculated based on a combination of empirical and constant viscosity approaches.
where c µ is the viscosity parameter and v * the shear velocity [m/s] and ν 0 is a basic viscosity term.
In the applied Hydro AS-2D model ν 0 is used to ensure numerical stability.In general, different ν 0 -values can be assigned to each cell, but in practice ν 0 is only adapted in sections.For the simulations in this case study no ν 0 adaption was employed at all, a viscosity coefficient of c µ =0.6 was used, the computation time was set to 32 000 s to achieve steady state conditions and the depth criterion for wetting/drying computations was set to 0.01 m.The hundred-years flood (Gars: 470 m 3 /s) was used as input discharge (boundary condition) being the most important design discharge for hazard zone mapping in Austria (RIWA-T, 2006).Separate roughness n-values were determined for the main channel and for overbank areas.Whereas main channel roughness coefficients were calibrated on the basis of field measurements (flow velocities and water surface elevations of a HQ30), only a single value (n=0.04) was used for the study reach floodplains to best exhibit the influence of different geometries on flood level and flood hydraulics.This value characterises grassland and/or undisturbed surface terrains (Habersack, 1995).Normally, calibrated twodimensional hydrodynamic-numerical models for hundredyear floods exhibit very high roughness values (e.g.n=0.2) for settlements (ARGE Kamp, 2005), but in this study the loss in energy (secondary currents . . . ) is calculated based on the high quality DTM data and, therefore, does not need to  be considered in the summarizing roughness n-value.Anyway, due to the lack of reference data it is not the aim of the presented paper to achieve validated flood stages for the hundred-years event by calibrating different overbank parameters (e.g.vegetation, land use, . . . ) but rather to show the impact of geometry details provided by LiDAR on the modelling results.

Results
Before running the hydrodynamic-numerical simulation we had to build up a precise LiDAR based DTM-W covering the wetted perimeter of a HQ100 at least and including the river bed and all flow prohibiting buildings.Our data basis comprised the LiDAR points aligned in a 1 m raster, the river bed cross sections (terrestrial survey) and 2-D building polygons (cadastral map).In a first step, we derived a DTM of the overbank areas without gaps.Since the LiDAR terrain points at hand were sparsely distributed due to prior filtering, we merged all terrain and DSM points and analysed the point cloud with SCOP++ (SCOP++, 2009), which implements the hierarchic robust interpolation (c.f.Sect.3).Thereafter, we derived the water-land-boundary as described in Brockmann and Mandlburger (2001) serving as clipping polygon to cut out all laser echoes within the river bed.The bed surface itself was subsequently constructed by densifying the cross sections using an interpolation approach, which considers the curved progression of the river axis (Mandlburger, 2000).To reflect the natural flow conditions in our geometry model as realistically as possible, we decided to include the buildings as block models in the DTM-W.Therefore, we used 2-D building polygons from the cadastral map and estimated an average height for each building.The resulting 3-D polygons were introduced to the DTM-W in a final step, where we interpolated all preprocessed data (i.e.classified LiDAR terrain points and densified river bed points) using linear prediction (Kraus, 2000).Thus, random measurement errors in the height component in the order of ±15 cm were considered in the DTM computation.A grid with an edge length of 1 m containing additional breaklines was derived and the resulting hybrid DTM-W is illustrated in Fig. 3. Virtually all the tasks mentioned above were carried out using the SCOP++ software developed at the Institute of Photogrammetry and Remote Sensing, Vienna.
Additionally, data reduction and conditioning was performed as described in Sect. 4. We chose a maximum height tolerance z max =±20 cm and started the TIN refinement process with a basic triangular point pattern (edge length 0 =16 m).The irregular TIN of the river bed was automatically replaced by cell elements aligned to the flow direction.The resulting professionally thinned and conditioned DTM-W grid, which can be used as hydraulic computation mesh straightforwardly, is shown in Fig. 4.    To estimate the impact of additional topographic details provided by LiDAR on instream hydraulics, numerical simulations for HQ100 were carried out for two different geometry variants.We compared a low resolution geometry I to a very detailed description of the terrain (geometry II).Geometry I consists of a regular 20 m grid.The spatial resolution conforms to DTED Level 2, which can be achieved using spaceborne remote sensing techniques like InSAR (e.g.SRTM).Anyway, the deployed 20 m grid was derived from LiDAR via static grid reduction (low pass filtering).In a preprocessing step the coarse grid was augmented by additional breaklines from terrestrial survey to improve the description of the river banks and weir channels.By contrast, geometry II was derived from the precise 1 m LiDAR DTM-W as a hybrid computation grid (floodplain: unstructured grid, river bed: structured grid) as described above.The two different grid variants are illustrated in Fig. 5.
For the investigated 2 km river reach geometry I featured a total number of 6602 nodes (13 091 triangular or quadrilateral surface elements) whereas geometry II consists of 96 225 nodes (191 229 surface elements).The 1 m LiDAR DTM-W contains a total number of 1.16 million grid posts and 20 600 breakline points (building terrain intersections and eaves).Thus, for geometry I only 0.5% of the available data source are utilised to describe the terrain geometry and, therefore, many surface details are lost.For geometry II 8% of the LiDAR DTM-W posts are necessary to describe the terrain surface with a maximum approximation tolerance z max =±20 cm.
The most important results of hydraulic simulations are water depths, flow velocities and flow direction vectors.Based on these parameters water surface levels, flood extents and bottom shear stresses are derived.Figure 6 contains a map of water depths and flood boundaries for the same section and the same geometry variants as shown in Fig. 5.A remarkable difference of the flood extent between model I and II can be observed, especially close to the railway dam and the road dams.These geometric details are contained in model II only.They restrict the inundated area whereas the flooded area exceeds beyond the dams in model I.The consequences of the restricted run-off area in model II are illustrated in Fig. 7.It can be seen that the flow velocities in the main channel are higher in model II than in model I compensated by higher overbank flow velocities in model I compared to model II.
Additionally, overbank areas are often characterised by non-parallel water flow.2-D-HN modelling is the method of first choice for deriving distinct flow directions in each cell element.Figure 8 shows the flow vectors for a detailed area around the Dungl hotel (black window in Fig. 7b).It can clearly be seen that secondary currents around the hotel are well represented in model II whereas unrealistic parallel flow vectors even crossing the building can be observed in model I.Moreover, Fig. 8 reveals the above mentioned differences in overbank velocities more clearly.Model I fea-tures flow velocities within a range of 0.5-0.8m/s, whereas model II shows values between 0-0.2 m/s.
A more thorough quantitative comparison of the velocity differences of both model variants was carried out for a 215 m long section along the center line of the river Kamp close to the Dungl hotel.The profile outline is plotted in Fig. 7, whereas the velocity profile itself is illustrated in Fig. 9.The mean velocity in model I is 2.17 m/s (sample size n=45) and 3.17 m/s (n=97) in model II.The difference of the mean velocities between model I and II is statistically significant with respect to a 95% confidence interval (p=0.0). Figure 9  Finally, the cross section plotted in Fig. 7 was examined in detail and the results concerning water depth, flow velocity and bottom shear stress are illustrated in Fig. 10.The water depth of both model variants exhibit only small deviations in the range of a few centimeters.However, an increase of the flow velocity can be observed in model II within the main channel, which is compensated by lower flow velocities in the over bank areas.Model I shows a reciprocal progression with lower velocities in the main channel and higher values in the foreland.At station 150 m, for instance, the velocity adds up to 0.76 m/s in model I and only 0.21 m/s in II.Furthermore, an interesting detail can be observed at the location of the Dungl hotel.Since the buildings are contained in the computation grid of model II in a generalised way, the water depth, velocity and bottom shear stress go down to zero.These details are not reflected in the coarser model I.

Discussion
Due to the fact, that LiDAR based hydraulic modelling can be seen as state of the art in hazard zone mapping, it should be highlighted that precise and reliable modelling results (i.e.flood extents . . . ) can only be achieved by exploiting the full potential of the LiDAR data concerning accuracy and resolution.The progress in computing power and hydraulic modelling software enables the application of more complex geometries as provided by LiDAR.The results described in the previous section made clear, that the consideration of a more detailed geometry has significant impact on the hydraulic modelling results, not only concerning the flood extent but, above all, regarding the distribution of flow velocities and bottom shear stresses.The latter are equally important for estimating the risk potential in hazard zone mapping.In the following the results stated above are discussed critically in detail.

Discussion
Due to the fact, that LiDAR based hydraulic modelling can be seen as state of the art in hazard zone mapping, it should be highlighted that precise and reliable modelling results (i.e.flood extents . . . ) can only be achieved by exploiting the full The potential discharge area in the village Gars is restricted by the railway dam (right bank) and road dams (left bank).For the 100-year flood event the dams act as a flow barrier in model II since they are reflected sufficiently in the computation grid whereas the coarser geometry approximation of model I allows overtopping of the dams and, therefore, inundating of large flood plain areas beyond the dams.Proper consideration of the dams requires a certain ground resolution of the topographic data acquisition.A resolution of 20-30 m, that can be achieved using spaceborne In-SAR techniques is, therefore, not sufficient for precise hazard zone mapping.Airborne LiDAR, in contrast, is well suited both concerning the required planimetric resolution of approx. 1 m and height accuracy better than 20 cm.The relevant flow restricting features can be considered in the computation grid either by means of a higher point density at the base and top of the dams or by including the dam breaklines as constraints in the mesh.It has to be mentioned, that in case of the reported 2002 flood the entire area of Gars including the flood plains beyond the dams were inundated, but the discharge of this extreme event was by far higher than the HQ100 discharge for which the simulations at hand were carried out.Furthermore, it should be stressed that the dams were treated as watertight in the simulations, i.e. no culverts were considered.However, openings in dams are nowadays more and more safeguarded by mobile flood control measures.
The restriction of the discharge area by dams and buildings as provided by geometry II entails a discharge concentration in and near the main channel.As a consequence, higher flow velocities and bottom shear stresses can be observed there.The significant increase of the bottom shear stress (+63% compared to model I) represents a severe issue for the structural integrity of projected flood safety measures as well as existing buildings.Many cases of dam breaks and undermining of building foundations were reported in the recent past (Habersack and Krapesch, 2006).
The findings of the study clearly reveal that hydraulic characteristics (e.g.flow vectors) are less realistic in geometry I compared to the LiDAR derived geometry II.The crosssectional comparisons of flow velocity and bottom shear stress have shown a decrease of these parameters in the precise model variant II compared to the coarse model I for overbank areas.Maximum bottom shear stress values of up to 10 N/m 2 can be observed in model I, whereas the values are five times lower (approx.1-2 N/m 2 ) in II.Low bottom shear stress entails a sedimentation of transported suspended loads.Siltation processes, in turn, can be observed at virtually every flood event and provoke a high vulnerability particularly for residential areas.Therefore, a correct estimation of flow velocities and bottom shear stresses in hydraulic models can be regarded as equally important as a precise estimation of the flood extent.This becomes even more relevant when using flood hydraulics as boundaries for suspended-and bed load transport modelling.The comparison of the simulation results at hand has shown the high impact of the detailed Li-DAR based geometry description II especially on flow velocity and bottom shear stress.It can be stated that a correct representation of these hydraulic modelling parameters can only be achieved on the basis of a precise and detailed description of the geometry.

Conclusions
In this article a method for the generation of a hydraulic grid was presented.It starts from a digital terrain model in the form of a regular or hybrid grid, which is derived from airborne LiDAR point clouds by a method considering the random measurement errors.An irregular thinning scheme allows deriving an irregular triangular mesh considering a userdefined height tolerance compared to the original DTM.In the example presented, the point set describing the terrain was reduced to less than 8% in comparison to the original point cloud and the regular DTM grid.The obtained triangular mesh is conditioned by inserting vertices in order to assure that the hydrodynamic numerical simulations can be solved.
Airborne LiDAR allows identifying vegetated areas and buildings.The properties of LiDAR even enable retrieving terrain point measurements below the forest cover and, thus, enable obtaining accurate terrain models.Especially for low vegetation areas, full waveform laser scanning is a step forward to increase the reliability of the ground point classification.The identified building blocks, in turn, can be inserted into the grid for flood event simulation.
A flood event simulation on a grid according to the proposed method was compared to a simulation on a regular hydraulic grid featuring a resolution in the order of 20 m.Eleva-tion data with similar resolutions are available publicly (e.g.SRTM).New spaceborne satellite missions aim at resolutions of a few meter, but the accuracy of airborne LiDAR will not be reached.Additionally, as short wavelengths (X-band) are used, no terrain below forest canopy can be derived.
The two event simulations were similar with respect to their water surface levels but showed remarkable differences concerning the flooded area, flow velocity, etc. Comparisons to real floodings (e.g.suspended sediment deposition, side erosion processes) indicate, that the simulation on the basis of the low resolution elevation information is not realistic.
We conclude therefore that (i) airborne LiDAR is the method of first choice for topographic data acquisition of floodplains, and (ii) that the method of first reducing random measurement errors and subsequently applying irregular thinning in overbank areas and cell arrangement with respect to the flow direction in the river bed leads to meshes, that need little conditioning such models are well suited as input for HN analysis.This asserts that the geometric information captured by airborne LiDAR is in its essence available for hydrodynamic-numerical simulations.

Fig. 3 .
Fig. 3. Perspecitve view of the final 1 m DTM-W of Gars am Kamp including the Dungl hotel (image center), buildings are considered as block models and integrated as breaklines into the hybrid grid.

Fig. 4 .
Fig. 4.Perspective view of a professionally thinned and conditioned DTM-W grid, which also serves as hybrid hydraulic computation mesh (c.f.geometry II); overbank area: irregular distribution of mesh nodes (i.e.unstructured grid), river bed and river bank: cells aligned to the principal flow direction (i.e.structured grid).

Fig. 5 .Fig. 6 :
Fig. 5. Computation grid of two geometry variants I and II for HN modelling, Gars/Kamp, Lower Austria; (a) Geometry I: regular 20×20 m 2 grid, additional breaklines for description of river bank and weir channel, river bed: structured grid; (b) Geometry II: hybrid grid derived from 1 m LiDAR DTM-W via adaptive TIN refinement and subsequent mesh conditioning, river bed and river bank: structured grid (post spacing 5 m in flow direction and 2 m perpendicular to it), over bank area: unstructured grid: maximum height tolerance z max =±20 cm, buildings in red.

Fig. 6 .
Fig. 6.Flood extents (red boundary line) and water [m] for a simulated HQ100 (470 m 3 /s) in Gars/Kamp; (a) extended flooding of residential area beyond the dams in model I; (b) inundation area restricted by dams (e.g.railway dam) in model II.
illustrates the progression of the flow velocity and bottom shear stress (i.e.parameter for erosion and deposition of sediments) along the longitudinal section.Both parameters show a similar progression in both geometry variants with greater values in model II than in model I.The maximum bottom shear stress values are 129 N/m 2 in I and 211 N/m 2 in II which represents a relative increase of +63% of model II compared to I.

Fig. 8 .Fig. 9 .Fig. 10 .
Fig. 8. Flow vectors of a detailed area around the Dungle hotel for a simulated HQ100 (470 m 3 /s) in Gars/Kamp; (a) model I shows unrealistic parallel flow vectors within and around the building; (b) flow vectors of model II reflect the secondary currents around the building correctly.