Flood risk assessment through large-scale modeling under uncertainty

The complexity of flood risk models is intrinsically linked to a variety of sources of uncertainty (hydrology, hydraulics, exposed assets, vulnerability, coping capacity, etc.) that affect the accuracy and reliability of the analyses. Estimating the uncertainties associated with the different components allows us to be more confident in the risk values on the ground, thus providing a more reliable assessment for investment, insurance and flood risk management purposes. In this study, we investigate the flood risk of the entire Central Apennines District (CAD) in Central Italy using the laRgE SCale inUndation modEl – Flood Risk (RESCUE-FR), focusing on the interaction between the uncertainty in the hydraulic Manning parameter and the risk variability. We assess the coherence between the quantile flood risk maps generated by our model and the official risk maps provided by the Central Apennines District Authority (CAD Authority) and focusing on three specific zones within the CAD region. Thus, RESCUE-FR is used to estimate the expected annual damage (EAD) and the expected annual population affected (EAPA) across the CAD region and to conduct a comprehensive uncertainty analysis. The latter provides a range of confidence of risk estimation that is essential for identifying vulnerable areas and guiding effective mitigation strategies.

1 Introduction

Floods are one of the world's most devastating natural disasters, posing a significant threat to human life, infrastructure, economy and the environment (Doocy et al., 2013; Rentschler et al., 2022; Llasat et al., 2009). Defined as the temporary inundation of land not normally covered by water (Sayers et al., 2013), floods have shaped human civilization throughout history and continue to challenge societies worldwide. In recent years, the frequency and severity of flood events have increased due to factors such as climate change, land use change and population growth (Schilling et al., 2014; Blöschl et al., 2017; Cred, 2020).

Understanding flood risk, defined as the probability of a flood event combined with its potential adverse consequences, is paramount to effective disaster preparedness, response and mitigation strategies (EU Floods Directive, 2007). Flood risk assessment plays a crucial role in assessing the vulnerability of communities, infrastructure and ecosystems to flood events. It involves estimating the likelihood and potential impact of flooding (defined as hazard), taking into account the exposure (the presence of people, assets and systems in hazard zones) and the vulnerability (the susceptibility of these elements to flood damage) (UNISDR, 2009). Accurate flood risk assessment informs land use planning, emergency management and infrastructure development, helping to identify high-risk areas and prioritize mitigation measures (Falter et al., 2014; Convertino et al., 2019).

Despite its importance, flood risk assessment presents several challenges, particularly at large spatial scales that are fundamental for comparative analyses. Traditional methods for large-scale analysis often rely on empirical or simplified models that may overlook important spatial variability and uncertainty in flood hazard and impact estimates such as a simple aggregation of local flood depth maps and damage assuming homogeneous return periods (Alfieri et al., 2016; Vorogushyn et al., 2018). In addition, the dynamic nature of flood risk, influenced by socio-economic changes and climate variability, adds further complexity to long-term assessments (Field, 2012).

The latest advances in hydrology and hydrodynamic modeling have enabled more accurate and detailed assessments of flood risk at larger spatial scales. These models provide improved representations of flood hazard, ranging from empirical methods to detailed hydrodynamic simulations (Manfreda et al., 2011; Schumann et al., 2013). Among these, LISFLOOD stands out as a widely used model for large-scale flood inundation mapping. It provides valuable insights into flood extent and depth (Bates et al., 2010). In addition, the advances in large-scale flood modeling have contributed significantly to the assessment of flood risk at the continental scale. The Catchment-based Macro-scale Floodplain (CaMa-Flood) model, developed by Yamazaki et al. (2011), is an example of this progress, using advanced hydrodynamic simulations to capture the complex interactions between precipitation, land surface characteristics and river systems. High-resolution topographic data, coupled with advanced modeling techniques, allow for more precise delineation of flood extent and depth, enhancing the reliability of flood risk assessments (Sampson et al., 2015; Yan et al., 2015).

Despite these improvements, flood hazard estimates are still quite uncertain. Uncertainty generally affects the reliability of risk estimates and arises from a variety of sources, including the climate projections, the model parameters and also the socio-economic factors determining the estimate of losses (Dankers et al., 2014; Ward et al., 2014).

The estimation of large-scale flood losses is a complex issue due to the different interactions between systems that are complex, such as urban, suburban and agriculture. In fact Merz et al. (2010b) and Meyer et al. (2013) highlight different types of damage. These include, for example, direct damage (damage caused by contact between water and structures), indirect damage (damage to customers/suppliers who cannot access the flooded area), damage due to business interruption (e.g., due to broken tools as a result of the flood), intangible damage (loss of life, epidemics, environmental damage, etc.) and risk mitigation costs (installation of mitigation infrastructure, maintenance, etc.).

Due to the difficulty of representing all of these types of damage in models, modelers generally limit themselves to assessing direct damage (Merz et al., 2010b) because it can be associated with various hydraulic variables provided by flood hazard models, such as water depth and velocity. These types of direct flood damage models can have different levels of complexity in terms of the variables considered and the spatial scale. They can be univariate, assuming that flood damage is influenced by only one variable (usually water depth), and multivariate, assuming that flood damage is influenced by several variables including flood depth, inundation duration, flow velocity and resistance parameters (e.g., building type, construction material) (Gerl et al., 2016). In terms of spatial scale, they can be micro-scale (the damage is assessed at the building scale), meso-scale (the damage is assessed at the spatial aggregation) and large-scale (the damage is assessed at the municipal, regional or national scale) (Merz et al., 2010b).

