Probabilistic Flood Hazard Maps from Monte Carlo Derived Peak Flow Values—An Application to Flood Risk Management in Zamora City (Spain)

Featured Application: The present manuscript shows an easy and achievable approach to improving ﬂood hazard maps by delimiting the uncertainty derived from ﬂood frequency analysis, and how this propagates through to the hydrodynamic modelling results on which these mappings are based. Abstract: All ﬂood hazard and risk assessment suffer from a certain degree of uncertainty due to multiple factors, such as ﬂood frequency analysis, hydrodynamic model calibration, or ﬂood damage (magnitude–damage functions) models. The uncertainty linked to the ﬂood frequency analysis is one of the most important factors (previous and present estimation point to 40%). Flood frequency analysis uncertainty has been approached from different points of view, such as the application of complex statistical models, the regionalization processes of peak ﬂows, or the inclusion of non-systematic data. Here, we present an achievable approach to deﬁning the uncertainty linked to ﬂood frequency analysis by using the Monte Carlo method. Using the city of Zamora as the study site, the uncertainty is delimited by conﬁdence intervals of a peak ﬂow quantile of a 500-year return period. Probabilistic maps are derived from hydrodynamic results, and further analysis include ﬂood hazard maps for human loss of stability and vehicle damage. Although the effect of this uncertainty is conditioned by the shape of the terrain, the results obtained may allow managers to achieve more consistent land-use planning. All those Zamora city results point out the probable underestimation of ﬂood hazard (the higher hazard areas increase around 20%) and risk when the uncertainty analysis is not considered, thus limiting the efﬁciency of ﬂood risk management tasks.


