Uncertainty and sensitivity analysis of coastal flood damage estimates in the west of the Netherlands

Uncertainty analyses of flood damage assessments generally require a large amount of model evaluations. This is often hampered by the high computational costs necessary to calculate flood extents and depths using 2-dimensional flow models. In this paper we developed a new approach to estimate flood inundation depths that can be incorporated in a Monte Carlo uncertainty analysis. This allows estimation of the uncertainty in flood damage estimates and the determination of which parameters contribute the most to this uncertainty. The approach is applied on three breach locations on the west coast of the Netherlands. In total, uncertainties in 12 input parameters were considered in this study, related to the storm surge, breach growth and the damage calculation. We show that the uncertainty in flood damage estimates is substantial, with the bounds of the 95 % confidence range being more than four times smaller or larger than the median. The most influential parameter is uncertainty in depthdamage curves, but five other parameters also contribute substantially. The contribution of uncertainty in parameters related to the damage calculation is about equal to the contribution of parameters related to the volume of the inflowing water. Given the emphasis of most risk assessments on the estimation of the hazard, this implies that the damage calculation aspect deserves more attention in flood risk research efforts. Given the large uncertainties found in this study, it is recommended to always perform multiple calculations in flood simulations and damage assessments to capture the full range of model outcomes.


Introduction
Flood management in Europe has traditionally mainly been concerned with measures designed to prevent flood events. More recently, a paradigm shift can be observed in flood risk management, moving towards an approach not only focussed on flood prevention (i.e. the hazard), but also on measures that reduce economical damage and casualties (DKKV, 2004;Tunstall et al., 2004;Vis et al., 2003). The interest in a risk-based approach has been amplified by severe flooding events in recent decades, as well as the observation that exposure to flooding has increased dramatically through population growth and urban development in vulnerable areas (Munich Re, 2005;Nicholls et al., 2008;. The prospect of future flood risk is continuing to increase because of these socio-economic drivers combined with the prospect of climate change (e.g. Te Linde et al., 2010;Milly et al., 2002) added further to this interest in risk-based flood management.
When focusing on flood consequences to support riskbased flood management, different models have been designed to assess direct economic damages to, for instance, buildings and infrastructure (Penning-Rowsell et al., 2010;Thieken et al., 2008;Kok et al., 2005). However, research shows there is quite some uncertainty in such flood damage simulations, which is rooted in the uncertainty of the underlying data and model assumptions Merz and Thieken, 2009;Apel et al., 2008;. These sources of uncertainty can relate to both epistemic uncertainty (incomplete knowledge), or aleatory uncertainty (natural variability) (see e.g. Apel et al., 2004). Research suggests to improve our understanding of the uncertainties in flood damage modelling, in order to make the results defensible and reliable (Saltelli et al., 1998), build public confidence in the output (Brugnach et al., 2007), communicate implications of limitations in data and scale (Hall and Solomatine, 2008), and ultimately provide the best possible basis for effective decision making (Ascough et al., 2008;Downton et al., 2005). Moreover, studying uncertainties forces researchers to scrutinize their own models , improving understanding of the system and the interpretation of the results.
Uncertainties can be addressed in several ways, for instance by using scenarios to create a range of possible outcomes. This is often done in studies concerning assessments of future situations which are inherently uncertain (Hall et al., 2005a;Aerts et al., 2008;Klijn et al., 2007;Bouwer et al., 2010), but can also be done to explore current uncertainties using, for instance, different inundation scenarios (Klijn et al., 2007;RWS-DWW, 2005; Van der Most and Wehrung, 2005). There are various studies that do not use scenarios but rather look at uncertainty in various components of a damage assessment. Some studies looked at uncertainties in a specific component, like the direct damage calculation part of the assessment Egorova et al., 2008) or the estimation of inundation depths (e.g. Hall et al., 2005b;Noack and Yoruk, 2008). More recently, studies have focused on the combined uncertainty in various components of flood damage assessments. Merz and Thieken (2009), for instance, used different flood frequency curves, inundation models and damage models to estimate uncertainty bounds around flood risk curves. A similar approach was used by De Moel and  who explored uncertainties in damage models, inundation depth and land use.
A more thorough analysis of uncertainty would require a statistical analysis of uncertainty in individual parameters, which are used as input for the various models and how they propagate through the entire modelling chain (Heuvelink, 1998). This is usually done using a Monte Carlo (MC) framework in which individual parameter values are (quasi-)randomly sampled and then evaluated in the model(s). This requires a large number (hundreds) of model evaluations. In flood damage assessments this is mainly hampered by the estimation of inundation depths as the models used for this (2-dimensional hydrodynamic flow models) are computationally very demanding (Gouldby and Kingston, 2007;Apel et al., 2008). The study by Apel et al. (2008) is a good example of applying MC to estimate uncertainty in flood risk modelling. To overcome the computational burden of 2-dimensional hydrodynamic models, they ran such models a priori for various dike breach scenarios to create lookup tables for inundation characteristics, which were subsequently used in a MC uncertainty analysis. They varied many different parameters in their study, though most were associated with the hazard part, and not so much the damage estimation. The damage estimation was performed after the MC analysis using three different depth-damage curves for buildings.
In addition to analysing the uncertainty in the output of a model or modelling chain, calculating a large number of model evaluations also makes it possible to perform sensitivity analyses. Sensitivity analyses enable distinguishing which uncertainty sources impact the output the most. Such insights allow identifying parameters, which are justified to keep constant, and can be used to direct research efforts to address parameters that heavily impact the output in order to make a more robust assessments. This paper aims to perform an uncertainty and sensitivity analysis of flood damage estimates, including uncertainty in input parameters associated with both the hazard part, as well as the damage calculation part of the assessment. Flood damage is quantified using a combination of three models: a breach growth model, an inundation model and a damage model. These models are embedded in a MC framework in order to determine the sensitivity of the model chain to different assumptions in the input parameters of these models, and to assess the uncertainty surrounding the resulting damage estimates. The following research objectives are identified: 1. To create a new and less computationally demanding modelling approach to estimate flood inundation depths after dike failure that can be incorporated in a MC analysis.
2. To identify uncertainty ranges for various input parameters of the different models.
3. To perform Monte Carlo based uncertainty and sensitivity analyses of flood damage estimates for three dike breach locations.
These analyses are performed for a case-study area with dikes and low lying polders in the western part of the Netherlands that is subject to storm surges. In Sect. 2 the case study area is introduced and in Sect. 3 the different data and models used to estimate flood damage are described. Section 4 then focuses on the uncertainty and sensitivity analyses, including the uncertainty ranges attributed to the input parameters. The results are described in Sect. 5 and conclusions are drawn in Sect. 6.