Quantifying the overall uncertainty in both hazard and impact estimates is essential to understanding the robustness of results in flood risk assessment and management (USACE, 1992; Peterman and Anderson, 1999; Downton and Pielke, 2005; Dankers and Feyen, 2009; Alfieri et al., 2015, 2016). Incorporating estimation uncertainty allows for a probabilistic assessment of risk, providing a confidence interval for the results obtained rather than a deterministic risk value. This is useful for informed decision-making (Kreibich et al., 2017; Merz and Thieken, 2009; Merz et al., 2010a), in particular where we do not have data to validate our own damage estimates (Figueiredo et al., 2018; Molinari et al., 2020). More generally, this type of estimation supports the decision-making needs of different stakeholders who may have different attitudes to risk or different cost–benefit ratios for risk reduction measures (Merz and Thieken, 2009).

In this paper, we show how the variability in the Manning roughness parameter can affect the risk estimate. We consider here the Manning parameter to be representative of most of the uncertainty characterizing the hydraulic modeling of flood risk. This is only one of the sources of uncertainty that need to be assessed in this context, such as the uncertainty in the estimates of the hydrological discharge, the vulnerability curve and the exposed-asset values. Assessing uncertainty from other sources could provide additional insights and improve the overall confidence, but achieving this requires a balance between the spatial resolution and the uncertainty assessment.

In this context, the use of a large-scale, simple and computationally inexpensive flood hazard model facilitates the parametric uncertainty analysis. To this aim we use the hydrologic–hydraulic modeling approach named RESCUE (laRgE SCale inUndation modEl), which was suited for the specific scope of the probabilistic assessment of flooded areas (Pavesi et al., 2022). Coupled with a damage model, this provides an estimate of the uncertainty in the flood risk. The RESCUE framework combined with the damage model is called RESCUE-FR (RESCUE – Flood Risk). The area of study is the Central Apennines District (CAD) in the center of Italy with a geographic surface of 42 298.22 km². The RESCUE-FR framework is applied to the 25 main river basins in the area, as defined by the Italian National Geoportal, with areas ranging from 125 km² to over 17 000 km². Although this article proposes a case study for Italy, where some damage studies have been carried out for the specific geographic region (see, e.g., Molinari et al., 2014; Amadio et al., 2019), we decided to adopt here the continental damage model of Huizinga (2007) and Huizinga et al. (2017) in order to be able to extend our analysis to a large transnational context in future works. In addition, we estimate the population affected by floods; within this scope, we make the common simplifying assumption that the population is fully affected if located in the flooded area, regardless of the water depth (Ward et al., 2013; Alfieri et al., 2015).

The innovative aspect of this work is twofold:

• A more complete view of flood risk for a large-scale assessment of the associated uncertainties is provided, given by the use of large-scale flood models. This approach is achievable, not through complex models requiring intensive computational processing and calibration parameters but through simpler large-scale models that illustrate the trade-off between quantile uncertainty and spatial resolution/detail.

• The RESCUE-FR framework is developed, which integrates simple, globally accepted hazard, vulnerability and exposure assessment methodologies to ensure wider applicability of flood risk estimates across regional boundaries.

The paper is structured as follows: in Sect. 2 we explain the methods of the RESCUE-FR framework, in Sect. 3 we present the study area, in Sect. 4 we present the main results and discussions, and in Sect. 5 we draw the main conclusions of our work.

2 Methods

2.1 Flood risk analysis

Hazard, vulnerability and exposure are the three fundamental modeling components of flood risk. In this work, we provide a probabilistic assessment of risk by accounting for the uncertainty due to the hydraulic modeling parameter, denoted here as n. To this aim we make some simplifying assumptions; we assume that the probability of flooding (i.e., the hazard) and the hydraulic characteristics of the inundation (e.g., water depth and velocity) are fully determined by the probability distribution function of the peak discharge p(Q) and of the Manning roughness parameter p(n) (where p(.) indicates the probability density function), where Q and n are independent of each other. Moreover, we assume that the damage D (which results from the product of vulnerability and exposure) depends only on the water depth in the flooded area h. Thus, risk R can be described by the following equation:

$$\begin{array}{}\text{(1)}& \begin{array}{rl}R& =\underset{\mathrm{\Omega }}{\int }{R}_{n}p\left(n\right)\mathrm{d}n\\ & =\underset{\mathrm{\Omega }}{\int }\left(\underset{Q^{\ast }}{\overset{\mathrm{\infty }}{\int }}D\left(h\right(Q,n\left)\right)p\left(Q\right)\mathrm{d}Q\right)p\left(n\right)\mathrm{d}n,\end{array}\end{array}$$

where Q^∗ is the minimum peak discharge value producing damages and $R_n=E_Q\left[D\right|n]$ (with E_X[.] denoting expectation relative to the probability distribution function of the variable X) is the expected annual economic loss in a given year, i.e., the expected annual damage (EAD), and for a specific value of n in the range of variability Ω. Thus, $R=E_n\left[R_n\right|n]$ is the average EAD (with respect to the probability distribution of n) and, like R_n, is expressed in monetary terms.

