Flood Damage Analysis : First Floor Elevation Uncertainty Resulting from LiDAR-Derived Digital Surface Models

The use of high resolution ground-based light detection and ranging (LiDAR) datasets provides spatial density and vertical precision for obtaining highly accurate Digital Surface Models (DSMs). As a result, the reliability of flood damage analysis has improved significantly, owing to the increased accuracy of hydrodynamic models. In addition, considerable error reduction has been achieved in the estimation of first floor elevation, which is a critical parameter for determining structural and content damages in buildings. However, as with any discrete measurement technique, LiDAR data contain object space ambiguities, especially in urban areas where the presence of buildings and the floodplain gives rise to a highly complex landscape that is largely corrected by using ancillary information based on the addition of breaklines to a triangulated irregular network (TIN). The present study provides a methodological approach for assessing uncertainty regarding first floor elevation. This is based on: (i) generation an urban TIN from LiDAR data with a density of 0.5 points ̈ m ́2, complemented with the river bathymetry obtained from a field survey with a density of 0.3 points ̈ m ́2. The TIN was subsequently improved by adding breaklines and was finally transformed to a raster with a spatial resolution of 2 m; (ii) implementation of a two-dimensional (2D) hydrodynamic model based on the 500-year flood return period. The high resolution DSM obtained in the previous step, facilitated addressing the modelling, since it represented suitable urban features influencing hydraulics (e.g., streets and buildings); and (iii) determination of first floor elevation uncertainty within the 500-year flood zone by performing Monte Carlo simulations based on geostatistics and 1997 control elevation points in order to assess error. Deviations in first floor elevation (average: 0.56 m and standard deviation: 0.33 m) show that this parameter has to be neatly characterized in order to obtain reliable assessments of flood damage assessments and implement realistic risk management.


Introduction
The European Directive on the assessment and management of flood risk (Directive 2007/60/EC) recognized for the first time at the European level that floods have the capacity to severely compromise economic development and undermine the economic activities of the community [1].From an economic point of view, the main purpose of flood risk management (FRM) is to reduce flood damage and this requires as much accurate elevation data as possible in order to obtain reliable digital surface models (DSMs), i.e., information on the height of the earth surface, including the objects on it, such as buildings, trees, bridges, etc.This is a critical issue in urban areas due to the presence of man-made structures (e.g., buildings, levees, roads, bridges, culverts) that make the land surface irregular, which may significantly modify spatio-temporal flood configuration and induce variable flow resistance characteristics.Another critical variable concerning FRM and related to topography is the first floor elevation of structures (i.e., the elevation at which floodwater enters the building) [2].It is essential for estimating flood damage based on depth-damage curves, which can be established either from field surveys or by referring to LiDAR.
Airborne light detection and ranging or laser induced direction and ranging (LiDAR), also referred to as airborne laser scanning (ALS), provides high resolution topographic datasets.This has enabled LiDAR-derived DSMs to be used for carrying out flood hazard analysis for both coastal and fluvial urban areas [3][4][5].Currently, LiDAR is the most accurate approach for obtaining reliable topographic data on a large-scale [6,7].Though the first laser instruments were built in the 1960s [8,9], it was not until 1993 that the first prototype of ALS was commercialized and applied for both technical and scientific purposes [6].Since then, LiDAR has been extensively used in the environmental sciences [10][11][12][13][14][15].
LiDAR has been used on a basin scale to characterize spatial and temporal changes regarding erosion and deposition [16], or to extract fluvial landforms [17].As far as flood mapping is concerned, high-resolution DSMs derived from LiDAR have increasingly been used to improve the performance of one-dimensional (1D) [18] and two-dimensional (2D) hydraulic models, which are suitable for use in urban areas where buildings and other man-made structures condition flood wave behaviour, with the result that for practical purposes, flow cannot be assumed to be 1D [5,19].In addition, LiDAR provides very valuable data for running 2D hydrodynamic models.Specifically, high resolution digital terrain models (DTMs) derived from LiDAR enable grain-scale surface roughness, an essential parameter in flood modelling, to be estimated [20].They also provide a very accurate topographic dataset of both the ground and the urban environment (e.g., buildings, roads, bridges).In the last stage of the hydrodynamic 2D approach, heights linked to the DSM are assigned to the elements of the mesh in which the Saint Venant equations are solved [21,22].
The spatial resolution of DSMs plays an important role in terms of accurateness of flood mapping in urban areas.DSMs derived from data sources other than LiDAR (e.g., national topographic maps) usually have coarser spatial resolution and lower vertical accuracy.As a result, flooded areas tend to be greater when compared with the results obtained through the use of DSMs, which has finer spatial resolution [23].This explains why LiDAR data predominate in flood assessment, as the higher the spatial resolution of the topographic data, the higher the reliability of flood mapping.The main limitation of LiDAR technology is related to its inability to characterize bathymetry in channels for practical purposes, due to the limited capacity of the signal for penetrating the water column, especially if the river flow has a suspended sediment load [6].Therefore, bathymetry data incorporated into DSMs are usually derived from topographic data sources other than LiDAR (e.g., echo soundings and electronic theodolite surveys) [24].
Less attention has been paid to the use of LiDAR data during the flood risk assessment process.This data source has been used to map the risk of flooding from storm surges and in sea level rise in coastal areas based on the characterization of socio-economic and ecosystem impacts [25].In the fluvial environment, LiDAR has been used to improve risk management strategies by considering flood scenarios [26] or assessing socio-economic impacts [27].However, uncertainty analysis is not addressed in the approaches cited above, even though risk involves exposure to the chance of loss [28].In this regard, estimation of flood damage using depth-damage relationships requires specification of the first floor elevation of the structure [29].This elevation may be obtained from LiDAR, which basically has two sources of error.The first is related to the accuracy of the method [30,31].The second is a consequence of transformations in the original raw LiDAR point cloud (e.g., classification, thinning, interpolation and breakline enforcement) in order to obtain a hydraulically consistent DSM [5].To this end, the LiDAR data should preserve all geometric elements as precisely as possible, especially if these have a man-made origin.
The aim of this paper is to characterize uncertainty in the first floor elevation of buildings prone to being affected by floods and located in communities where lowest floor elevations are at ground level.To this end, a 2D hydrodynamic model was used to obtain the 500-year flood zone in an urban environment.Next, a geostatistical approach was put into practice to determine the spatial distribution of errors between the DSM derived from LiDAR and first floor control points in buildings.The Monte Carlo method was then used to describe errors with a probability density function (PDF).