Case study area
Most of the flood-prone part of the Netherlands is divided into 53 so-called dike-ring areas (areas entirely surrounded by embankments or high ground) which have protection levels ranging from 1/1250 per year to 1/10 000 per year. Dikering area number 14, also called Central Holland, was selected for this study (Fig. 1). Dike-ring area 14 is the most densely populated and economically important part of the Netherlands. The three largest cities in the Netherlands (Amsterdam, Rotterdam, The Hague) are located in this dike ring, as well as the main airport (Schiphol Amsterdam). In total, dike-ring 14 covers 223 000 ha and contains 3.25 million inhabitants (RWS-DWW, 2005). The study area can be flooded from various sources. The western part of the dike-ring area, which is protected by dunes, can be flooded from the North Sea. The Rotterdam harbour area, just south of Dike-ring 14, is located behind a storm surge barrier (Maeslandt barrier), but flooding remains possible if the barrier fails. Though not considered in this study, the south eastern part of the dike-ring area can be flooded when a dike breach occurs in the dikes along the river Lek, a branch of the river Rhine.
Dike-ring area 14 consists of a large number of polders that can be regarded as smaller compartments. These polders are surrounded by smaller in-land levees, which are the by-product of a history of peat excavations and drainage of the resulting lakes. Some polders have elevations as low as 6.5 m below mean sea level. Recent flood simulations (e.g. Asselman, 2008) show that these smaller embankments act as barriers in case the main dike ring fails, and are very effective in slowing down the flooding process, and in limiting the flood extent (see e.g. Fig. 5 in Klijn et al., 2010). Because of its size, the topography and the smaller embankments, it is not expected that the entire dike-ring will flood after a flood defence fails.
For this study, three breach locations were chosen to estimate flood damages. These locations are Katwijk, Ter Heijde and Maassluis (Fig. 1). All three locations are vulnerable to flooding due to storm surges on the North Sea and correspond with places that are considered relatively weak points in the flood defence system (VNK, 2006) and are therefore currently the subject of dike-and dune-reinforcements.

Data and models
In order to assess potential flood damage, the following types of information or models are required: 1. Information on the hydraulic loadings (i.e. storm surge water levels) and their probability of occurrence (derived from hydraulic models or statistical analyses of time series).
2. Model to simulate breach growth after a dike fails.
3. Model to simulate the inundation resulting from the dike failure (hydraulic model).
4. Model to estimate the damage corresponding with the simulated inundation (damage model).
The methodology with the data and models used is schematically illustrated in Fig. 2 and the separate components are discussed in more detail in the following paragraphs. The entire modelling chain, including the uncertainty and sensitivity analysis, is coded in a set of Matlab functions and makes use of the existing Matlab toolboxes Simlab (Simlab, 2011) and the TopoToolbox (Schwanghart and Kuhn, 2010).