Similarly, we can estimate the average population exposed to the effects of flooding in a given year, i.e., the expected annual population affected (EAPA), which is instead expressed in number of inhabitants. We denote ${\mathrm{\Theta }}_n=E_Q\left[\mathit{\eta }\right|n]$ as the EAPA for a specific value of n in the range of variability Ω and the average EAPA $\mathrm{\Theta }=E_n\left[{\mathrm{\Theta }}_n\right|n]$. These are obtained by substituting the damage function D(h(Q,n)) in Eq. (1) for the function η(Q,n), counting the number of affected inhabitants as a function of the hydrological load Q and of the Manning parameter value n. As previously mentioned, we assume that the affected population does not depend on the water depth.

Note that Eq. (1) provides the average with respect to n variability in EAD or EAPA at any point of the flooded area. Results can be integrated in space at the desired aggregation scale. Given the main objective of this work, in the application of the framework to the case study (Sects. 3 and 4) we particularly focus on the probability distribution of R_n and Θ_n at different spatial scales.

Finally, we used a simplified approach to account for flood protection works. In large-scale flood hazard models, the resolution of the digital terrain model (DTM) is often insufficient to represent flood defense structures, and comprehensive databases detailing such defenses are often not available for manual integration. Previous studies have used various simplified approaches, such as assuming a uniform level of protection (Rojas et al., 2013) or using the gross domestic product (GDP) per capita (Feyen et al., 2012) or the exposure (Zanardo and Salinas, 2022) as a proxy for the level of protection. We adopt here a uniform level of protection by assuming the threshold value Q^∗ in Eq. (1) is equal to the designed return period of the flood protection works. While such an assumption may add an additional component of uncertainty, we surmise that it is less significant than the uncertainty in the Manning parameter as it mainly affects the lower part of the damage distribution, which has less of an impact on risk.

In the next subsections we describe in detail the methodology of the three components of the RESCUE-FR framework.

2.2 The hazard model: RESCUE

RESCUE is a large-scale inundation model that allows for the assessment of flood hazard over large areas (Pavesi et al., 2022). It is a hybrid model combining geomorphological and hydrologic–hydraulic approaches and consists of four distinct components.

1. Geomorphological analysis. Using a DTM, the river network is extracted and segmented into nodes and reaches. This segmentation helps to identify critical points such as breach nodes or channel heads. Each river segment is divided into equal lengths where possible.

2. Cross-section definition. Average cross-sections are defined for each reach identified in step 1. Key hydraulic properties such as the wetted perimeter and area are derived using the Height Above Nearest Drainage (HAND) model (Nobre et al., 2011) as described in Zheng et al. (2018). These properties, together with water levels, are used to calculate key parameters such as the mean wetted area and the hydraulic radius. The Manning equation is then used to generate the rating curve for each segment.

3. Hydrological load. The discharge in all river reaches for a given return period is calculated using the rational formula. This approach requires minimal information on the hydroclimatic conditions of the catchment and includes parameters such as the runoff coefficient and the critical rainfall intensity derived from an intensity–duration–frequency (IDF) curve (see the “Model setup” Section). This approach is used because it is simple and flexible, and although it is usually applied to small catchments, it can be applied to larger catchments if it is properly calibrated.

4. Flood map generation. Flood maps are generated by solving the 1D hydrodynamic model under steady-state conditions. Starting from downstream boundary conditions, the gradually varied-flow equation is solved along the river network to calculate water profiles for each segment. Water levels at segment nodes are interpolated and spatially propagated using the HAND map to identify flooded areas based on the terrain elevation and the water depth.

The framework produces the extent of the flooded area for specific peak discharge values, corresponding to given return period scenarios, and the related values of the water depth h(Q,n) in any DTM cell of the area. To solve the integrals in Eq. (1), that is to account for the probability distribution of the peak discharge and the uncertainty associated with the Manning roughness parameter n, a Monte Carlo analysis is performed, resulting in several flood maps for each return period and Manning value.

2.3 Probabilistic mapping of flood hazard scenarios

The general approach to probabilistic flood mapping involves the evaluation of the statistical distributions for both model inputs (the peak discharge Q) and the parameters. Samples are taken simultaneously from each distribution and used as inputs to numerous model simulations. For each simulation, the inundated area is delineated and the flood inundation is defined in probabilistic terms by analyzing the distribution of the results (e.g., the water depth). In RESCUE we incorporate only the statistical variability in the Manning roughness parameter n, considering it to be representative of most of the uncertainty characterizing the hydraulic modeling of flood risk. Indeed, floodplain roughness is the main factor of uncertainty impacting the hydraulic modeling (Annis et al., 2020). This is represented in RESCUE-FR by the Manning coefficient. It represents a significant proportion of the overall uncertainty in hydraulic flood risk modeling, as it incorporates factors such as variations in terrain elevation, the extent of flood-prone areas and potential channel restrictions due to the presence of infrastructure.

To account for the variability in the Manning coefficient, we perform a Monte Carlo analysis using Latin hypercube sampling (LHS), as described in McKay et al. (2000). Compared to the case of Pavesi et al. (2022), where a uniform distribution of the Manning parameter n was assumed to assess the adaptation of flood hazard in flood-prone areas to official maps, we identified two key points:

1. A uniform distribution for the roughness parameter is used, giving equal weight to all roughness values in a relatively wide distribution, resulting in a certain overestimation of the water depth in the channel and consequently in the floodplains.

2. The elevation of the channel cells obtained from the DTM is uncertain due to the presence of a waterbody. This affects the estimate of the channel bottom slope, which often results in null or very low values, thus introducing an additional source of uncertainty in the estimation of the rating curve.