Study Area
The municipality of Navaluenga is located in Central Spain on the banks of the Alberche River, between the Sierra del Valle (eastern range of Gredos, Spanish Central System) and the Sierra de la Paramera (40 ˝24'30"N; 4 ˝42'17"W; 761 m a.s.l.; Figure 1).
consistent DSM [5].To this end, the LiDAR data should preserve all geometric elements as precisely as possible, especially if these have a man-made origin.
The aim of this paper is to characterize uncertainty in the first floor elevation of buildings prone to being affected by floods and located in communities where lowest floor elevations are at ground level.To this end, a 2D hydrodynamic model was used to obtain the 500-year flood zone in an urban environment.Next, a geostatistical approach was put into practice to determine the spatial distribution of errors between the DSM derived from LiDAR and first floor control points in buildings.The Monte Carlo method was then used to describe errors with a probability density function (PDF).

Study Area
The municipality of Navaluenga is located in Central Spain on the banks of the Alberche River, between the Sierra del Valle (eastern range of Gredos, Spanish Central System) and the Sierra de la Paramera (40°24′30″N; 4°42′17″W; 761 m a.s.l.; Figure 1).
According to the National Statistics Institute of Spain, Navaluenga has a population of 2027 inhabitants (data corresponding to 2014), though there is a significant increase in population during the summer months, when the number of inhabitants rises to roughly 20,000.The total housing in the municipality is estimated to be 4311 dwellings, of which 3392 are considered to be second homes (the vast majority of them have the lowest floors located at the ground level).The Alberche River rises at approximately 1800 m a.s.l. and flows for 70 km before reaching Navaluenga.Its time of concentration (Tc) is around 8.5 h and it drains an area of 717 km 2 .In this reach the flow regime is totally natural.The Alberche has several weirs, i.e., natural pools that are used by local people as recreation areas in summer.In this area, several torrents flow into the Alberche River.Especially noteworthy is the Chorrerón Stream (Tc ≈ 3 h), which crosses the urban area perpendicularly from north to south to flow into the Alberche (Figure 1).According to the National Statistics Institute of Spain, Navaluenga has a population of 2027 inhabitants (data corresponding to 2014), though there is a significant increase in population during the summer months, when the number of inhabitants rises to roughly 20,000.The total housing in the municipality is estimated to be 4311 dwellings, of which 3392 are considered to be second homes (the vast majority of them have the lowest floors located at the ground level).The Alberche River rises at approximately 1800 m a.s.l. and flows for 70 km before reaching Navaluenga.Its time of concentration (Tc) is around 8.5 h and it drains an area of 717 km 2 .In this reach the flow regime is totally natural.The Alberche has several weirs, i.e., natural pools that are used by local people as recreation areas in summer.In this area, several torrents flow into the Alberche River.Especially noteworthy is the Chorrerón Stream (Tc « 3 h), which crosses the urban area perpendicularly from north to south to flow into the Alberche (Figure 1).
The village of Navaluenga has suffered flood events linked to the Alberche River and the Chorrerón Stream since at least the Early Middle Ages, as attested by documentary records that existed in the late fifteenth century.The most recent events in the 1990s and the early 2000s caused important economic losses and put part of the local population at risk [32].

DSM Construction
The LiDAR data used were provided by the Spanish National Geographic Institute (IGN) [33].Raw LiDAR data points were derived from airborne LiDAR systems with a density of 0.5 points¨m ´2.The altimetric precision (Z) obtained was therefore greater than 20 cm, taking the root mean square error (RMSE) as a statistic.The resulting point clouds were subsequently automatically classified by the IGN.The reference system for this topographic dataset is the ETRS89/UTM zone 30N (compatible with WGS 84).
One of the critical steps for generating DSM from LiDAR data is to separate the LiDAR points into ground (terrain) and non-ground (man-made structures).The point cloud was therefore filtered into two main categories: ground and buildings.This was done using LAStools [34] within the ArcGIS 10.1 environment.Regarding bathymetry, a field survey was carried out using a differential GPS (Trimble 5700).Bathymetry was performed with an average density of ~0.3 points¨m ´2.
With respect to the production of DSM for use in flood modelling, mass points derived from LiDAR alone do not adequately capture changes in topographic breaks (e.g., riverbanks, roads, streets buildings).Significant topographic changes in terrain slope can be represented as breaklines [35].These complement the mass points and have the role of enforcing the topographic breaks of a DSM within triangulated irregular networks (TINs).From the number of possible approaches for treating breaklines in LiDAR datasets, we chose local knowledge, GIS, orthophotos and cadastral mapping as the most appropriate for identifying the location of critical linear features.We then took the average z-values from LAS data to obtain 3D breakline strings [36].With respect to 2D modelling, the breaklines enabled the floodplain (i.e., river banks) and urban areas (i.e., roads and streets) to be defined.Because of the critical hydraulic importance of these features, they should be specifically included in the DSM.
Breaklines were manually digitized to suitably represent the banks of both the Alberche River and the Chorrerón Stream.In urban areas, the streets obtained from cadastral mapping were used as breaklines.The original 2D features were converted to 3D by adding the 2D features and LAS files to the Arcmap's 3D analyst extension, taking the medium height of the LiDAR points intersected by the 2D features as reference.To minimize errors, streets were segmented, using the sections delimited by slope breaks as criteria.The TIN was edited with the point features (i.e., ground and buildings derived from LiDAR) as mass points.The line features were integrated as breaklines in the case of streets and river banks, while isobaths were added as softlines.In the final step, the TIN was transformed to a raster by applying the linear method as interpolation choice.The resulting digital elevation model (DEM) was defined with a spatial resolution of 2 m in order to be consistent with the point spacing of the LiDAR data.This DEM enabled both accomplishing 2D hydrodynamic modelling and estimating uncertainty in first floor elevation.

2D Hydrodynamic Modelling
Hydrodynamic simulation was carried out by applying Iber two-dimensional hydrodynamic software [37].Iber is a numerical tool for 2D simulation of unsteady turbulent free surface flow and sediment transport in water-courses.To solve the hydrodynamics, Iber uses the finite volume method.This method is suitable for flows in mountainous rivers [38], where shocks and discontinuities may occur, giving very sharp hydrographs.
The hydrodynamic model was based on a rectangular triangulated irregular network (RTIN) obtained from the DSM of the study site.The maximum and minimum lengths of the sides of the triangles were 4 and 2 m, respectively.Tolerance (maximum vertical distance between the DSM and the geometry created) was considered to be 0.2 m, which coincides with the altimetric accuracy of the LiDAR data.With respect to boundary conditions, discharges corresponding to the return period of 500 year for the Alberche River (i.e., 2,006 m 3 ¨s´1 ) and the Chorrerón Stream (i.e., 167 m 3 ¨s´1 ) were considered.The roughness coefficient was obtained from official land cover mapping on a scale of 1:25,000 [39] and then the corresponding value of Manning's n was assigned to every land cover unit.
Respecting hydraulic infrastructures (see Figure 1), there are 7 culverts along the reaches studied.To calculate flow transfer through them, Iber uses the Manning's equation [40].Likewise, there are two bridges and a weir in the reach of the Alberche that crosses Navaluenga.Iber simulates bridges as internal conditions by changing the equations with which flow is calculated on the edges of the elements affected.The internal condition was defined in the upstream cross sections of the bridges.The discharge coefficients used were 0.6 for free flow pressure under deck and 0.8 for flooded flow pressure under deck.The discharge coefficient chosen in the case of over-deck flow was 1.7.Finally, the weir discharge coefficient was also taken as 1.7.
As additional general information on the model, a second-order numerical scheme was used, which was more accurate than that of the first order [41].In addition, 0.45 was considered as the Courant-Friedichs-Lewy (CFL) condition and 0.01 m as the wet-dry limit.As regards the drying method, the default option provided by Iber was used.This is a robust method that preserves the water mass throughout the domain and has the advantage that the calculation time virtually does not depend on the dry-wet process.

Geostatistical Uncertainty Analysis
First floor errors in the urban area of the LiDAR-derived DSM were assessed from 603 elevation control points with an error lower than ˘15 cm, provided by the Spanish cadastre [42].The remaining 1294 points (with the same vertical accuracy as the one mentioned above) were located away from the streets.To assess the errors, the control points were subtracted from the DSM.The spatial distribution of the first floor errors in the urban area was explored by performing Monte Carlo simulations based on geostatistics.Each simulation resulted in a first floor error map in which the error at control points never changed in the simulation process (conditioning the simulations by known errors).The stochastic simulations provided a series of random plausible maps that could be used to represent uncertainty about true first floor elevation error [43].
Sequential Gaussian Simulation (SGS) was chosen as this is the most common technique for generating conditional stochastic simulations for spatially continuous variables.SGS uses the Kriging geostatistical interpolator to simulate values [44].Kriging is the generic name for a family of generalized least-squares regression algorithms based on regionalized variable theory that assume that spatial variation of the variable is statistically homogeneous throughout the region [45].In contrast to deterministic interpolators, Kriging uses a model of the spatial correlation or structure of processes.The spatial structures are characterized by a variogram, which is estimated from the sampled data.The variogram is then used to estimate the Kriging weights used for data interpolation.The spatial structure of the first floor error in the urban area was modelled using a spherical variogram.SGS is conditional in the sense that it will honour the variogram model as well as the existing observations.SGS starts by defining a random path for visiting each node of the grid (i.e., 2 m of spatial resolution) once.At each node of the dataset, a specified number of original observations and previously simulated nodes are selected for conditioning purposes, and Kriging is used to determine the locations-specific mean and variance of the conditional cumulative distribution function (CCDF).Finally, a random value is drawn from the CCDF; the value is added to the dataset and the procedure is repeated until all nodes have been visited [44].
Variogram modelling and the simulation process were performed entirely with the Stanford Geostatistical Modelling Software (SGeMS) package [46].A total of 1000 realizations of first floor elevation errors in the urban area were generated to guarantee a stable result [47].Uncertainty was computed by evaluating the statistics associated with the range of simulated errors.These calculations allowed mapping of the first floor error mean and standard deviation and the probability of exceeding an error threshold value of ˘0.5 m.

500-Year Flood Zone
The river Alberche has a steeper lateral slope along its right bank than on its left bank.This causes a tendency for the river to overflow mostly on its left bank.This overflow merges with flooding from the Chorrerón Stream.Moreover, most of the population is located on the left bank of the river, which puts it at significant risk of being affected.Given the low capacity of the channel for the return period studied (500 years), water overflows even before it enters the village.The overflow situation is made worse by the obstruction caused by the first medieval bridge located at the entrance to the village (B1, see Figure 1).This creates a hydraulic depth of 9.31 meters on the left bank in the village and an additional flow parallel to the channel.The situation is exacerbated by a more recently built second bridge a hundred meters downstream (B2, see Figure 1), which contributes to these high depths being maintained (Figure 2).

500-Year Flood Zone
The river Alberche has a steeper lateral slope along its right bank than on its left bank.This causes a tendency for the river to overflow mostly on its left bank.This overflow merges with flooding from the Chorrerón Stream.Moreover, most of the population is located on the left bank of the river, which puts it at significant risk of being affected.Given the low capacity of the channel for the return period studied (500 years), water overflows even before it enters the village.The overflow situation is made worse by the obstruction caused by the first medieval bridge located at the entrance to the village (B1, see Figure 1).This creates a hydraulic depth of 9.31 meters on the left bank in the village and an additional flow parallel to the channel.The situation is exacerbated by a more recently built second bridge a hundred meters downstream (B2, see Figure 1), which contributes to these high depths being maintained (Figure 2).
After the second bridge, water overflows before the entry of the Chorrerón, due to the low capacity of the river channel and its urban culverts.Therefore, both overflow situations combine at the confluence.Upstream from the mouth of the Chorrerón, the minimum outfall cross-section for its expected flow (167 m 3 •s −1 for the 500-year flood period) is 4 × 0.5 m 2 .The confluence of the Chorrerón and the Alberche also plays an important role, especially if both are simultaneously provoking a flood as in the scenario described here.Downstream of the Chorrerón, the existing weir partially maintains this level and also accelerates the flow downstream.This situation does not cause a reduction of water depths due to the fact that the section of the river channel is reduced on its last stretch.The flow generated on the left margin generates a by-pass in the populated area that avoids the weir.Flow subsequently merges into the river Alberche downstream.As a result of this, the highest depth in the channel is 10.39 m, while in the flooded area it is 9.31 m.Outside the area of influence of the points of conflict described above, maximum depths are 8.15 m in the channel and 7.78 m in the flooded area.In the urban area, depths are between 0.01 m and 2.6 m.Fluctuations in flow velocities occur where there are changes in the flow regime (i.e., from subcritical to supercritical).These are mainly due to narrowing of the channel cross-sections caused by the bridges on the urban reach of the Alberche.In the river channel, maximum flow velocities of 9.81 m•s −1 are given, while the minimum velocities (i.e., 1.10 m•s −1 ) are found close to the river banks After the second bridge, water overflows before the entry of the Chorrerón, due to the low capacity of the river channel and its urban culverts.Therefore, both overflow situations combine at the confluence.Upstream from the mouth of the Chorrerón, the minimum outfall cross-section for its expected flow (167 m 3 ¨s´1 for the 500-year flood period) is 4 ˆ0.5 m 2 .The confluence of the Chorrerón and the Alberche also plays an important role, especially if both are simultaneously provoking a flood as in the scenario described here.Downstream of the Chorrerón, the existing weir partially maintains this level and also accelerates the flow downstream.This situation does not cause a reduction of water depths due to the fact that the section of the river channel is reduced on its last stretch.The flow generated on the left margin generates a by-pass in the populated area that avoids the weir.Flow subsequently merges into the river Alberche downstream.
As a result of this, the highest depth in the channel is 10.39 m, while in the flooded area it is 9.31 m.Outside the area of influence of the points of conflict described above, maximum depths are 8.15 m in the channel and 7.78 m in the flooded area.In the urban area, depths are between 0.01 m and 2.6 m.Fluctuations in flow velocities occur where there are changes in the flow regime (i.e., from subcritical to supercritical).These are mainly due to narrowing of the channel cross-sections caused by the bridges on the urban reach of the Alberche.In the river channel, maximum flow velocities of 9.81 m¨s ´1 are given, while the minimum velocities (i.e., 1.10 m¨s ´1) are found close to the river banks (mean = 3.1 m¨s ´1; standard deviation = 1.71 m¨s ´1).Flow velocities within the urban area are between 0.02 m¨s ´1 and 1.99 m¨s ´1 (mean = 0.43 m¨s ´1; standard deviation = 0.44 m¨s ´1).With regard to Froude numbers, values are high (i.e., above 2.5) for both the Alberche River and the Chorrerón Stream.Froude numbers in the flooded areas present lower values, though these are slightly higher on the left bank of the Alberche because of the preferential flow in that area (Figure 2).

First Floor Elevation Uncertainty
The distribution of first floor elevation errors was slightly biased towards positive values when analysed for the whole study site.When it was defined by the 500-year flood zone only, errors were biased towards negative values.Average first floor elevation values were similar in both analyses (i.e., 0.54 ˘0.32 m in the municipality as a whole and 0.56 ˘0.33 m in the 500-year flood zone).With respect to ranges, errors varied from ´2.2 m to 3.6 m in the first case, and from ´2.2 m to 1.4 m when the analysis was focused on the area prone to flooding (Figure 3).In the non-urban area, mean elevation errors were slightly lower compared to measurements of first floor elevation errors (i.e., mean value = 0.48; standard deviation = 0.45 m), with a range between ´2.89 m and 2.22 m and a standard deviation of 0.45 m.If the data are analysed jointly, it appears that errors vary depending on the type of land cover (Figure 4).Therefore, the highest errors occur in urbanized areas (mean = 0.56m; standard deviation = 0.38) while the lowest ones appear in forests and woodlands (mean = 0.19 m; standard deviation = 0.56 m).Froude numbers in the flooded areas present lower values, though these are slightly higher on the left bank of the Alberche because of the preferential flow in that area (Figure 2).

First Floor Elevation Uncertainty
The distribution of first floor elevation errors was slightly biased towards positive values when analysed for the whole study site.When it was defined by the 500-year flood zone only, errors were biased towards negative values.Average first floor elevation values were similar in both analyses (i.e., 0.54 ± 0.32 m in the municipality as a whole and 0.56 ± 0.33 m in the 500-year flood zone).With respect to ranges, errors varied from −2.2 m to 3.6 m in the first case, and from −2.2 m to 1.4 m when the analysis was focused on the area prone to flooding (Figure 3).In the non-urban area, mean elevation errors were slightly lower compared to measurements of first floor elevation errors (i.e., mean value = 0.48; standard deviation = 0.45 m), with a range between −2.89 m and 2.22 m and a standard deviation of 0.45 m.If the data are analysed jointly, it appears that errors vary depending on the type of land cover (Figure 4).Therefore, the highest errors occur in urbanized areas (mean = 0.56m; standard deviation = 0.38) while the lowest ones appear in forests and woodlands (mean = 0.19 m; standard deviation = 0.56 m).We have characterized the spatial structure of first floor elevation errors with the variogram shown in Figure 5.The variogram illustrates the variance between pairs of values (y-axis) at different spacing (x-axis), and, most importantly, increasing variance with the spacing of data pairs.Regarding the above, no anisotropy was detected in the errors, and the variogram revealed that errors appear to become spatially unstructured at lag distances >122 m.It was possible to fit a theoretical variogram to the experimental variogram.Visual interface was used for variogram curve fitting and the best fit model variogram found is shown in Figure 5.
As regards uncertainty, the magnitude and spatial distribution of error in the DSM elevations values were evaluated in relation to the elevation control points in order to assess the potential effect of DSM error on first floor elevation.The uncertainty error analysis was carried out using the Monte Carlo method, simulating 1000 urban DSM error images with the variogram modelled to the experimental errors.This method offered flexibility at the cost of heavy computational load (almost 24 h were required to generate the 1000 simulations).The set of 1000 simulations established the bounds within which it can be affirmed that the true map lies (Figure 6).
All the realizations show patterns of elevation error that are similar on a gross scale but differ in details from map to map.The large-scale similarities reflect the spatial distribution of measured errors, which persists from simulation to simulation, whereas the small-scale differences reflect the randomness added by Gaussian Simulation of values between the points at which elevation error was measured.We have characterized the spatial structure of first floor elevation errors with the variogram shown in Figure 5.The variogram illustrates the variance between pairs of values (y-axis) at different spacing (x-axis), and, most importantly, increasing variance with the spacing of data pairs.Regarding the above, no anisotropy was detected in the errors, and the variogram revealed that errors appear to become spatially unstructured at lag distances >122 m.It was possible to fit a theoretical variogram to the experimental variogram.Visual interface was used for variogram curve fitting and the best fit model variogram found is shown in Figure 5.
As regards uncertainty, the magnitude and spatial distribution of error in the DSM elevations values were evaluated in relation to the elevation control points in order to assess the potential effect of DSM error on first floor elevation.The uncertainty error analysis was carried out using the Monte Carlo method, simulating 1000 urban DSM error images with the variogram modelled to the experimental errors.This method offered flexibility at the cost of heavy computational load (almost 24 h were required to generate the 1000 simulations).The set of 1000 simulations established the bounds within which it can be affirmed that the true map lies (Figure 6).
All the realizations show patterns of elevation error that are similar on a gross scale but differ in details from map to map.The large-scale similarities reflect the spatial distribution of measured errors, which persists from simulation to simulation, whereas the small-scale differences reflect the randomness added by Gaussian Simulation of values between the points at which elevation error was measured.The average probability of first floor elevation errors within the ˘0.5 m range was 0.57, and 0.52 when the analysis took the 500-year flood zone into account (Figure 7).This probability increased to 0.75 when both first floor elevation and elevation errors measured in the floodplain and other non-urbanized areas were included in the 500-year flood zone.In terms of spatial representativeness, 68% of the street surfaces within the 500-year flood zone had a 0.4 probability of first floor elevation being between ˘0.5 m.This percentage increased to 77% when the whole street area was considered (Table 1 The average probability of first floor elevation errors within the ±0.5 m range was 0.57, and 0.52 when the analysis took the 500-year flood zone into account (Figure 7).This probability increased to 0.75 when both first floor elevation and elevation errors measured in the floodplain and other non-urbanized areas were included in the 500-year flood zone.In terms of spatial representativeness, 68% of the street surfaces within the 500-year flood zone had a 0.4 probability of first floor elevation being between ±0.5 m.This percentage increased to 77% when the whole street area was considered (Table 1).Table 1.Area occupied by deciles of probability considering that first floor elevation error is between −0.5 m and 0.5 m.

2D Hydrodynamic Modelling
One-dimensional (1D) hydraulic models are unable to represent and correctly simulate physical and hydrodynamic conditions, which are critical for understanding systems as complex as urban areas.Instead, 2D hydraulic modelling is a much more consistent approach because of its ability to

2D Hydrodynamic Modelling
One-dimensional (1D) hydraulic models are unable to represent and correctly simulate physical and hydrodynamic conditions, which are critical for understanding systems as complex as urban areas.Instead, 2D hydraulic modelling is a much more consistent approach because of its ability to represent overbank flow [48] and high topography complexity, even in urban areas where flow can be represented on the scale of individual buildings [49].
2D models require the geometric characterization of bathymetry and its surrounding area as a continuous surface in order to obtain reliable flood inundation extents [50].The precision of digital topography, expressed both in terms of spatial resolution and vertical accuracy, is a key factor for flood mapping.Specifically, in urban areas, high-resolution topographic data enable flood modelling based on a 2D approach to be addressed.This is because the topography of man-made structures is adequately represented, thanks to advances in airborne LiDAR systems, among other topographic data sources, which make it possible to acquire high quality terrain data.However, this source of topographic data by itself is not able to properly represent fluvial systems that are partially occupied by urban areas, where continuous and linear topographic features associated with significant changes in gradient are quite common.The LiDAR data used here have a relatively high horizontal density of points (i.e., 0.5 points¨m ´2).However, continuous linear features such as river banks, and especially streets, would only be resolved if the data had been sampled to a higher spatial resolution, which, had this been possible, would have resulted in large and unwieldy files.
As a result, a DSM obtained from LiDAR data only was not sufficiently reliable for use in a 2D hydrodynamic model.The addition of breaklines (also defined as structure lines or skeleton lines) was therefore required (Figure 8).These constitute linear features that describe surface changes [51].Their integration in the DSM of the study site contributed significantly to obtaining a reliable, morphologically correct and enhanced DSM [52].
Remote Sens. 2016, 8, 604 11 of 17 represent overbank flow [48] and high topography complexity, even in urban areas where flow can be represented on the scale of individual buildings [49].2D models require the geometric characterization of bathymetry and its surrounding area as a continuous surface in order to obtain reliable flood inundation extents [50].The precision of digital topography, expressed both in terms of spatial resolution and vertical accuracy, is a key factor for flood mapping.Specifically, in urban areas, high-resolution topographic data enable flood modelling based on a 2D approach to be addressed.This is because the topography of man-made structures is adequately represented, thanks to advances in airborne LiDAR systems, among other topographic data sources, which make it possible to acquire high quality terrain data.However, this source of topographic data by itself is not able to properly represent fluvial systems that are partially occupied by urban areas, where continuous and linear topographic features associated with significant changes in gradient are quite common.The LiDAR data used here have a relatively high horizontal density of points (i.e., 0.5 points•m −2 ).However, continuous linear features such as river banks, and especially streets, would only be resolved if the data had been sampled to a higher spatial resolution, which, had this been possible, would have resulted in large and unwieldy files.
As a result, a DSM obtained from LiDAR data only was not sufficiently reliable for use in a 2D hydrodynamic model.The addition of breaklines (also defined as structure lines or skeleton lines) was therefore required (Figure 8).These constitute linear features that describe surface changes [51].Their integration in the DSM of the study site contributed significantly to obtaining a reliable, morphologically correct and enhanced DSM [52].In general, the accuracy or current errors of LIDAR data are not relevant for many applications (e.g., cadastral work, hydrological studies, agronomy and urban planning).However, in hydraulic modeling, small differences in geometry can generate significant changes in levels and flow rates, due to the basic equations of fluid mechanics.These changes are not localized but will spread to In general, the accuracy or current errors of LIDAR data are not relevant for many applications (e.g., cadastral work, hydrological studies, agronomy and urban planning).However, in hydraulic modeling, small differences in geometry can generate significant changes in levels and flow rates, due to the basic equations of fluid mechanics.These changes are not localized but will spread to other areas of the territory and will affect the results of the hydraulic model.Traditionally, cross-sections for studies on flooding were done with traditional topographical instruments such as the theodolite.These cross-sections were obtained "in situ" from traditional surveying, with expert criteria applied directly to the territory.This way of working required visiting the study area and looking for culverts, walls or local constrictions that could affect hydraulic simulation.This involved several days' work, sometimes with accessibility problems near the river.
Theoretically, LiDAR data provide good accuracy and questions are not usually raised about their reliability, or whether they accurately represent reality.However, LiDAR may not contain important details from a hydraulic point of view because of its systematic sampling procedure.There is much room for improvement in the representation of many elements of LiDAR-derived DSMs, despite the high density of points represented.For example, LiDAR-derived DSMs are unable to reflect thin structures, like levees or continuous walls, which are an obstacle to water.Also, they do not allow small structures such as irrigation ditches to be represented, possibly because such information has not been captured correctly.
These structures can provoke overflows and varying hydraulic variables, such as water depth or velocity.It is therefore necessary to include them in the model so that hydraulic outputs can be an accurate representation of reality.The importance of making corrections to LiDAR-derived DSMs is justified to ensure they correctly represent both the natural and man-made physical geometries that influence hydraulic modelling.The case study presented here is a clear example of this effect.The confluence of the Alberche River and the Chorrerón Stream is located within the urban area.Moreover, the existence of hydraulic structures as bridges and culverts directly contributes to flooding.Small differences in the topography of these conflictive points will affect results by changing water surface elevation or velocity.As a final state, it is recommended to take advantage of LiDAR data, but in combination with classic topography at conflictive points, enabling the suitable representation of elements that may be essential to flow behaviour in the final DSM and thereby reducing errors in flood assessment.

Geostatistical Analysis
Flood risk implies exposure to the possibility of injury or loss.As a result, there is a need for uncertainty to be both described and dealt with, a situation that is highly complex, since numerous sources of uncertainty have to be characterized, specifically: (i) uncertainty of return periods in flood frequency analysis; (ii) uncertainty in flood mapping as a result of using simplified hydraulic models, as well as errors in DSM and in the parameterizing of roughness; and (iii) uncertainty in flood damage analysis, due to the paucity of information on the relationship between depth and inundation damage and the lack of accuracy in estimating structure and contents values [53].
Therefore, estimating first floor elevation constitutes an additional source of uncertainty to be considered within this complex scheme.This factor plays a critical role in flood damage analysis, as it determines the depth-damage relationship and, as a consequence, economic losses in buildings [54].Basically, uncertainty related to characterization of first floor elevation is tied to survey errors, errors in the interpolation of the LiDAR point cloud, and simplifications inherent in adding breaklines to urban areas and the floodplain in order to obtain a consistent DSM from a hydraulic point of view.
In the approach presented here, we hypothesize that errors in the non-urban area are lower than in the urban area, since there is no forced transformation of the DSM by adding breaklines.In addition, errors also show differences that depend on the land cover (Figure 9).Therefore, measurements of average errors in non-urban areas are similar to those obtained by other authors, e.g., [55].Positive outliers detected in the sampling area may result from suspended objects being hit by laser beams.As regards negative outliers, their presence could be due to laser beams being reflected among buildings several times before they are detected.These specular reflections result in the laser beam having a longer travel time, and a lower elevation is therefore calculated during post-flight processing [56].In the urban area, with particular reference to the first floor elevation of buildings, estimates of average errors, standard deviations and range are slightly higher than those measured in the non-urban area.This is because LiDAR derived-DSM presents as additional source of error derived from adding breaklines.Incorporating breaklines into the TIN is based on assuming that breaklines adopt as the z-value the average height of the LiDAR point cloud intersected by a given breakline.In order to reduce the resulting uncertainty, streets were divided into segments.The existence of breaks in the slope were taken as criteria, which enabled the error to be kept to a minimum but not overridden.
Remote Sens. 2016, 8, 604 13 of 17 measured in the non-urban area.This is because LiDAR derived-DSM presents as additional source of error derived from adding breaklines.Incorporating breaklines into the TIN is based on assuming that breaklines adopt as the z-value the average height of the LiDAR point cloud intersected by a given breakline.In order to reduce the resulting uncertainty, streets were divided into segments.The existence of breaks in the slope were taken as criteria, which enabled the error to be kept to a minimum but not overridden.The advantage of using a geostatistical method is that it allows errors associated with the DSM values to be studied and the spatial structure to be taken into account by using the variogram.Also of benefit is the use of geostatistical simulations of equiprobable spatial distribution to determine uncertainty in the distribution of error values.It is possible to assess DSM error with global statistics on its accuracy (such as the commonly used root mean square error), but this does not include the spatial variation observed in the variogram.Measured ground or bare-earth LiDAR elevations can be considerably less accurate depending on the land cover type or terrain variables, such as slope, convexity and roughness [55].In addition, elevation errors are also spatially correlated, meaning that error values tend to be similar at nearby locations, and this can also have a considerable impact on the output of spatial analyses that use DSM as input.Monte Carlo simulations have been frequently criticized because of their high computing demand as the number of simulations has to be 500 at the very minimum [57].However, increased computing capacity in the last decade has facilitated its practical application in small-to medium-sized studies such as this one, in which 1000 simulations were performed.
To perform the geostatistical survey it is necessary to have a sufficient number of points (i.e., considerably more than 30) to be able to compute a robust experimental variogram [58].This can be a great drawback in undeveloped areas where the small-scale spatial structure of the errors may remain hidden due to sparse reference data.Obviously, the data quality of these control points also affects the error study.It is therefore essential to have well-surveyed data points with higher precision and a clear definition of data sources.Communicating uncertainty in terms of probability distributions to a diverse audience is always challenging, but it is becoming more popular and effective with the use of mapping techniques.In our case, analysis of the map of probabilities The advantage of using a geostatistical method is that it allows errors associated with the DSM values to be studied and the spatial structure to be taken into account by using the variogram.Also of benefit is the use of geostatistical simulations of equiprobable spatial distribution to determine uncertainty in the distribution of error values.It is possible to assess DSM error with global statistics on its accuracy (such as the commonly used root mean square error), but this does not include the spatial variation observed in the variogram.Measured ground or bare-earth LiDAR elevations can be considerably less accurate depending on the land cover type or terrain variables, such as slope, convexity and roughness [55].In addition, elevation errors are also spatially correlated, meaning that error values tend to be similar at nearby locations, and this can also have a considerable impact on the output of spatial analyses that use DSM as input.Monte Carlo simulations have been frequently criticized because of their high computing demand as the number of simulations has to be 500 at the very minimum [57].However, increased computing capacity in the last decade has facilitated its practical application in small-to medium-sized studies such as this one, in which 1000 simulations were performed.
To perform the geostatistical survey it is necessary to have a sufficient number of points (i.e., considerably more than 30) to be able to compute a robust experimental variogram [58].This can be a great drawback in undeveloped areas where the small-scale spatial structure of the errors may remain hidden due to sparse reference data.Obviously, the data quality of these control points also affects the error study.It is therefore essential to have well-surveyed data points with higher precision and a clear definition of data sources.Communicating uncertainty in terms of probability distributions to a diverse audience is always challenging, but it is becoming more popular and effective with the use of mapping techniques.In our case, analysis of the map of probabilities resulting from geostatistical research could lead to a series of decisions being taken.These might include carrying out general checks in those areas with high probability errors or doing additional sampling where the estimated error is greater than the threshold considered here (i.e., ˘0.5 m).
In the method applied here, analyses of other variables that can also be related to DSM errors are lacking.If spatial variation of the errors could be linked to terrain complexity, distance to streams, or vegetation cover, the origins of DSM errors could be better explained.If this were the case, estimation of the first floor elevation of buildings could be optimized, thereby generating more realistic values.Hence, one the challenges for further refinement of the technique would be proper use of auxiliary data.
A practical implementation of Monte Carlo simulation is the addition of a repeated number of realizations of the error model to the original data.The new set of modified DSM values can be used for further research, such as analysing the effects on derived topographic parameters, or as input in a hydrological model.Another future application of the present study could be to input the 1000 modified DSM values in the hydraulic model and analyze their influence on the accuracy of inundation mapping.

Conclusions
Determining first floor elevation is critical for characterizing flood risk.To our knowledge, this is the first time that the errors and uncertainty associated with this parameter have been estimated using LiDAR-derived DSMs as source data and Monte Carlo simulations based on geostatistics.The implementation of a probabilistic scheme for characterizing first floor elevation errors is very suitable, as it conveys the uncertainty inherent in the errors, which are spatially distributed.First floor elevation uncertainty is not only due to the global vertical accuracy of LiDAR data or errors of interpolation.A considerable role is also played by breaklines which, when integrated into TINs, represent complex terrain, such as urban areas, much better than TINs created with mass points only.As a counterpart, breaklines may to some extent reduce vertical accuracy of the original LiDAR as a result of the inherent simplifications to be assumed during the process of converting 2D breaklines data into a new feature class that contains z-values.Errors derived from the approach used here were similar to those described in other papers already published (i.e., in terms of measures of central tendency and spread), both in the urban area, where streets were integrated in the TIN as breaklines, and in the non-urban area where breaklines were not added, except on river banks.These errors and the associated uncertainty could be significantly improved by increasing LiDAR's point data.The denser the data, the greater the accuracy that can be achieved for determining first floor elevation.The approach deployed here is of paramount importance, particularly with regard to decision-making during the flood risk assessment and management process.This is because it not only enables flood damage to be assessed more reliably but also identifies the parts of the area prone to flooding that require improved topography, aspects that both contribute to a better characterization of hydrodynamic and economic losses.

Figure 1 .
Figure 1.Location of the study area.Coordinate system: ETRS89, UTM Zone 30 N. There are several points of conflict in the reaches studied: (i) B1 and B2 correspond to bridges on the Alberche River; (ii) W is a weir also on the Alberche River; and (iii) C1 to C7 represent culverts on the Alberche River and the Chorreron Stream.

Figure 1 .
Figure 1.Location of the study area.Coordinate system: ETRS89, UTM Zone 30 N. There are several points of conflict in the reaches studied: (i) B1 and B2 correspond to bridges on the Alberche River; (ii) W is a weir also on the Alberche River; and (iii) C1 to C7 represent culverts on the Alberche River and the Chorreron Stream.

Figure 2 .
Figure 2. Some of the Iber two-dimensional hydrodynamic software outputs.(A-C) show depths, velocities and Froude numbers in the study area, with the existence of hydraulic structures (see Figure 1) considered as a hypothetical scenario; (D-F) show mapping of the same parameters but without hydraulic structures in the model.

Figure 2 .
Figure 2. Some of the Iber two-dimensional hydrodynamic software outputs.(A-C) show depths, velocities and Froude numbers in the study area, with the existence of hydraulic structures (see Figure 1) considered as a hypothetical scenario; (D-F) show mapping of the same parameters but without hydraulic structures in the model.

Figure 3 .
Figure 3. Histograms of first floor elevation errors for both the whole study are and the 500-year flood zone.

Figure 3 .
Figure 3. Histograms of first floor elevation errors for both the whole study are and the 500-year flood zone.

Figure 4 .
Figure 4. Box-whisker plots of elevation errors for the first floor of buildings and the different land covers present at the study site.

Figure 4 .
Figure 4. Box-whisker plots of elevation errors for the first floor of buildings and the different land covers present at the study site.

Figure 5 .
Figure 5. Variogram model fitted for the DSM error at control points located at streets.Triangles represent the experimental omnidirectional variogram.The black line is a spherical model of sill of 0.08, range 122 and a nugget effect of 0.035.

Figure 6 .
Figure 6.Two of the 1000 elevation error realizations generated using Sequential Gaussian Simulation.Each realization is equally probable and has the same statistical properties as the measured errors.The area and location of the simulations is outlined in red on the orthophoto.A corresponds to simulation 114 and B to simulation 772.

Figure 5 . 17 Figure 5 .
Figure 5. Variogram model fitted for the DSM error at control points located at streets.Triangles represent the experimental omnidirectional variogram.The black line is a spherical model of sill of 0.08, range 122 and a nugget effect of 0.035.

Figure 6 .
Figure 6.Two of the 1000 elevation error realizations generated using Sequential Gaussian Simulation.Each realization is equally probable and has the same statistical properties as the measured errors.The area and location of the simulations is outlined in red on the orthophoto.A corresponds to simulation 114 and B to simulation 772.

Figure 6 .
Figure 6.Two of the 1000 elevation error realizations generated using Sequential Gaussian Simulation.Each realization is equally probable and has the same statistical properties as the measured errors.The area and location of the simulations is outlined in red on the orthophoto.A corresponds to simulation 114 and B to simulation 772.

Figure 7 .
Figure 7. Probability that first floor elevation errors are between −0.5 m and +0.5 m.

Figure 7 .
Figure 7. Probability that first floor elevation errors are between ´0.5 m and +0.5 m.

Figure 8 .
Figure 8.General view of the TIN of the study site.(A) shows a zoom-in of the same urban area including breaklines; (B) shows a zoom-in of the urban area without breaklines.

Figure 8 .
Figure 8.General view of the TIN of the study site.(A) shows a zoom-in of the same urban area including breaklines; (B) shows a zoom-in of the urban area without breaklines.

Figure 9 .
Figure 9. Classified land cover map with overlaid elevation errors for the 1987 field measurements provided by the Spanish cadastre.

Figure 9 .
Figure 9. Classified land cover map with overlaid elevation errors for the 1987 field measurements provided by the Spanish cadastre. ).

Table 1 .
Area occupied by deciles of probability considering that first floor elevation error is between ´0.5 m and 0.5 m.