Hydraulic load and breach growth
In order to estimate the volume of water flowing into the dike-ring after dike failure, the hydraulic load, in this case the water level during a storm surge event, and size of the breach have to be estimated. This is done in the same way as in 2-dimensional flow models, using time dependant water levels and dynamic breach growth (Kamrath et al., 2006;Van Mierlo et al., 2007). The water level outside the embankments of the dike-ring is assumed to be a combination of the astrological tide (almost sinusoid), with a trapezoid storm surge superimposed (see Fig. 3), as is common practice in Dutch water management (Ministerie van V&W, 2007). The maximum increase in water level associated with the storm surge is associated with the return period on which the embankments are designed (for this dike-ring a water level with a probability of exceedance of 1/10 000 per year). For example, for the breach location Ter Heijde this corresponds with a water level of +5.7 m NAP (Dutch standard for mean sea level). This type of information is derived from the official documents of the Ministry of Transport, Public Works and Waterways in the Netherlands (e.g. Ministerie van V&W, 2007;RWS Waterdienst, 2008) and related websites (www.waterbase.nl and www.waternormalen.nl). The volume of inflowing water is also determined by the size of the breach and how fast it grows. For this, we use the Verheij-VanderKnaap formula (Verheij, 2002(Verheij, , 2003, which is also used in the 2-dimensional hydrodynamic flow model Sobek (see e.g. Van Mierlo et al., 2007). The growth of the breach depends on the difference in water level outside and inside the dike ring and the material of the dike itself. The water level difference over the dike is determined by looking up the water level just inside of the breach and subtracting that from the water level outside. With the width of the breach calculated for various time steps, the discharge flowing through the breach is calculated using the widely used Poleni's formula (see e.g. Kamrath et al., 2006). The total volume entering through the breach is then calculated by integrating the discharge through the breach over time (Fig. 3).

Inundation model
As a full 2-dimensional hydrodynamic model to simulate inundation after a dike failure is not so easy to incorporate in a Monte Carlo framework, because of the computational burden, a new approach to estimate inundation depths has been developed. This model is designed to directly calculate the end result of a regular 2-dimensional flow model, the maximum inundation depth of an area after flooding, given a specific volume (see Sect. 3.1). It has specifically been designed to simulate flooding in areas enclosed by dikes or other defence structures (bathtub flooding), as they are a distinct feature in the west of the Netherlands. Furthermore, it has been set up in such a way as to allow a large amount of model evaluations with differing volumes. In order to save computing time, we follow an approach similar to Apel et al. (2008), using a pre-process part that is performed before the Monte Carlo simulations, and a lookup part that is integrated in the Monte Carlo framework.

Pre-process: determining micro-basins and flooding sequence
In the pre-process step, the digital elevation model (DEM) of the area in question is split into numerous (thousands) individual depressions or micro-basins. The micro-basins are determined using functions from the TopoToolbox (Schwanghart and Kuhn, 2010). In this study the Algemeen Hoogtebe-stand Nederland data (AHN: general elevation database of the Netherlands) was used for elevation data, which has a spatial resolution of 5 m by 5 m and gives elevations in centimetres above NAP. Any "NoData" values (mainly water bodies) were filled using the minimum elevation of the cells surrounding it. Because of the size of the data set, the original 5 m by 5 m grid cells were aggregated to 15 m by 15 m grid cells, taking the mean of the aggregated cells. As this would result in a lowering of the elevation of the dikes surrounding the case study area, the cells corresponding to the primary dikes were aggregated using the maximum value of the cells during aggregation. Furthermore, as the AHN data is not filtered in urban areas, a technique was used that filters out cells based on the relative elevation difference, with adjacent cells to retrieve ground level elevations in urban areas.
In total, the dike-ring area was divided into almost 35 000 micro-basins, varying in size depending on the local topography (see Fig. 2 for an example). Generally the size varies between 100 m 2 (in areas with a lot of relief in the DEM, like urban centres) and 10 km 2 (in large flat polder areas).
After a breach location is defined, micro-basins are flooded one at a time starting with the one at the breach location. When a micro-basin is flooded, the lowest point around the inundated area is determined. This point gives the elevation up to which the micro-basin is filled and is also the point where water overflows into the next micro-basin. This next micro-basin is subsequently flooded in the same way. When during the process of flooding micro-basins the elevation of the overflow point is higher than the water level in previously flooded micro-basins, the water level in these previously flooded micro-basins is increased up to the elevation of the overflow point. The process of subsequently filling micro-basins continues until the flood water would flow out of the dike-ring area, representing the maximum inundation possible from the specified breach location. During this process, a lookup-table is created. After each micro-basin is flooded this table records the id-number of the basin, how high the water level is at that moment and the cumulative flood volume reached up to that point.

Lookup part: determining inundation map corresponding to volume
After the lookup-table has been created during the preprocess, it can be used to estimate the flooded area corresponding to specific volumes. The volume closest to the specified volume is looked up in the table, which provides the micro-basins that should be flooded, as well as their water levels. The water level in the last micro-basin is then slightly adjusted in order to reach the exact specified volume. While the pre-process of inundating the entire DEM to create the lookup-table is computationally demanding (numerous hours, depending on size and resolution), looking up the desired volume and creating the corresponding inundation map using the table is rapid (approximately one minute). In this way, many inundation maps, each corresponding with different inflow volumes of water (caused by differences in storm surge water level, flood duration, material of the dike, etc.), can be produced and used for subsequent damage calculations in a Monte Carlo framework.

Comparison with 2-D simulations
As no empirical data were available for model calibration, some of the inundation results were compared with inundation maps resulting from 2-dimensional hydrodynamic simulations carried out with Sobek. Sobek itself is also not calibrated for the case-study area under investigation due to the lack of empirical data. However, Sobek has been successfully used to simulate observed flood events in other areas (see e.g. Hesselink et al., 2003;Alkema, 2007) and is the principal model for flood inundation in the Netherlands. Figure 4 shows that our model produces flood extents that are more compact than those computed with the Sobek model. There are two main reasons for this. Firstly, the Sobek model also accounts for flow through local waterways (e.g. canals), through which water can propagate quickly. Areas located at a large distance from the breach location can be flooded by overflow of these canals. Secondly, our approach inundates one micro-basin at a time and completely fills it before looking for the next micro-basin. In Sobek, water can propagate in four directions at the same time.
Another reason for differences in flood extents simulated by our model and Sobek is related to the resolution of the DEM used. grid is used, whilst our model uses a 15 m by 15 m grid size. As a result, the elevation of line elements within the dikering area is more smoothed out in the 100 m grid, as opposed to the 15 m grid. To correct for this, the 100 m DEM used by Sobek is manually adjusted using auxiliary information on secondary embankments. In the simulations illustrated in Fig. 4 it can be seen that a small part of a secondary embankment was not (yet) adjusted (see the yellow ellipse in Fig. 4). As a result, water propagated through that area in the Sobek simulation, whilst in our model (using a 15 m grid) the water was stopped at the line element, resulting in different inundation extents.
Overall, there are some shortcomings associated with our inundation approach as opposed to full 2-dimensional hydraulic models (though the higher resolution of the DEM also brings some advantages). For the purpose of analyzing uncertainties and sensitivities, it is however more important to be able to perform many model calculations in a reasonable amount of time than simulating a single most realistic flood extent. For the objectives of this study the advantages of our newly developed approach thus outweigh the disadvantages.

Damage model
The damage model used in this study is the Damage Scanner, a meso-scale model which has originally been developed by De Bruijn (2006) for studying the effect of future land-use and climate change scenarios on flood damages (Klijn et al., 2007;Aerts et al., 2008). Since then, it has been used in various other studies that assessed the development of flood risk over time and uncertainties in flood damage assessments (e.g. Bouwer et al., 2010;. The Damage Scanner uses so-called depthdamage curves (Smith, 1994;Merz et al., 2007), which calculate potential flood damage by combining land-use information ("the assets at risk") and inundation depth. In the Damage Scanner, 13 land-use categories are distinguished, including three residential categories, infrastructure, industry and various agricultural land-use categories. For a more detailed description of the Damage Scanner the reader is referred to De Bruijn (2006), Bouwer et al. (2009) and De Moel and Aerts (2011).

Uncertainty and sensitivity analysis
Uncertainty analysis (UA) concerns the determination of the range in model output as a result of imprecisely known input parameters (see e.g. Helton, 1993;Crosetto et al., 2000). Sensitivity analysis (SA) is basically an extension of an uncertainty analysis, aiming to determine the contribution of uncertainty in individual input parameters on the range in output (Helton, 1993;Crosetto et al., 2000;Kann and Weyant, 2000). Performing sensitivity analysis enables to identify parameters that have little effect on the output and are thus justified to keep constant. In addition, parameters that have a large effect are also identified, which should thus receive additional attention in order to cope with the uncertainty they introduce. UA and SA are usually performed within a Monte Carlo modelling framework, as they both require a large amount of model evaluations. We follow the MC approach as described by Crosetto et al. (2000) and Helton (1993) who distinguish the following steps: (1) assigning distributions to input parameters, (2) generating samples of different combinations of input parameters, (3) evaluating the model using the generated combinations of input parameters, and (4) analyzing the results for uncertainty and sensitivity (Fig. 5).
In practice, there are two different approaches to SA; local and global Pappenberger et al., 2008). In the local approach, input parameters are varied one factor at a time while keeping the other ones fixed (Saltelli et al., 1999;Saltelli, 1999;Turanyi, 1990). The simplest implementation is to manually vary input parameters one at a time in order to judge the effect on the output. There are, however, also various techniques using random sampling and regression analysis (Turanyi, 1990). By varying parameters one at a time, only the sample space around the baseline situation is explored and interaction between input parameters is not accounted for. The above shortcomings are not present in global approaches. In a global SA, the entire sample space is explored and the total variation around the output is apportioned to the different input parameters, taking interaction between these parameters into account (Ratto et al., 2001;Archer et al., 1997;Saltelli et al., 1999;Pappenberger et al., 2006). In this study we therefore performed a global sensitivity analysis using Simlab functions for Matlab (Simlab, 2011) to sample different combinations of input parameters using the method of Sobol' (see e.g. Sobol, 2001). This method generates quasi-random samples that cover the entire sample space and gives first order and total sensitivity indices after the model evaluations have been processed. The first order indices indicate the relative importance of input parameters without any interaction, whereas the total sensitivity indices take all possible interactions into account.

Uncertainty ranges in input parameters
Distribution types and ranges for the different parameters considered in this study have been based on existing literature as much as possible. In cases where little to nothing is known on the shape of the distribution, a uniform distribution has been used. For most parameters some indication of the possible range could be derived from literature. For two parameters (Uc and Hfp) there was no literature readily available so the range was estimated by the authors using relevant data (the soil and elevation map, respectively). In the following paragraphs the distribution type and range will be discussed for each input parameter.

Water level
Uncertainty distributions were assigned to five parameters related to the storm surge water level (Table 1). These include the amplitude of the astrological tide (Amplitude), the elevation of the storm surge on top of the astrological tide (Max dH), the duration of the storm (DurationStorm), the duration of the peak of the storm (DurationPeak), and the timing of the peak of the storm surge with respect to the peak of the normal high tide (Offset).
To create a curve representing the astrological tide for the case study area, a time series was used for the year 2008. Using this time series, an average tidal curve was constructed and the natural variation in amplitude was determined (i.e. the Amplitude parameter). Analysis of the time series showed that variation of the amplitude of the tide can be characterized by a normal distribution with a standard deviation of 20 cm. The mean amplitude of the tide differs slightly per location (167-175 cm) and was derived from governmental databases. Superimposed on this astrological tide is a trapezoid shaped storm surge which is characterized by a total duration (DurationStorm), a duration of the peak which is somewhat flatter (DurationPeak), and a maximum height (Max dH; see Fig. 3). By default, a total storm surge duration of 35 h, with the peak lasting four hours, is used in the Netherlands for our case-study area (Ministerie van V&W, 2007). These numbers were therefore used as the means of normal distributions characterizing the uncertainty in these two parameters. The standard deviation related to the duration of the storm was retrieved from Asselman et al. (2010), who give an overview of various studies related to storm duration. While estimates of the standard deviation differ between studies (and locations), the standard deviation for the storm duration for the locations investigated in this study was about 10 h, and is therefore used in this study. The standard deviation for the duration of the peak of the storm surge was set arbitrarily to 1 h, roughly in line with the duration of the storm as a whole (about a quarter of the average duration). The maximum increase in water level associated with the storm surge was based on the design water levels as they are used in Dutch water management. In engineering practice, however, another 30 cm is added on top of this design water level to account for uncertainty in the calculation of this design water level and to be on the safe side (Ministerie van V&W and ENW, 2007). To account for this uncertainty, we therefore assumed a uniform (because of a lack of knowledge on the characteristics) distribution ranging from 30 cm below the design water level to 30 cm above the design water level in our calculations.
The timing of the peak of the storm surge with respect to the astrological high tide (Offset) finally determines the shape of the storm surge water level curve. The Ministry of Transport, Public Works and Water Management in the Netherlands assumes a phase difference between the two of 4.5 h, with the peak of the storm surge occurring 4.5 h after astrological high tide . This results in a prolonged period of high water levels after passage of the flood peak. This is probably done in order to be on the safe side, with respect to the consequences of a flood, as prolonged high water levels after a breach would enhance the inflow of water into the dike-ring area. There is, however, no reason to assume any link between the timing of the peak of the storm surge and the level of the tide. In our analysis we therefore set this timing completely random using a uniform distribution ranging from six hours before high tide to six hours after high tide.

Breach growth
For the breach formation part of the model there are also five parameters to which uncertainty distributions were assigned (Table 2). These parameters are the water level outside the embankment at which the breach starts to grow (StartLevel), the initial width of the breach that is needed as input for the breach algorithm (InitialWidth), a factor representing the erodibility of the material of which the embankment consists (Uc), the elevation of the land behind the embankment (Hfp), and the time it takes before the breach has reached this elevation (T0). Generally, a breach is assumed to occur when the water level outside the embankment reaches the design level. For the case of Ter Heijde this is, for instance, at +5.7 m NAP, corresponding to the 1/10 000 per year water level. Research has shown that the probability of failure of defences, especially for the three locations investigated in this study, can actually be lower (RWS-DWW, 2005). Based on failure probabilities found for various locations in the Netherlands, we assumed that breaching can occur starting from a water level associated with a probability 10 times higher than the design level (i.e. 1/1000 per year for our locations). This water level can be derived using so-called decimerings levels, which is a change in water level that is associated with a change in probability of 10 times. These decimerings levels have been calculated for all flood defences in the Netherlands (RWS Waterdienst, 2008) and amounts to 70 cm near Ter Heijde. Lacking knowledge on the characteristics of this uncertainty, a uniform range from +5.0 to +5.7 m NAP has therefore been assigned to the StartLevel parameter at Ter Heijde.
To model the width of the breach over time, the Verheij-VanderKnaap formula (Verheij, 2002(Verheij, , 2003 was used, which is based on empirical data. In this approach there is first vertical growth of the breach given a specified width (Initial-Width) until it reaches the lowest point behind the embankment (Hfp) in a specified amount of time (T0). After this initial vertical growth, horizontal breach growth occurs depending on the slope of the water level through the breach (which determines the flow velocity through the breach) and the material of the embankment (Uc). The default values used in the Sobek 2-dimensional hydrodynamic model for which the Verheij-VanderKnaap formula was developed are an initial breach width of 10 m, and a T0 of 0.1 h. Little is empirically known on the variation of these parameters, though Verheij (2003) quotes some maximum ranges that are quite large (1-100 m and 0.1-12 h, respectively). Given that the default values are quite low in this range, we have chosen to represent these two parameters with uniform distributions, ranging from 1-20 m for the InitialWidth parameter and 0.1-0.6 h for the T0 parameter. Values for the erodibility (Uc) of the embankment are usually ∼0.2 m s −1 for sand and ∼0.5 m s −1 for clay. The embankments at the two westernmost locations (Ter Heijde and Katwijk) are predominantly sandy. As little was known on the uncertainty of this parameter, we used a uniform range from 0.1-0.3 to represent Uc for these two locations. For the Maassluis breach location we assumed that the embankment is more clayey. This was based on the soil map of the Netherlands (Ten Cate et al., 1995), which shows that the soil in the vicinity of Maassluis is more clayey. As embankments are usually constructed from local material, it is likely that the embankment contains more clay as well. Here we used a uniform range from 0.25-0.45 (same range but higher values). The elevation and range of the land behind the breach (Hfp) was estimated from the digital elevation model for each location.

Damage estimation
For the damage estimation part of the modelling chain, uncertainty bands were designed for two parameters: the maximum value at risk and the shape of the depth-damage curve (Table 3). Uncertainty estimates for the maximum value at risk for the different land-use categories were mainly derived from Briene et al. (2002), supplemented with interpretations of Egorova (2004) and insights learned from De Moel and . Triangular distributions were used to describe the uncertainty in the maximum value at risk, following Egorova et al. (2008). To describe the uncertainty in the depth-damage curves, beta distributions (distributions between 0 and 1 based on two positive parameters) were used, also in accordance with Egorova et al. (2008). The exact shape of the beta distribution depends on the location on the depth-damage curve (i.e. the water depth) and a factor denoting the magnitude of the uncertainty in the damage curve ("k" in Egorova et al., 2008). Egorova et al. (2008)

Parameter
Distribution Range Value at Risk triangular depending on category (Briene et al., 2002) Depth-damage Curve beta depending on place in curve and "k-factor" (Egorova et al., 2008) advise to use a k-value smaller than 0.3 and use 0.1 themselves (the higher the k-value, the higher the uncertainty). We adopted the value of 0.1 as it was found that with such a value the highest damage factor (or alpha value) is two to three times that of the lowest damage factor (given a 90 % confidence range). This is in line with the observation of De Moel and Aerts (2011) who show that the damage curves used in the Rhine Atlas (ICPR, 2001) are about a factor two lower (less steep) compared to the depth-damage curves used in the Netherlands and Flanders (Kok et al., 2005;Klijn et al., 2007;Vanneuville et al., 2006). While the triangular and beta distributions are unique for each damage category, it was assumed that they are dependent on each other and the categories were not sampled individually but rather collectively. This means that for each Monte Carlo iteration, a p-value (between 0 and 1) was sampled for the value at risk and a p-value was sampled for the damage curves. These p-values were then applied on all values at risk and on all depth-damage curves for that Monte Carlo iteration using their respective distributions.

Results and discussion
For each breach location, 896 model evaluations were performed (this number depends on the Sobol' sampling scheme and the number of dynamic parameters). As a control, a calculation with 3584 evaluations was also performed for one location (Katwijk), yielding practically the same results. This implies that a total number of 896 model evaluations is sufficient to give statistically sound results. The results are shown in Table 4 and illustrated in Fig. 6. Table 4 shows mean and median damages for the breach locations and the 2.5 % and 97.5 % percentiles to illustrate the uncertainty. The standard deviation is not used to illustrate the uncertainty because of the non-Gaussian distribution of the output (see histograms in Fig. 6). This clear non-Gaussian distribution also means that the median is a better indicator for the average as compared to the arithmetic mean. Damage estimates using default values for certain input parameters as they are commonly used in the Netherlands (see Sects. 4.2.1 and 4.2.2) are also given for comparison. Between the three locations, a breach near Ter Heijde results in the largest direct flood damage (median of 3 billion Euro), closely followed by Katwijk (median of 2.6 billion Euro). A breach near Maassluis will result in substantially 3.5 × 10 9 4.1 × 10 9 0.44 × 10 9 Median economic damage ( C) 2.6 × 10 9 3.0 × 10 9 0.36 × 10 9 2.5 % percentile (C) 0.34 × 10 9 0.62 × 10 9 0.08 × 10 9 97.5 % percentile (C) 11.3 × 10 9 14.3 × 10 9 1.3 × 10 9 damage using "default" values ( C) 3.7 × 10 9 3.8 × 10 9 0.52 × 10 9 a Numbers between brackets correspond with the number of model evaluations that have been excluded as no breach occurred there.
lower damages (median of 0.36 billion Euro). The estimates using the default input parameter values are clearly higher than the medians of the distributions found in the Monte Carlo analyses, which shows that using these default parameter settings results, indeed, in a relatively high estimate of potential flood damage. VNK (2006) also calculated potential flood damages for dike-ring 14, using ten different breach locations. VNK (2006) calculated damages with and without the failing of line elements, resulting in two damage estimates per breach scenario. The estimates which assume no failing of line elements are conceptually the most comparable to our results, as in our study line elements in the DEM are assumed not to fail either. While the scenarios are not completely comparable, the results from VNK (2006) seem to be a bit higher than our estimates. Up to a certain degree this is not surprising, as VNK (2006)  . This is probably related to differences in the flood extent between this study and VNK (2006), as illustrated in Fig. 6. Results of the uncertainty (histograms) and sensitivity (pie charts) analyses for the three locations considered in this study. The histograms show the variation in the damage estimates resulting from the Monte Carlo analyses. Note that a different scale is used for the xand y-axis for Maassluis as compared to Ter Heijde and Katwijk. The pie charts represent the total variance of the output and indicate how much the uncertainty in each of the input parameter contributes to the variance of the damage estimates. Fig. 4 (VNK, 2006, also used Sobek, though not the exact same simulations as shown in Fig. 4). In our inundation estimate the flood waters accumulate more in low-lying polder areas that are mainly used for grazing. The Sobek calculation, on the other hand, floods a larger area near the coast including some more urban centres and horticulture areas, resulting in higher damages. The bar graphs of Fig. 6 and the 2.5 % and 97.5 % percentiles in Table 4 show the spread of flood damage estimates, and illustrate the uncertainty in such estimates. Overall, it can be observed that the uncertainty of flood damage estimates is substantial with boundaries of the 95 % confidence interval that are easily four times lower and higher than the median. The distributions are right-skewed and show a long tail going into high flood damage estimates. This seems to be a bit less for the Maassluis location, which is probably the result of the assumption that material of the flood defence is more clayey there, resulting in higher Uc values. These higher Uc values for Maassluis result in slower breach growth, resulting in a smaller range of inflowing water and thus a smaller range in resulting damage estimates.
The effect of the different Uc values is also apparent in the sensitivity analysis. Figure 6 (pie charts) shows that the influence of the Uc parameter on the total uncertainty is quite substantial for the locations Ter Heijde and Katwijk, but is much less important at the Maassluis location. Overall, there is a generally consistent picture across the three locations with respect to which input parameters contribute the most uncertainty. The order is not exactly the same everywhere, but in all locations the uncertainty in the shape of the depth-damage curves seems to be most influential. This is then followed by uncertainty in the erosivity of the material of the embankment (Uc), the duration of the storm (Dura-tionStorm), the timing of the storm surge peak (Offset), the maximum storm surge height (MaxdH) and the value of elements at risk (Value). Together, the variation in these six parameters causes over 90 % of the total uncertainty. When grouping these parameters into parameters related to the estimation of the inflowing volume (Uc, DurationStorm, MaxdH and Offset) and the damage calculation (Curves and Values), it seems that both groups account for about equal amounts of uncertainty around the final damage estimate. These findings fit well with the conclusions by Apel et al. (2008), who also find that uncertainties in the damage calculation for river floods are about equal to those in the Q-H relation (which determines the inflowing volume) for the lower Rhine region in Germany. Similar results were also found by Apel et al. (2009) in a different setting.
In risk assessments, substantially less effort often goes into estimating the potential damage as compared to the estimation of the flood hazard . Our results indicate, however, that uncertainty in the damage estimation is just as important as in the hydraulic boundary conditions. This implies that the damage estimation in flood risk assessments, but also in flood risk research in general, deserves more attention than it often currently gets. Nevertheless, a significant part of the uncertainty is related to the inflowing volume of water. This uncertainty does not only affect flood damage estimates, but also flood simulations in general. In practice though, many of the hydraulic parameters identified in this study as being influential are fixed in 2dimensional flood simulations, thus disregarding the uncertainty in the inflowing volume and resulting flood extent and depth information.
It is not straightforward to reduce the uncertainty in the parameters identified as being influential in this study. For the estimation of the inflowing volume, the value of Uc can be checked, up to some degree, by field surveys of the potential breach location. For the estimation of the duration of a storm surge event and the height of the storm surge this may be more difficult as they are derived from observed time series. These time series are limited in length (generally in the order of a 100 year record), meaning that large extrapolation is necessary because of the high safety levels in the Netherlands (1/10 000 per year in this case), which inevitably introduces substantial uncertainty (Te Linde et al., 2010). For the timing of the storm surge with respect to astronomical tide (offset), it is even theoretically impossible to reduce the uncertainty, as there is no reason to assume storm surges are somehow related to the lunar cycle. With respect to the damage estimation part, there is uncertainty in the generalization in damage categories (in our case land-use classes), but there is also a stochastic component (also known as aleatory variability) in that one does not know a priori which assets will be hit by floating debris (and sustain more damage) and which will not be hit. This latter uncertainty is impossible to minimize but damage estimations can theoretically be improved by using more detailed information on the assets at risk for the location in question. Subdividing residential land-use classes further to more detailed categories and/or using information on the state of individual buildings would allow us to define more detailed depth-damage curves and to better differentiate values at risk. A prerequisite is of course that such detailed information is available from, for instance, insurance companies (for values at risk) or national surveys on building quality. For industrial land use a further distinction into different industrial types (chemical industry, oil refineries, manufacturing, etc.) would also be worthwhile, but may be more difficult to assess than the residential assets because of the wide variety of industries and their susceptibility (Seifert et al., 2010). Bottom-up approaches, such as the reference installation approach developed by Geldermann et al. (2008), seem most appropriate to capture this heterogeneity if the available data allows it.
Finally, it has to be noted that these findings concern coastal storm surge flooding in a case-study area characterized by low-lying polder areas, partly below mean sea level. In other geographic locations the uncertainty and influence of modelling parameters may be different.

Conclusions
In this study, uncertainty and sensitivity analyses have been performed on flood damage estimates related to coastal storm surges in the west of the Netherlands. To facilitate this, a new approach to calculate inundation depths that can be incorporated in Monte Carlo simulations has been developed and successfully applied. Similar to the rapid inundation model developed by Ward et al. (2011), the purpose of this new modelling approach is to complement -rather than replace -existing 2-dimensional hydrodynamic approaches to answer questions that require high amounts of model simulations. The new computational efficient approach estimates inundation depths satisfactory and allows performing uncertainty and sensitivity analyses using Monte Carlo simulations. The combination of uncertainty and sensitivity analyses gives more insight into which parameters are most crucial for the uncertainty in the output of the modelling chain and thus allows for tailored follow-up activities.
We show that considerable uncertainty is associated with flood damage estimates related to coastal storm surges in the west of the Netherlands. The upper and lower estimates of the 95 % confidence range are easily four times smaller or larger than the median. The most influential parameter contributing to this uncertainty found in this study was uncertainty in the shape of depth-damage curves. Overall, the contribution of uncertainty in parameters related to the damage calculation (depth-damage curves and values at risk) is about equal to the contribution of parameters related to the volume of the inflowing water (storm surge duration, height and timing of the storm surge and the material of the embankment). This shows that attention in flood risk assessments and flood risk research efforts should be more balanced between the two, as it is currently often skewed towards the hazard assessment.
To increase the accuracy of flood damage assessments, the uncertainties in the parameters identified as being influential should be reduced. This is, however, not always easy to accomplish and in some cases even impossible. Assuming that it is unlikely that many of the uncertainties identified in this study will reduce considerably in the near future, it is important to still account for them in flood damage assessments. It is therefore recommended to make a range of estimates (e.g. high-middle-low) by adjusting some or all of the parameters that significantly influence the damage estimate. Given that a substantial amount of the uncertainty relates to the estimation of the inflowing volume, this does not only hold for damage assessments, but also flood simulations in general. These ranges should then transparently be communicated to follow-up activities (like cost-benefit analyses) and decision makers.