In general the discussion of the reliability of the rating curves and the key parameters influencing the rating curves is based on the model proposed by Zheng et al. (2018) and has been addressed in the literature by various authors (Garousi-Nejad et al., 2019; Godbout et al., 2019; Johnson et al., 2019; Ghanghas et al., 2022). In particular, from the results of Ghanghas et al. (2022), changes in the Manning roughness and channel slope predominantly affect rating curves during high flows, while the absence of bathymetric data has a significant influence at low flows, with the effect decreasing as discharge increases. In RESCUE, our primary focus is on high flows, allowing us to neglect the bathymetric factor. In addition, the study by Johnson et al. (2019) suggests that wider ranges for varying the Manning parameter can be considered to improve the accuracy of the rating curve estimate and the effectiveness of flood-mapping models.

To address issues related to the channel bed slope parameter, particularly in low-relief areas characterized by very low or zero slopes, we decided to extend the lower limit of the Manning range parameter from 0.03 to 10⁻³ ${\mathrm{m}}^{-\mathrm{1}/\mathrm{3}}\mathrm{s}$, deviating from the approach used in Pavesi et al. (2022). In addition, to overcome the challenges associated with the uniform distribution of the Manning parameter, as discussed in point 1, we chose a normal probability distribution of the Manning parameter. In general, there is no widely accepted probability distribution for the Manning coefficient, as its estimation relies heavily on empirical formulae and expert judgment, making it inherently subjective (Stephens and Bledsoe, 2020). Comparative studies such as that of Papaioannou et al. (2017) show that normal, lognormal, gamma and beta distributions can all be used with a minimal difference in results, leaving the choice of distribution to the expert judgment of the modeler. For the probabilistic modeling of the Manning parameter, we assumed for n a normal distribution with the parameters of μ = 0.0505 ${\mathrm{m}}^{-\mathrm{1}/\mathrm{3}}\mathrm{s}$ and σ = 0.01921 ${\mathrm{m}}^{-\mathrm{1}/\mathrm{3}}\mathrm{s}$. These parameters were chosen to ensure that 99 % of the sample is positive and that the sample values around the mean are consistent with literature values for the Manning parameter n in floodplains as reported in Chow et al. (1988).

2.4 Flood damage model

The results from the hazard model, the water depth in each DTM cell of the flooded area h(Q,n), can be combined with the damage model D(h) to produce probabilistic estimates of risk, according to Eq. (1). The damage model follows the procedure described in Huizinga (2007), which is valid for all European countries. The choice of this type of damage model is guided by the fact that the hazard model is implemented and solved at a large scale and therefore requires a consistent damage model at the same scale (Merz et al., 2010b). In addition, the RESCUE model in combination with this damage model can also be used in an international context in future works. This is introduced considering the product of the vulnerability, which strictly depends on the water depth, and the exposure, i.e., the exposed-asset value.

Besides, with regard to estimating the affected population η(Q,n), we consider the resident population affected for each positive water depth as in Alfieri et al. (2015) and Ward et al. (2013).

After calculating D(h(Q,n)) and η(Q,n) maps, in all the DTM cells of the inundated areas generated by the hazard model, we can derive the R_n and Θ_n maps by numerically solving the inner integral in Eq. (1) cell by cell; then, the empirical cumulative distribution functions (ECDFs) of R_n and Θ_n are calculated for each cell.

In the following subsections the two components of the damage model and the population database used in the analyses are described in detail.

2.4.1 Vulnerability curves

The vulnerability assessment is based on relative (i.e., expressed in percentage terms) depth–damage functions at the European scale. As documented in Huizinga (2007) and Huizinga et al. (2017), these curves are derived from a comprehensive literature review of flood damage data and damage functions for 11 European countries. Five depth–damage curves are established for distinct economic sectors, namely residential, commercial, industry, infrastructure and agriculture. For those countries for which no information on flood vulnerability is available, the vulnerability curves are the mean value between the curves available for the European continent.

2.4.2 Exposure

The exposure is available for each of the 11 European countries and at the European scale; the latter is the average maximum asset damage among all the countries for which we have information. For countries with no information on the overall value of the exposed assets, a methodology is applied to scale the average maximum losses (obtained from countries with available information) at a national scale across different economic sectors using GDP (gross domestic product) per capita PPS (purchasing power standard) obtained from Eurostat. For further explanation, see Huizinga (2007).

In this way, the maximum-damage values are known for all countries. Since these values refer to the year 2007, it is necessary to take inflation into account and to update the maximum-damage values for each country. This adjustment is made by calculating the ratio between the consumer price index (CPI) of the current year and that of 2007. For further clarification the reader is referred to Huizinga et al. (2017).

This procedure is widely accepted as the standard approach for large-scale damage assessment in Europe, as shown by Rojas et al. (2013) and Merz et al. (2010b). The outlined procedure provides clear damage estimates for individual countries at the national level but lacks differentiation at the sub-regional level (NUTS 3). The NUTS classification (Nomenclature of Territorial Units for Statistics) is a hierarchical system used to subdivide the economic territory of the EU for the purposes of the collection, development and harmonization of European regional statistics and the socio-economic analysis of regions, whereby NUTS 1 represents large socio-economic regions, NUTS 2 represents basic regions for the implementation of regional policies and NUTS 3 represents small regions for specific diagnosis. Our study is particularly focused on this differentiation, and therefore we have rescaled the maximum-damage values to the NUTS 3 level. This rescaling was achieved by assigning weights derived from the ratio between the national GDP and the GDP of the respective province.