Introduction
Floods are probably the most frequently recurring natural phenomenon affecting society (human and goods) in terms of space and time, regardless of their geographical location or socioeconomic development as shown by the data collected by the International Disasters Database for the period 1900-2018 [1]. This is the main reason why flood risk management has become an essential tool from both an economic and a social perspective, with the main objective of reducing their associated losses. Analyses of economic damage are always easier than those of social damage, due to the not-easily-valuable prize of human life. Assessing direct tangible damages (those direct damages resulting from the physical contact of floodwater with property and its contents, and the most common ones), economic flood losses have been increasing throughout the past half-century. In the last decade (2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018), economic losses associated with floods exceeded USD 35 billion [1] and, as an example within this period, flood losses exceeded USD 19 billion in 2012 alone [2][3][4].
Clearly, any assessment of such flood risks must commence with an evaluation of the magnitude of peak flow events with particular return periods. Management of flood risk then proceeds by setting some design recurrence interval that drainage network and urban water systems must be able to deal with. However, an evaluation of peak flow values is not an easy task, due to the uncertainty associated with the limited data availability from gauge stations. Thus, even if there are data series of maximum annual flow with lengths greater than 30 records, a value considered sufficiently long by Aronica et al. [5] for the application of the cumulative probability distribution function (CDF) to the annual maxima series, the representativeness of this time series with respect to the global behaviour of the peak flow variable is limited. This lack of knowledge causes the existence of the so-called epistemic uncertainty, resulting from imperfect knowledge of the system (e.g., [6][7][8]), or from simplifications associated with the selected modelling approach and parametrizations.
Due to this limited knowledge, the use of a probabilistic approach to flood hazard mapping should be preferred to deterministic ones for at least three main reasons [9][10][11]: (1) hydrological and hydraulic assessments are always affected by uncertainty; (2) a correct representation of the results of a hydraulic analysis should also quantify its associated uncertainty, and this goal only can be accomplished within a probabilistic framework; and (3) the end users (stakeholders and decision makers) can only guide and support the correct definition of flood mitigation strategies under the best knowledge of flood hazard and risk.
These probabilistic maps try to delimitate the uncertainty in flood hazard analysis, which could be related to many variables. Many previous studies analysed the effect of uncertainty associated with the roughness surface (Manning's n value) of hydraulic models (e.g., [5,[12][13][14][15]). Additionally, other studies (e.g., [5,[15][16][17][18]) analysed the uncertainty in upstream and/or downstream boundary conditions. Other authors considered additional uncertainties in flood hazard by including extreme value statistics [19][20][21], in many cases, by relating the changes in rainfall or flow results that may occur as a result of climate change (e.g., [22][23][24][25][26]). The importance of an uncertainty analysis in flood risk management was pointed out previously [7,27] in order to understand the implications for decision makers of limited data availability, model uncertainties, or other uncertainty sources. Monte Carlo simulations (MC) as a technique to reduce uncertainty in flood hazard analyses have been widely used. Previous studies have used MC methods (e.g., [15,17,19,[27][28][29]) to analyse all variables: surface roughness, boundary conditions, and return period as flow frequency variable.
Several studies have already analysed uncertainty in the flood hazard and risk analysis, how this uncertainty affects the results obtained, and consequently, the effect it has on flood risk management. So, the goal of this study is not to reinforce these past findings, but to try to quantify the effect of the coupled peak flow return period data in urban areas, especially when multiple cultural heritage elements are located in them. In order to achieve these objectives, the information recorded at two gauging stations was used to characterise (in a first phase) the peak return period ratio for the Douro River in Zamora. In a second phase, Monte Carlo techniques were applied to analyse the uncertainty of the peak flow/return period ratio, caused by the limited length of the recorded data series and, therefore, a partial knowledge of the general behaviour of this ratio. From the Monte Carlo method results, the peak flow return period ratio and its associated uncertainty were delimited through the estimation of the confidence intervals of this ratio. Using a 2D hydrodynamic model, flood hazard scenarios were then analysed for different elements exposed in the city of Zamora. All those results allow for a more consistent representation of flood-prone areas, as well as for hazard level classification. These flood hazard scenarios offer us, and may offer flood risk managers in the future, the possibility of an estimation of direct economic damage and its related uncertainty. On the other hand, changes in the distribution of pedestrian hazard stability for flooded urban areas [5,30] were also analysed as have those related to car stability, with the results showing a significant increase (up to 20%) of high level hazard areas from the benchmark scenario (those represented by the results associated with the estimated peak flows, without considering its uncertainty as confidence intervals). Finally, all these results were integrated and mapped to analyse changes in flood risk management associated with uncertainty in the estimation of flood hazard.

Study Area
The study area comprises the Douro River reach in Zamora city (Spain), which is 7 km long with a mean width of 225 m and a riverbed slope of 0.0006 mm −1 ; it is characterized by a relatively deep, clayey sand channel with some fully vegetated sandy bars and moderate sinuosity. Each of the riverbanks show a different shape (Figure 1), so the left bank shows a flat floodplain while the right bank shows a steep slope, where the urban centre of Zamora city is mainly located, although some neighbourhoods are located on the left bank, which is more prone to flooding as the historical record shows. The city of Zamora, which was declared a Historic-Artistic site in 1973, has numerous cultural heritage sites, and although these are mainly located on the right bank of the Douro River, some are located on its left bank. The urban nature of the study area means that there is a large population (the city of Zamora has 61,500 inhabitants), which may suffer from the effect of flooding. In addition, the presence of many cultural assets increases the importance of potential flood losses, due to their non-evaluable value, given the consideration of the impossibility of restoration or recovery of the possible damages caused by natural hazards.

Materials and Methods
Daily peak flow records and surface topography are the two main data sources ( Figure 2) used in the flood hazard and risk analysis in the city of Zamora, Spain. The Douro River flow data are recorded in two flow gauge stations near the study area, one located just upstream of Zamora city (gs2121) and another one about 37 km upstream in the city of Toro (gs2062). Between the two gauging stations, there is no large river mouth. The daily peak flow record of gs2121 is only 14 years long, while the gs2062 has a 44-year-long record. So, with the aim of improving the gs2121 peak flow record, annual maximum peak flow values of gs2121 and gs2062 were plotted, as was a linear regression fit with a R 2 value of 0.9807 (Figure 3a). Linear regression Equation (1) was used to improve and extend the gs2121 maximum annual peak flow record.
The linear regression equation shows the existence of a flood wave lamination effect between both stations, with a net loss of hydrograph peak flow between them.
On the other hand, the study area surface topography is publicly available as point cloud Light Detection and Ranging (LiDAR) data [31], which are filtered and classified following the American Society for Photogrammetry and Remote Sensing (ASPRS) standard classification. The vertical accuracy of the LiDAR data used is reported to be 20 cm [31].
The Terrain and Building LiDAR class types were used to derive the Digital Elevation Model (DEM). Due to the absence of LiDAR data in the river channel (since water has a low reflectance), it was modelled from bathymetry data obtained in the form of crosssection measurements from boat-mounted echo sounders with D-GPS (Differential Global Positioning System) data acquisition. Both LiDAR and bathymetry topographic data were combined into a 1 m spatial resolution DEM. The use of Flood Frequency Analysis (FFA) is common to estimate flood design values by fitting a probability distribution to the historical (recorded) flood data series [32]. Commonly used probability distributions for FFA include the Gumbel, lognormal, log-Pearson type III, and the generalized extreme value (GEV) distributions; precisely the latter is the one proposed as the best option for the study area [33]. Distribution fitting is probably usually highly dependent on hydrologic data availability. So, longer data series can improve the hydrologic frequency calculation. However, in practice, the lengths of measured hydrologic data are generally not long enough, and probably not fully representative of the inherent variability of peak flows over time. In fact, the presence of wet and dry periods can affect the maximum annual peak flow records, and this leads to uncertainty in the design and management of water resources. This uncertainty is usually illustrated by a confidence interval (e.g., [32,34,35]). The computing of confidence intervals for design floods from probabilistic models was previously investigated (e.g., [21,34,[36][37][38]) by different methods, such as the Delta, profile likelihood function (PLF), and the Bayesian Markov chain Monte Carlo (MCMC) methods. In the present study, an easier approach to computing of confidence intervals is used as is shown at Figure 2.
During the first step, an FFA for the extended gs2121 maximum annual peak flow record is carried out with AFINS software, which uses the GEV probability distribution with a Maximum Likelihood (ML) estimator for parameter (location-x 0 , scale-α, and shapeβ) adjusting. During the second step, the Monte Carlo method is used for the uncertainty analysis by computing confidence intervals. For that purpose, 500 synthetic series of 500 (uniformly distributed) probability values were generated in the MatLab software environment. The use of Uniform PDF is based on the fact that this distribution does not assume any prior knowledge of parameter distribution, providing an equal probability for any possible parameter value, so it is appropriate in the absence of verifiable data [14,39,40].
Each probability value of each series was transformed into a peak flow value by using the GEV-ML parameters previously obtained, so 500 annual maximum peak flow synthetic series of 500 record data were available for FFA by using GEV probability distribution. As a result of the second FFA, five hundred values of peak flow values related to each return period considered were available for classical inference statistics, which were used to compute the mean and confidence intervals for the range of return periods defined (2, 5, 10, 25, 50, 100, 200, 500, and 1000 years). For computing the confidence intervals of each T year return period peak flow value, the variance of GEV distribution was obtained following the approach of Lu and Stedinger [41].
The computed confidence intervals are considered as an uncertainty measure of the T year return period peak flow values. Flood frequency analysis is one among a wide range of factors that can lead to uncertainty in flood risk assessments. This range of factors could be divided into two main groups or types [6,8]: (a) the random uncertainty linked to the natural variability of natural process (which cannot be reduced by more detailed information), and (b) the epistemic uncertainty related to the partial or limited knowledge of the natural process (which can be reduced by acquiring more knowledge in the form of more data or more sophisticated analysis). This epistemic uncertainty is taken into account by a probabilistic approach to flood maps (against the deterministic approach, which usually does not take into account the uncertainties [13]), such as those produced in the Zamora city flood hazard assessment presented here.
So, from the data obtained in the previous steps, we proceeded to the 2D hydrodynamic modelling of the peak flow related to the 500-year return period, as well as some of the values associated with the different confidence intervals considered (50%, 80%, 90%, and 99%). For hydrodynamic computing, Iber software [42] was used. As a result of this hydrodynamic modelling, different flow depth and flow velocity maps were obtained for each of the modelled peak flows.
Finally, the flood hazard and flood risk assessment was carried out in the location of Zamora city. For the flood hazard analysis (mainly related to the possible human stability lost due to flow parameters), two approaches were used to categorize hazard levels, one proposed by Aronica et al. [5], which is in the style (although a little more elaborate) of that proposed by the Spanish legislation [43], and the other one proposed by Russo et al. [30]. The proposal of Aronica et al. [5] categorizes the flood hazard into four levels, while the proposal of Russo et al. [30] uses five levels of flood hazard, from "Low hazard for adults and children" when the "flow depth * flow velocity" does not exceed the value of 0.4 m 2 s −1 to "Extreme hazard for adults" when it is over 1.2 m 2 s −1 , the flow depth is over 1.2 m, or the flow velocity is over 3 m s −1 .
As specific flood damage models are not available for Zamora city, general purpose United States Army Corps of Engineers (USACE) [44] and Joint Research Centre from the European Commission [45] flood damage models were selected for flood risk analysis. Those flood damage models were applied for the hydrodynamic results related to the 500-year return period peak flow and its confidence intervals. Monetary economic damage values, differences, and uncertainties were used in aiming to improve flood risk management in Zamora city. A damage model for vehicles [46] was also considered in the manner of completing the flood hazard assessment in the urban area of Zamora city; in this case, the most common vehicle types (Sedan and SUV) in Zamora, from the USACE model vehicle options, were considered.

Results
Flood hazard probabilistic maps for Zamora city were derived from the MC analysis. The first time, the improved and extended maximum annual peak flow record for gs2121 gauge station at Zamora city was used for FFA by using GEV probability distribution with a Maximum Likelihood (ML) estimator for parameters (Figure 3b). This FFA gave us the GEV parameter values, which were used further in the MC analysis. The values of the GEV parameters were then used to transform the 500 series of probability values (according to a uniform distribution) into series of instantaneous maximum annual flow values. For each of these 500 series of peak flows, a GEV distribution function was fitted again to obtain the maximum instantaneous peak flow values associated with the different return periods considered in the study, as well as the confidence intervals of this estimated value. The results of the MC analysis for Zamora city are shown in Figure 3c and Table 1. Using the 500-year return period values as an example for flood hazard analysis, 2D hydrodynamic models were computed at the study site for the quantile value and some of the upper side confidence intervals (50%, 80%, 90%, and 99%). From hydrodynamic computing some flooding variables were defined as the flood-prone area, flow depth, and flow velocity. The results about flood-prone areas allow us to define the extension and probability of flooding areas in Zamora city, and at the same time, to identify the urban areas where, with a probabilistic confidence of 99%, there is an absence of flood hazards related to the 500-year return period. Figure 4 shows the probability map of flood-prone areas, where the left bank of the Douro River denotes a higher uncertainty due to its flatter shape. A clear example of this uncertainty is observed at the urban development of Zamora city on the left bank of the Douro River, which is partially affected by flooding when the quantile peak flow value is used but is fully affected when considering the 99% confidence interval. Flow depth and velocity are the main flood variables used to define hazards due to human loss of stability and direct economic damages. For the first one, the models proposed by Aronica et al. [5], and Russo et al. [30] were computed for Zamora city, using both the quantile and confidence interval values. For the 500-year return period flood, both models show ( Figure 5) a predominance of the highest hazard level, from the results associated with the hydrodynamic modelling of the quantile value associated with that return period (extending this maximum hazard zone as we consider the different confidence intervals). Differences between the two models considered can be observed, so the Aronica et al. [5] model shows a greater extension of mid-hazard levels than the Russo et al. [30] model. These differences may be due to the fact that the Aronica and co-authors model takes into consideration both variables, the flow depth and velocity, to define hazard levels. The model proposed by Russo and co-authors does so as well, but for the two higher levels, they use a Boolean selection, which is based on two or three options: (1) a value of "flow depth * flow velocity"; (2) a value of "flow depth"; or (3) a value of "flow velocity". Both models of human loss of stability show minor differences between the quantile and confidence interval hydrodynamic results for most of the urban area of Zamora city, mainly on the right bank of the Douro River (Figure 5a,e). More differences are appreciated on the left bank of the Douro River, where mid-and low-hazard levels dominate for the 500-year quantile model, while the high-hazard level increases its extension for the different confidence interval models (Figure 5b-d,f-h). This trend is homogeneous for both the Aronica et al. and the Russo et al. model.
The direct economic damages for housing related to the 500-year return period quantile and confidence intervals were computed by using the USACE [44] and JRC [45] models. The first observation that can be made of the results is associated with the maximum level of damage reached in each model. This factor is conditioned by the different maximum percentage of damage considered in each model: 71% in the USACE model versus 100% in the JRC model. However, there is a certain parallelism in the distribution of areas for both models ( Figure 6).
As in the previous analysis (loss of human stability), the greatest differences are seen in the part of the Zamora city located on the left bank of the Douro River. In this area, there is both an extension of the affected area and a significant increase in the percentage of damage to buildings.
The damage model proposed by the USACE assigns percentages of damage (for quantile peak flow value, Figure 6a) that fluctuate between 10% and 50% for the urban area of Zamora. These percentages are concentrated between 25 and 50% when we consider the value of the confidence interval of up to 80% (Figure 6b,c). For the confidence interval of 99%, the area with damage percentages above 50% becomes dominant, preferably for the urban area located on the left bank of the Duero River (Figure 6d).  The damage model proposed by the JRC generally assigns higher damage percentages. In the urban area of Zamora, both for the 500-year-old quantile and for its confidence intervals of up to 80% (Figure 6e-g), the damage percentages fluctuate between 25 and 75%, with the average value increasing in direct relation to the confidence interval. For the 99% confidence interval (Figure 6h), the area with expected damage between 50 and 75% is dominant, and the area with damage greater than 75% becomes noticeable.
Like people and houses, vehicles are another property type exposed to flooding, which are affected in two ways: firstly, by direct damage due to flooding, and secondly, by the increased flood vulnerability to other exposed property when the vehicle is moved by the flood as a suspended load. To partially analyse this situation, the vehicle damage model [46] proposed by the USACE was used; specifically, the models associated with sedan-and SUV-type vehicles (as they are probably the most common vehicles in Zamora city). The hydrodynamic results for most of the urban area of Zamora city show high percentages of vehicle damage (Figure 7). These damage percentages are higher for sedantype vehicles than for SUV-type vehicles since the USACE damage model considers (for flow depth values of up to 2 metres) the damage to sedan vehicles to be between 10 and 15% higher than the damage to SUV vehicles. The urban area of Zamora city located on the left bank of Douro River is again the area that shows the greatest variation between the results associated with the 500-year return period quantile and those of its confidence intervals. For sedan-type vehicles, there is a clear increase in the percentage of damage (towards values above 75%) in the areas common to all models (Figure 7a-d), while the areas affected only by the models associated with the confidence intervals present damage percentage values of over 50% (Figure 7d). In the case of SUV-type vehicles, the trend is similar to that noted above, with the difference that for the model associated with the quantile peak flow, there is a large area with damage percentages of between 25 and 50% (Figure 7e). This area is reduced for the 50% confidence interval (Figure 7f) and almost disappears for that of 80% (Figure 7g). Finally, the model associated with the 99% confidence interval shows damage percentages higher than 75% ( Figure 7h); in the uncommon areas with the model associated with the cliff, damage percentages lie between 25 and 50%. All these high values of damage to vehicles are very likely to result in vehicles being washed away, thus increasing the vulnerability of the elements located downstream.
If we focus the results on the consolidated urban area of Zamora city (both sides of the Douro river, 5.6 km 2 in extension), the flood-prone area increases from 0.27 km 2 (500-year return period quantile) to 0.31-0.34-0.36-0.38 km 2 (for the confidence intervals of 50-80-90 and 99%). These results point out an increase of 41% in the flood-prone area, related to the peak flow quantile and 99% confidence interval. All these changes in the flood zone are also reflected in the distribution of the different flood hazard levels associated with the human loss of stability, housing economic losses, and damage to vehicles, as can be seen in Table 2. In all cases, two trends can be identified: the first trend shows a downward trend (from the quantile to 99% confidence interval) in the area (%) occupied by one hazard level, but the second trend shows an upward tendency only for the higher hazard level of the damage model. Those results show us that the increase in the flood-prone area from the quantile to 99% confidence interval is not evenly distributed, but it promotes the spatial expansion of areas with a higher flood hazard level. Clear examples of this can be found in the results of the two damage models associated with building losses, where a higher flood hazard level rises from 14% to 36% (USACE model) and from 9% to 27% (JRC model).
In general, the increase in the higher flood hazard level area (%) could be estimated at around 20%, that is, half of the total increase of the flooding area is related to the higher flood hazard level. The two flood damage models related to human loss of stability tend to concentrate the results in the lowest and highest hazard levels, with a very low representativeness of the intermediate hazard levels. This effect could be either due to the results of the flow depth and flow velocity variables specific to the study area, or a more general effect linked to the models themselves. This trend, however, is not observed in the results linked to the economic damage models for housing and vehicles.
The use of different damage models developed for the same variable (human loss of stability, housing economic losses) highlights the uncertainty in the results associated with the damage model selection. This is clearly perceptible in the case of the city of Zamora when analysing the results associated with economic damage to buildings. Thus, the model proposed by the USACE [44] does not give results for damages higher than 75% (since the model itself does not consider losses higher than this value, regardless of the value of flow depth), and the majority of losses are between 25 and 50%. However, the model developed for the JRC by Huizinga et al. [45] does consider losses of up to 100%, and this is clearly reflected in the results, where the most frequent loss range is 50-75%.

Discussion
Conceptually, all flood hazard and risk maps are probabilistic, although classical maps are sometimes referred to as "Deterministic Flood Maps" [47]. As these are not representations of tangible objects or strictly deterministic variables (as they are in a topographical, geological or vegetation map), the elements and immaterial variables represented (hazard or risk) are of a potential future nature and, therefore, stochastic. This stochastic component can be derived both from the hydro-meteorological methods of estimating flood discharges and from the flood frequency analysis, which use probability distribution functions and parameters such as the return period or the probability of exceedance. However, it can also come from the spatial-temporal probability of the location of exposed elements (people, material goods, commercial and service flows) and their vulnerability, which do not remain constant in time and space, and therefore must also be assigned a probability of occurrence.
However, "Probabilistic Flood Maps" (PFM; [11]) are those which, in addition to the intrinsic stochastic component of hazard and risk, have other probabilistic elements incorporated into the cartographic representation, such as the confidence intervals of the frequency distribution function of the values represented. Therefore, probabilistic mapping of hazard and flood risk is characterized by a zoning of the values (in discrete classes or continuous models), not only of the probability of the phenomenon or of the exposure and vulnerability of the elements at risk, but also of the uncertainty or error in the estimation of the variable represented. This uncertainty can be methodological, epistemic or associated with data sources and boundary conditions [10].
Probabilistic cartography has a number of advantages over cartography that could be called classic or conventional. The main advantage is its versatility, because it makes it possible to incorporate a wide range of errors and uncertainties onto the probabilistic intervals in the estimation of the variables involved in the cartographic production. This means that it is easy to adapt the cartographic result to future changes or hypotheses, such as climate change scenarios, geomorphological evolution or land-use evolution models [15]. Therefore, probabilistic cartography is ideal for future forecast maps, such as the "Flood Factor" maps of the First Street Foundation [48], which generate risk cartography for past, present and future time horizons (within 30 years). This versatility allows the incorporation of novel methodologies based on decision making under uncertainty theories, in particular "Value of Information" (VOI) as proposed by Alfonso et al. [11].
From an applied point of view, this advantage of probability mapping allows this versatility to be transferred to flood risk management and the adoption of mitigation measures and strategies. For example, in the field of insurance and reinsurance, probabilistic cartographies reach their maximum development, as they allow the correlation of risk in coverage with the probability of risk of loss [49], and the adjustment of insurance policies with confidence intervals, following more conservative or more entrepreneurial policies, replacing the classic FEMA "Flood Insurance Rate Map" (FIRM) [50]. Another field of application where probability maps could be very useful is in territorial and urban planning since the different land classifications (urban, rural, and protected) and land-use restrictions (residential, industrial, recreational, and services) could be associated not so much with the probability of the phenomenon occurring as with the uncertainty of its estimation. This would avoid the so-called "border or edge effects" in landscape and urban planning, associated with abrupt zoning limits, which are so damaging to risk management, and would encourage a greater progressive gradation in classes, types of uses and use restrictions in the territory [51][52][53]. Finally, a flood risk analysis for critical infrastructures, such as nuclear power plants, large reservoir dams, industrial poles and the like, requires cartographies that make it possible to remain on the safety side, using the maximum magnitude and minimum probability confidence intervals when selecting their location, design and dimensioning [54,55].
However, not all are advantages in the use of probability mapping over conventional flood hazard and risk mapping. In the legislative field, and specifically in administrative and criminal law, the diffuse limits associated with statistical probabilities generate problems of interpretation and application [56,57]. If technical and legal problems already exist when legislation establishes the zoning of activities based on return periods-as in Spain, where the floodable zone is defined as that flooded by the 500-year return period [58]-the legal difficulties are extended by cartographies that add uncertainties and confidence intervals, which are difficult to incorporate into legislation and regulatory development.
Despite these legal and other methodological limitations, the different comparative analyses between deterministic and probabilistic maps for case studies [47] have demonstrated the usefulness of the latter maps, especially if they are used as an element of relative comparison, as is the case in this work. The probabilistic cartography produced for Zamora city can also be used as a peak flow uncertainty measure. So, for the 500-year return period, this uncertainty is about 36.5% when you compare the quantile and the 99% confidence interval peak flow values. This uncertainty (although the approach to their magnitude is different in both studies) is similar in magnitude to that obtained in the town of Navaluenga (Spain) by Ballesteros-Cánovas et al. [59], which was just above 40%. This uncertainty magnitude, associated with the flow frequency law, was found to be much higher than other sources of uncertainty, such as the damage model, or the hydrodynamic model's calibration [59]. The importance of the correct estimation of the flow frequency distribution function has been recognized previously as one of the main sources of uncertainty in the damage estimation associated with floods (e.g., [6,60]). However, other authors [61] obtained different results, where the flood damage model is the greatest source of uncertainty for correct flood risk assessment and management, much more than the peak flow frequency distribution. There is no unanimity in the hierarchy of the sources of uncertainty on flood risk assessments, and thus, other studies highlight the importance of other variables, such as the calibration of hydrological-hydraulic models [62], the influence of the quality of the topography [63] and even the correct estimation of the relative height of the main floor of the buildings with respect to the public road [64]. In any case, the results obtained in the present study, similar to those obtained previously by other authors [58], show the importance of correct estimation of the peak flow frequency distribution.
All these uncertainty sources could develop into a lack of confidence in flood hazard and risk maps, therefore making flood risk management tasks more difficult. As we show with our results for Zamora city, uncertainty in peak flow value for a 500-year return period produces important differences in flood-prone areas and in the spatial distribution of flood hazard levels for human loss stability, as well as in direct economic damage for buildings and vehicles. The use of an easy, feasible approach by Monte Carlo analysis allows us to delimitate the uncertainty linked to flow frequency analysis and the absence of a complete knowledge of maximum annual peak flow distribution. This uncertainty, and the variability in the spatial distribution of flood-prone areas is also conditioned by the morphology of the terrain (slope), as can also be seen in our study area ( Figure 6). Thus, the right bank of the Duero River, where the urban area of the city of Zamora is mainly developed, shows a steep slope, and the effect of this uncertainty in the estimation of the peak flow is less evident. On the contrary, the left bank, with a gentler slope and a better developed flood plain, shows more clearly the effects of uncertainty in the estimation of the peak flow, which is reflected in a greater spatial variability of the distribution of flood zones and flood hazard.
Probabilistic flood maps could be a powerful tool in flood risk management since they allow the development of probabilistic Expected Annual Damage (EAD) curves, which in turn would allow better sizing of flood insurance policies. In the same way, the development of probabilistic EADs could improve any further cost-benefit analysis of flood mitigation measures, including self-protection. All this uncertainty reduction could make it possible to improve the homeowner involvement in terms of both participation and financial implications. Regardless, the homeowner's involvement in flood risk management is a hard task, firstly because residents consider flood protection to be an almost exclusive government responsibility [65], and secondly because the profit margin of insurance companies may be reduced due to the increased expenditure in the analysis and control of the self-protection measures implemented and the quality of materials used [66], which may limit the interest of insurance companies in these tasks.
Human loss of stability analysis in Zamora city shows that if we do not take into account the uncertainty of the results, we are probably underestimating the high-hazard areas ( Figure 5). The use of probabilistic maps could allow the flood risk managers to improve the security of population evacuation routes, as well as to delimit the higher hazard areas so that people cannot access them when the occurrence of a flood event is foreseen, or for the evacuation of the population from those areas (thus preventing them from going out into the street). One of the main problems associated with this hazard management is related to the disparity (or variability) of the results obtained, depending on the model used and its complexity. The two models ( Figure 5) used in this analysis can be regarded as being of medium complexity, and they show significant differences in the results, which could be increased if lower [67] or higher [68] complexity models are used. Probably, as is pointed out by Milanesi et al. [68], the analysis of the human loss of stability (which include such factors as slipping, toppling and drowning) is a complicated task, due to the multiple factors involved (surface slope, water density, suspended load, such as tree branches, etc.), but the use of very complicated models may not be absolutely necessary for flood risk management tasks, such as land use limitations.
Finally, regarding hazards for vehicle stability and damage, there are some models too, although the literature here is less mature than that for damage curves for buildings; there are only three studies focused on damages to vehicles in contact with floodwater [69]. According to Martínez-Gomarriz et al. [69,70], the model proposed by USACE [45] and Shand et al. [71] are the more robust and testing. Damage to vehicles can include different processes, such as sliding, toppling, and floating, which could be the origin of vehicle impact against other structures and their resulting damage, although vehicle damage does not require these movements and can be caused only by flooding and the deterioration of the materials of which it is composed. The same holds for human stability (manoeuvrability, abilities, or psychological aspects); vehicle characteristics (weight, shape, or design) determine their stability level [70], and all of them introduce a certain degree of uncertainty and complexity into this type of analysis. Similar to the scenario about human loss of stability, non-consideration of uncertainty in flood hazard by flood risk managers may increase flood damages (as can be derived from the Zamora city results, Figure 7). Once again, the availability and use of probabilistic flood maps could improve the flood damage assessment and management by means of the limitation of vehicle parking areas, both temporary and permanent in time.

Conclusions
Probabilistic flood maps arise as a valuable option to improve flood hazard and risk assessment and management, as they include the uncertainty associated to all of those analyses. In Zamora city (Spain), probabilistic maps were developed with the aim of considering the uncertainty linked to flood frequency analysis, which were calculated from the Monte Carlo method and quantified at around 40% (similar to the previously quantified value in Navaluenga in central Spain, but derived from a different approach).
The flood frequency analysis result uncertainty will be propagated through the hydrodynamic model to the flood-prone area extension and the spatial distribution of flow parameters (depth and velocity), as well as the flood hazard of human loss of stability and vehicle damage.
The results of Zamora city point towards some key points. First, how the terrain surface shape (mainly slope distribution) may increase or reduce the effect (flood-prone area extension) of peak flow uncertainty, so this uncertainty produces few changes in the flooded areas on the right bank of the Duero River, while such changes are much more noticeable on the left bank, with its gentle slope and extensive development of the floodplain.
Furthermore, the changes in flood-prone areas and flood hazard levels of human loss of stability and vehicle damage are not uniform, and they are not equally distributed for all hazard levels. From our results in Zamora, we conclude that changes produce a noticeable increase in the extension and spatial distribution of higher flood hazard areas (which increase around 20% when we compare the results for the quantile and 99% confidence interval), whereas all other hazard levels show a downward trend.
Finally, all the results for Zamora city draw different probability scenarios about flood hazard, which allow flood risk managers to improve urban-use planning from a flood risk point of view. So, managers could define the spatial distribution and restriction of vehicle parking from the probability that one street shows high hazard for vehicle instability (and the probability that those conditions occur), or improve the evacuation paths for the population living in areas where hazard to human loss of stability is higher too. The making of all these decisions by risk managers (who may be politicians who do not have a deep knowledge of flood dynamics) will be facilitated by the availability of this type of probabilistic mapping in which, in a simple way (from a colour scale), one can observe the distribution of the hazard in the urban city environment.