Figure 1 Area of interest: the Central Apennines District (Italy) divided into municipalities (black boundary); the main cities in terms of population and economic activity are reported (black dots).

The assessment of exposed assets is carried out using the CORINE (Coordination of Information on the Environment) Land Cover (CLC) map, which includes 44 different land cover classes. We adopt the methodology outlined in Huizinga (2007) to assess exposure on land cover grid data. This methodology is based on the integration of data from the European Land Use/Cover Area frame Survey (LUCAS). Following the methodology of Huizinga (2007), statistical mapping aligns the observed LUCAS data with the CLC classes, resulting in a cross-tabulation. This cross-tabulation provides a comprehensive allocation of land use percentages within each CLC cell. Finally, we calculate a damage function for each land cover class as a linear combination of the corresponding damage functions for different land use classes. The weights in the sum are determined from Huizinga (2007).

2.4.3 Population

To assess the population at risk, we chose to use the European baseline dataset HANZE 2.0 (Paprotny and Mengel, 2023). The HANZE (Historical Analysis of Natural Hazards in Europe) dataset was released in 2017, revised and expanded as HANZE 2.0 in 2023. It was the first comprehensive exposure dataset with a resolution matching pan-European flood hazard maps, namely 100 m, covering the years 1870 to 2020, designed specifically to enable the analysis of exposure and land use change within flood assessment studies. This decision was made in favor of the finer grid resolution compared to the 1 km resolution available in other European datasets.

3 Case study: Central Apennines District

To illustrate and discuss the advantages of a probabilistic approach for large-scale flood risk mapping, we apply RESCUE-FR to the area pertaining to the Central Apennines District (CAD). The CAD is a geographical area of 42 298.22 km² located in the heart of Italy (yellow area in Fig. 1). It stretches from the Tyrrhenian Sea to the Adriatic Sea, including the Apennines, and comprises seven regions: Emilia-Romagna, Tuscany, Umbria, Lazio, Marche, Abruzzo and Molise. As shown in Fig. 1, Rome, Perugia, Pesaro, Ancona, Ascoli Piceno, Teramo, Pescara, Chieti, L'Aquila, Viterbo and Latina are among the most important cities (the boundary of each municipality is reported in the figure) in the district in terms of population and economic activity.

The CAD climate is Mediterranean, with rainfall concentrated in the autumn and winter months. The morphology is characterized by mountains, hills and valleys, shaped by the erosive action of the rivers and streams that cross the area. The most vulnerable areas are those located in the floodplains of rivers and their tributaries.

The Central Apennines District Authority (CAD Authority) provides the official flood risk maps for the whole district through the PGRAAC (Piano di Gestione del Rischio Alluvioni dell'Appennino Centrale, i.e., the flood risk management plan of the Central Apennines). These maps are qualitative risk maps that can be considered a reference for the assessment of the RESCUE-FR model. The PGRAAC maps are obtained by superimposing flood hazard maps (which generally result from different hydrologic methods and detailed hydrodynamic models) and a detailed land use vector map. These risk maps do not represent risk in economic terms but are based on the identification of four qualitative risk classes: R1 (low risk), R2 (moderate risk), R3 (high risk) and R4 (very high risk). The CAD Authority maps are based on a risk matrix that considers different return period flood maps and evaluates the land use within the flooded area. Regions with lower return periods and urban land use, with the potential for a large amount of damage to buildings and harm to people, have higher risk (R4 or R3). Conversely, areas with higher return periods and agricultural or forest land use have moderate (R2) to low risk (R1). The four classes, from R1 to R4, are in order of increasing risk severity.

Model setup

Identifying the exact location of rivers is crucial to flood risk assessment. Relying solely on an unconditioned DTM for river extraction can lead to inaccuracies, particularly in planimetry, which may lead to incorrect overlapping with land use types such as residential or others that are not appropriate. To address this issue, we used the MERIT Hydro dataset (Yamazaki et al., 2019), which is hydrologically conditioned globally, ensuring a more accurate representation of river networks. This choice increases the reliability of the analysis by minimizing inaccuracies in river delineation and subsequent overlaps with land use categories for damage estimation.

The RESCUE model setup mirrors the configuration used in Pavesi et al. (2022) in terms of geomorphological analysis and hydrological discharge estimation. A threshold area for channelization of A_c = 10 km² and an average reach length of λ = 1 km are used for geomorphological analysis. The rational method uses the IDF curves at the regional scale based on the two-component extreme value (TCEV) distribution provided by the VAPI project (Rossi and Villani, 1994). The hydrologic model also includes the modified Soil Conservation Service curve number (SCS-CN) approach by Hawkins et al. (2010).

Only in the Adriatic coastal areas, where the SCS-CN method does not accurately reproduce the peak discharge, is the rational approach discarded. Therefore, to ensure reliable discharge estimation in these catchments, we integrated official basin-scale studies obtained from hydrologic–hydraulic reports available online (AUBAC, 2024).

We also deviate from Pavesi et al. (2022) for the hydraulic Manning distribution. We choose a normal distribution for the Manning parameter n with a mean of μ = 0.0505 ${\mathrm{m}}^{-\mathrm{1}/\mathrm{3}}\mathrm{s}$ and a standard deviation of σ = 0.01921 ${\mathrm{m}}^{-\mathrm{1}/\mathrm{3}}\mathrm{s}$. These parameter choices were made on the basis of the discussions outlined in Sect. 2.3.

The analysis was carried out on the 25 main catchments defined by the Italian National Geoportal. The simulations were run separately for each of these basins. The basin to be simulated is then defined a priori, and the area in which the simulation is carried out coincides with the area of the basin. In order to assess risk and damage scenarios in each catchment, we have discretized the peak discharge return period into chosen intervals ranging from a 100- to 1000-year return period. In this regard, for each return period, we conducted 11 simulations using Latin hypercube sampling (LHS) to ensure that each sample was representative of a certain probability of occurrence of the Manning parameter. In a single simulation, the same Manning coefficient is used across the entire spatial domain.

To account for flood protection works, we assumed a uniform level of protection in all areas. In Italy the areas exposed to higher potential damage are generally protected by flood defenses up to a return period value of 200 years. Assuming a uniform flood protection return period of 200 years in the whole area may lead to a significant underestimation of the risk in those areas, whose extent is not negligible, that are effectively unprotected. For this reason, we look for a trade-off between these two issues and choose a uniform flood protection return period of 100 years. The choice of this return period is in agreement with the studies of Ward et al. (2013) and Rojas et al. (2013). Hence, we assume $Q^{\ast }=Q_{\mathrm{100}}$ in Eq. (1) to mitigate the overestimation of flood risk associated with the lower return periods, where Q₁₀₀ is the $\mathrm{1}-\mathrm{1}/\mathrm{100}=\mathrm{0.99}$ quantile of the probability distribution of the peak discharge (the discharge for a return period of 100 years).

As described in Sect. 2.4.1, we adopted the European-scale vulnerability curves following the procedure outlined by Huizinga (2007). The NUTS 3 level was considered to be the most appropriate scale for our case study, which focused on the Central Apennines District. It was therefore essential to assess differences between different provinces rather than regions (NUTS 2 level). To ensure consistency with the methodology presented in Huizinga (2007), where a cross-tabulation between LUCAS and CLC data is presented, we have updated the maximum-damage values for each category from Huizinga (2007) to current values adjusted for inflation. In addition, we have chosen to use the CLC 2000 dataset (Büttner et al., 2002) and refer to the land use percentages in each cell as reported in Huizinga (2007).

Figure 2 Comparison between (a, c, e) the PGRAAC maps and (b, d, f) the RESCUE-FR 50th quantile of R_n risk maps: zoom of the (a, b) upper Tiber River near the city of Perugia, (c, d) middle Tiber Valley near Rome and (e, f) entire Musone River in the region of Marche. See the text for the definition of the PGRAAC and RESCUE-FR risk classes.

4 Results and discussion

In this section, we present and discuss the results of the analysis of flood risk estimation uncertainty carried out in the CAD case study. First, we focus on evaluating the consistency between the flood risk maps generated by RESCUE-FR and the official risk maps provided by the CAD Authority. This comparison is essential to assessing the reliability of RESCUE-FR; results are illustrated in Fig. 2 for the whole CAD region.

In Fig. 2, we compare the PGRAAC maps (Fig. 2a, c and e) with the RESCUE-FR risk maps at the 50th quantile of risk (the median of R_n distribution) (Fig. 2b, d and f), expressed in euros per year. Note that agricultural land use has significantly different maximum-damage values compared to the other land use types, varying by about 2 orders of magnitude because the maximum-damage value for agriculture in European countries is on average lower by 2 to 3 orders of magnitude than that associated with residential, commercial and industrial land uses (Huizinga, 2007). To illustrate these differences, we have chosen a logarithmic scale for risk visualization (per order of magnitude), in line with the qualitative approach used by the CAD Authority. To simplify the representation and favor the comparison with PGRAAC maps, we assign three different risk classes based on the loss entity: C1 risk class for a value up to EUR 100 per year, C2 up to EUR 1000 per year and C3 above EUR 10 000 per year. Our classes, no matter how they are defined, do not correspond to the R1–R4 risk classes used by CAD that are not defined in numeric values.

Figure 2a and b focus on an area of the upper Tiber River at the confluence with the Chiascio River near Perugia. Here, the two flooded areas (RESCUE and PGRAAC) differ significantly, especially along the Chiascio River. In particular, RESCUE tends to underestimate the flood extent for higher return periods compared to PGRAAC; this is due to errors in the DTM representation of the valley (results not shown). However, in the case of forest/agricultural land use, despite a significant underestimation of the flooded area, the error in risk assessment due to underestimation of the flooded area by the hazard model is not very significant because the flooded area underestimated by RESCUE is in the R1 low-risk class in the PGRAAC risk map. This would be different if the land use was of a different type (residential, industrial, etc.), which would result in a higher-risk class and make the local underestimation more relevant.

Figure 2c and d show a zoomed-in section of the middle Tiber Valley, revealing very good agreement between the two maps in the transition from moderate-/low-risk to high-/very high-risk zones. However, we do notice some differences in high-risk areas, which could be underestimated, particularly in zones adjacent to the river channel. This discrepancy is primarily due to the coarser resolution of the RESCUE-FR model. The CLC dataset, with a spatial resolution of 100 m, cannot capture detailed land use variations within a single cell compared to the land use vector map used by the CAD Authority. Therefore, if a cell has a small percentage of residential or industrial land use but the majority is agricultural, it will be classified as agricultural damage. However, the resulting underestimation of risk compared to the total risk is almost negligible. Finally, in Fig. 2e and f we show risk maps of almost the entire Musone River, where we can draw similar conclusions to those drawn from Fig. 2c and d. The land use of the CAD area at risk is shown in Fig. 3, where Fig. 3a–c are focused, respectively, on the area at risk covered by Fig. 2a and b, Fig. 2c and d, and Fig. 2e and f.

Figure 3 Land use of the Central Apennines District (Italy) area at risk: zoom of the (a) upper Tiber River near the city of Perugia, (b) middle Tiber Valley near Rome and (c) entire Musone River in the region of Marche. Classes are given in the legend.

In general, RESCUE-FR shows consistent agreement with the official CAD maps, with the exception of areas that cannot be accurately mapped by the hazard model due to limitations of the DTM and the spatial resolution of the CLC, which may miss localized areas of land use. In the following we analyze and discuss the results of the uncertainty analysis.

Figure 4a shows an illustrative median R_n risk map of a segment within the Aterno-Pescara basin, derived from the RESCUE-FR model analysis. Similar to Fig. 2, the legend is log-scaled to highlight the variation in risk between different land use types within cells. The purpose of this figure is to highlight the variability due to parameter uncertainty that is reflected in the R_n probability curve and to understand the added value of the uncertainty assessment. In general, three patterns of risk distribution were identified throughout the analysis. These are shown in Fig. 4b–d. To illustrate these patterns, we select three specific cells with different land uses, different empirical cumulative distribution functions (ECDFs) and different distances from the channel along a cross-section of the risk area.

Figure 4 R_n empirical cumulative distribution function (ECDF): (a) zoom of the RESCUE-FR median risk map of the Aterno-Pescara River, (b) local distribution in a cell near the main channel (commercial/residential land cover), (c) local distribution in a cell far from the channel (agricultural land cover) and (d) local distribution of the risk that a cell is very unlikely to get wet (commercial/residential land cover). For a description of the CORINE Land Cover (CLC) classes see, e.g., Büttner et al. (2002).

In Fig. 4b, we show the R_n distribution of cell 1, which is a typical cell close to the main river channel and therefore always wet, with water depths often above 6 m (maximum damage). The distribution trend is asymptotic, and there is no variation in the inter-quantile range. The risk values are the values of a cell which remains in commercial/residential land use, such as CLC class 1.1.2. For a description of the CORINE Land Cover (CLC) classes see, e.g., Büttner et al. (2002).

In Fig. 4c, we see a different pattern of the risk R_n distribution, for a typical cell far from the main river channel that starts to be wet from about the 7th quantile and with a large variation in its inter-quantile range. The R_n values are the values of a cell which remains in agricultural land use such as CLC class 2.4.1. For this type of land cover class the risk value is very low compared to the land cover class representing the residential/commercial risk. The variability for agricultural land cover is therefore negligible in absolute economic terms.

Figure 4d, on the other hand, shows a typical cell which, despite being very close to the channel, is topographically much higher than the channel bottom and therefore only flooded when the water depth reaches very high levels. Here we can see that if we did not use such an approach to quantify the parametric uncertainty, we would greatly underestimate the risk associated with a certain area. In fact, in this case, the land use is commercial/residential, and therefore the relative risk is very high if the cell is flooded. Without this approach, we would probably have rated the risk in this cell as zero.

To summarize, there is considerable variability in the risk values provided by large-scale maps due to the uncertainty in the Manning parameter, embedding the parametric uncertainty in the hydraulic modeling approach. Thus, it is important to understand how these risk values R_n can be used effectively. A clear and comprehensive view of the extent of risk variability across the whole area is essential in order to identify the critical points across the area of interest and thus delineate the areas where risk reduction efforts should be focused and their priority.

Figure 5 Expected annual damage R_n distribution at the municipality scale: (a) median quantile and (b) 90th quantile. N/D: no data.

Figure 6 Distribution of the expected annual population affected Θ_n (hab, inhabitants) at the municipality scale: (a) median quantile and (b) 90th quantile.

Therefore, in Figs. 5 and 6 we present the main characteristics of the probability distribution of the expected annual damage R_n and the expected annual population affected Θ_n aggregated at the municipal level; for both, we show the median and the 90th quantile. Note that the boundary of each municipality is depicted in black as in Fig. 1. It is important to stress that, as mentioned in Sect. 3.1, the analysis was carried out only for the 25 largest catchments. Since within some municipalities there are areas pertaining to catchments that are not simulated with RESCUE-FR, we have excluded these municipalities from the discussion if the non-simulated area is 100 % (gray municipalities in the figures). In all the other cases, where only a fraction of the municipality area is not simulated, the risk (in terms of EAD or EAPA) in that municipality is generally underestimated. This is visible in the maps because the non-simulated area within the municipality boundary is in gray. The gray areas are concentrated in the coastal areas, where the catchments are generally small. In general R_n increases with increasing Manning roughness because greater roughness leads to deeper water depth. The water depth increases monotonically with the damage in the damage model adopted.

From the analysis conducted, it emerges that the median EAD in the CAD is approximately EUR 137 million, with an average exposed population of around 5350 individuals annually. When the 90th quantile of the total EAD and EAPA distributions is assessed, these values increase by 34 % and 27 %, respectively.

The values in Figs. 5 and 6 show a wide range of annual losses in terms of EAD and EAPA. As expected, the highest risk values are found in the major economic and populated cities (as reported in Sect. 3), where there are more developed economic activities and a larger resident population. In particular Fig. 5a shows that 63 % of the CAD exhibits a median risk of less than EUR 1000 per year or even null. For another 27 %, the risk appears to be less than EUR 100 000 per year, while for the remaining 10 %, it exceeds EUR 100 000 per year. The reason for this large percentage at low risk is that the highest risks are in cities passing through main rivers, which have larger discharges than tributaries.

Moreover, by comparing Fig. 5a and Fig. 5b, it emerges that a different number of areas have a significant increase in EAD at the 90th quantile compared to the median. This allows for quickly identifying areas of high uncertainty and carefully deciding where to focus for further assessment or for risk management actions. This provides an immediate indication of areas that may require priority action to reduce the impact of potential flood-related disasters.

As expected being the largest city in the area, the municipality of Rome stands out as having the highest median risk, estimated at around EUR 82 million per year. This corresponds to 60 % of the total risk in the CAD. This value is in line with the results of a previous study (Fiori et al., 2023), which provided an estimated risk of EUR 24.3 million per year, the same order of magnitude. The difference is justified considering that the risk estimate of Fiori et al. (2023) results from the implementation of a detailed 2D hydrodynamic model, limited to the historic center of the city of Rome rather than the entire municipality. For the same study area mentioned above, our model has an average value of around EUR 43 million. Despite a simplified representation of the phenomenon, the model is able to provide reliable estimates, also enabling the assessment of parametric uncertainty.

As for the population at risk, Fig. 6a shows that 93 % of the CAD has an average population affected of fewer than 10 persons per year or no risk at all. For another 6 %, the risk appears to be fewer than 50 persons exposed per year, while for the remaining 1 % it exceeds 2500 persons exposed per year. Furthermore, it is interesting to note that in almost 50 % of the simulated area the population at risk is equal to zero (light-green areas). This can be explained by the fact that in the catchments that do not cross large urban centers, such as cities shown in Fig. 1, the population is concentrated in the mountainous areas, at high altitudes above the rivers and not affected by floods. Finally, the 90th quantile does not show significant variability in the population exposed at the municipal level, since this varies only in relation to the flooded area and the intersection of this area with the population distribution, without taking into account the variation in water depth.

5 Conclusions

In this study, we present and discuss the assessment of the parametric uncertainty in flood risk estimation performed with the large-scale risk model named RESCUE-FR. The model couples hazard assessment with economic damage estimation through widely recognized methodologies and can be easily extended to the entire European continent. More importantly, thanks to its simplicity, RESCUE-FR allows for investigating uncertainty without excessive computational efforts. We analyze here the effect of the variability of the Manning parameter, being representative of most of the different sources of uncertainty related to hydraulic modeling. Further, we demonstrate how the resulting distribution of risk can be used effectively for risk assessment and management. In fact, the methodology provides the distribution risk/damage values in economic terms and not a deterministic value of risk/damage, as is often the case with existing risk maps that are based on detailed models and provide deterministic outputs. The proposed approach takes into account parametric uncertainty and uses a large-scale hydrologic–hydraulic analysis that ensures a comprehensive and homogeneous, albeit simplified, representation of processes across the region. This provides more informed results for stakeholders and decision-makers and allows for expanding to other areas where these maps are lacking. The risk analysis was carried out for the CAD region (Central Italy). From this analysis we draw the following main conclusions.

• The RESCUE-FR model shows consistent agreement with the qualitative risk maps provided by the CAD Authority. There are areas of discrepancy, where we generally observe an underestimation of the risk, which are mainly due to the accuracy of the DTM and the limitations in the detailed mapping of land use changes given by the coarse resolution of the CORINE Land Cover.

• Through the uncertainty analysis, the RESCUE-FR model provides a range of estimated flood risk (the expected annual damage), revealing significant variability in risk values depending on the land use considered.

• We identified three different patterns of risk distribution at the cell scale, depending on the distance from the canal, the orography and the different land use types. These patterns are useful to identify the critical points across the area of interest, where a single-value estimate of risk could significantly underestimate risk.

• In terms of the affected population, the variability is lower as it depends mainly on the extent of flooding, which is less influenced by the Manning parameter with respect to water depth. In this case, the median value could be sufficient for the assessment of the expected annual population affected. Conversely, the expected annual damage has a higher variability being directly linked to water depth, requiring a thorough examination of the distribution quantiles for a comprehensive risk assessment of the whole area.

• Aggregating risk values at the municipal level provides a comprehensive overview of flood risk across the CAD region. This approach identifies areas of vulnerability and variability in risk levels, helping to prioritize risk reduction efforts.

• The comparison with a previous study for the municipality of Rome demonstrates the reliability of the RESCUE-FR model in providing valuable information on flood risk, despite its simplified representation compared to more detailed hydrodynamic models.

We recall that the accuracy of RESCUE-FR estimates might be limited by the simplifying assumptions that, on the other hand, allow for an easy implementation of the uncertainty analysis. These include the spatial resolution of the maps used in the input (DTM and CLC). Although these maps allow for consistent use of information at a larger scale, the accuracy of the maps could be improved by a finer resolution at the cost of an increased computational load. Furthermore, the use of open databases containing the distribution of major defense structures, such as levees and dams, could improve the accuracy of risk estimates.

Finally, by assessing the uncertainty associated with individual parameters, we gain valuable insight into interpreting flood risk and understanding map results over large areas. However, it is important to recognize that uncertainty comes from multiple sources beyond hydraulic parameters. Therefore, extending our analysis to assess different sources of uncertainty, including those related to hydrologic and socio-economic modeling, will make the analysis more complete, as it would take into account and compare different phenomena. This is far beyond the scope of this work and will be the subject of future analyses.