1 Introduction

Flooding is considered one of the most disastrous phenomena that affect many parts of the world, causing losses of human lives and infrastructure and property damage (Teng et al. 2017; Mokhtar et al. 2018). Flood events can be induced by both short-duration and high-intensity or prolonged and low-intensity rainfall events. Due to climate change, local intense rainfall is becoming increasingly frequent, resulting in a rapid increase in river levels that contribute to fluvial/riverine floods. In urban environments, such events induce pluvial floods in the case of excessive runoff caused by exceeding the capacities of drainage systems (Lin et al. 2019). To define flood hazard (FH) and related zones, inundation maps depicting water depths related to the characteristic return periods of triggering events (i.e. rainfall) are usually considered. The boundary condition for this evaluation is a reference scenario generally corresponding to an extreme flood event (Woo and Waylen 1986). This reference extreme flood scenario is commonly identified based on hydrological historical data (i.e. rainfall or discharge data, e.g. Sutcliffe 1987), or, in their absence, on geomorphological field observations (e.g. Furdada et al. 2008; Magliulo and Cusano 2016; Montané et al. 2017). Comprehensive analyses of hazards that provide information on the magnitude and frequency levels of events allow for evaluating the risks and quantitatively estimating the expected annual damage induced by a specific event. For this purpose, in recent years, many countries have produced flood risk maps for their territories considering a specific return period range. In Italy, following the European Flood Risk Directive 2007/60/EC of 2007, the Ministry of Environment and Land Protection (now the Ministry of the Ecological Transition) has provided guidelines for implementing FH and risk maps at the hydrographic basin scale.

Recently, many methodologies have been developed for FH assessment to improve predictions of the spatial extent, depth of water level, and flood frequency (e.g. Cook and Merwade 2009; Ongdas et al. 2020; Kim et al. 2021). Empirical models are often based on various statistical and data-driven approaches; they can be categorized into several methods, such as multicriteria analysis (e.g. Chen et al. 2014), which is developed using multiple indicators based on statistical laws and expert knowledge; statistical methods (e.g. Nuswantoro et al. 2016; Guerriero et al. 2020a, b), which include bivariate (e.g. Youssef et al. 2015; Costache 2019; El-Magd 2019) and multivariate (e.g. Shafizadeh-Moghadam et al. 2018; Tien Bui et al. 2019) models; and machine learning methods (e.g. Wang et al. 2015; Al-Abadi 2018; Janizadeh et al. 2019; Costache 2019; Rahman et al. 2019). In general, the empirical approach may be associated with significant uncertainty due to the unavailability of sufficient historical data and high-quality field records (Lin et al. 2019; Tufano et al. 2019). Numerical hydrologic and hydrodynamic models (e.g. Liang et al. 2016; Wang and Yang 2020) simulate physical hydrological processes leading to flooding events through scenario simulation analyses (Bates and De Roo 2000; Di Baldassarre et al. 2010; De Risi et al. 2015; Gonçalves et al. 2015; Toda et al. 2017). Although numerical models require a significant amount of input data (Tufano et al. 2016, 2021) and are usually applied to small areas to replicate the physical processes behind water flow (Ji et al. 2012; Chen et al. 2016; Carmo 2020), they provide effective results that allow us to better understand the hydrological processes associated with flooding generated by rainfall (Rashid et al. 2016). In this context, hydrological input data are generally represented by runoff hydrographs. Hydrographs are often not directly available due to a lack of gauging stations, but multiple techniques have been developed to estimate synthetic unit hydrographs representative of river dynamics (Salami et al. 2017). One of the most common event-based rainfall–runoff methods is the soil conservation service curve number (SCS-CN) method (SCS 1956, 1964, 1972, 1993), for which the peak discharge of stream flow is obtained from a synthetic unit hydrograph for a given rainfall event through convolution. Due to the spatiotemporal variability of rainfall, the quality of measured rainfall data, and the variability of watershed characteristics (i.e. soil type, land use, surface condition, and antecedent soil moisture), the SCS-CN method can exhibit variability in runoff computations (Ponce and Hawkins 1996).

In large river basins, flood events induced by rainfall can be sustained by contributing areas of variable extent as functions of the spatial concentrations of the triggering meteorological events. Flood hazard assessment at the basin scale often considers flood events generated by the overall basin extent; therefore, the potential for events sustained by a fraction of the basin surface is not considered. This consequence might represent a limitation in the evaluation since, in such systems, the magnitude of a flood is related to both the return period and spatial extent of the triggering rainfall event. Deriving multiple flood hazard scenarios might be significant, especially when associated with a statistic of flood-contributing areas for the selected basin. On this basis, in this paper, a multiscenario flood hazard estimation procedure accounting for the potential variability levels of contributing areas in rainfall-induced flood development is presented and tested. The procedure is based on a probabilistic analysis of a long-term rainfall series of annual maxima for intensity–duration–frequency (IDF) curve generation, followed by surface runoff estimation by the SCS-CN rainfall–runoff method considering extents and characteristics of multiple contributing areas, hydrodynamic modelling through the hydrologic engineering centre river analysis system (HEC-RAS) model, and finally, flood hazard map production by inundation scenario overlay.

To test its potential, this procedure is applied to the archaeological site of Sybaris (southern Italy), which represents the most important Archaic Greek settlement of Magna Graecia in southern Italy and is in a floodplain susceptible to floods (Petrucci and Polemio 2007). This site has been involved in several destructive riverine flood events, such as those that occurred in 2013 and 2018. The choice of this site is based on the growing attention to the protection of cultural heritage sites from geological phenomena and the ongoing need for specific methodologies supporting risk evaluation and mitigation measurement planning for these cultural assets (Iriarte et al. 2010; Jigyasu et al. 2014; Rispoli et al. 2020; Valagussa et al. 2020). Especially in Italy, out of a total of 205,670 cultural heritage sites, approximately 16,000 sites (~ 8% of the total) are exposed to high flood hazards. Of these, 804 sites are in the Calabria region of southern Italy (Trigila et al. 2021). The results from procedure application indicate (1) its potential in supporting the understanding of the hydrological processes behind flood initiation, (2) the reliability of predictions of the rainfall–runoff processes, and (3) the significance of simulating the dynamics of overland flows in flood hazard estimation.

The paper is structured as follows: Sect. 2 provides an overview of the test site. The methods section describes the probabilistic rainfall and discharge analyses (which are used to derive different scenarios of the event), the SCS-CN method for runoff hydrographs estimation, the hydrodynamic model HEC-RAS, and the application of the approach to the study area. In Sect.  4, the obtained results are described and discussed. In Sect. 5, the obtained results are validated, and Sect. 6 provides concluding remarks.

2 Study site

The test site selected for the application of the developed procedure is the Sybaris archaeological site located approximately three kilometres from the Crati River mouth in the alluvial plain of Sibari, in southern Italy (Fig. 1a). From a geological point of view, the plain is between the Calabrian Arc and Southern Apennines within a complex geologic context characterized by the overlapping of thin Calabrian Arc units on the Apennine units. Calabrian Arc units, which are outcropped along the southern margin, consist of various igneous and metamorphic rocks (Liberi and Piluso 2009), while Apennine units are composed of a Mesozoic-Tertiary carbonatic succession and flysch deposits overlapped by a Plio-Pleistocene siliciclastic succession (Selli 1962; Caruso et al. 2013). The Sibari Plain originated from tectonic activity (Pliocene–lower Pleistocene; Ambrosetti et al. 1987) and was subsequently filled by fan deltas and marine terrigenous sediments. In particular, the upper part is occupied by Holocene alluvial deposits, mostly originating from rivers that cross the plain (e.g. Crati River), with a maximum thickness estimated up to ~ 100 m (Cianflone et al. 2018). The consolidation process is responsible for the strong subsidence (~ 3 mm/year) that affects the plain, both in historical and present times (Ferranti et al. 2011; Cianflone et al. 2015). In the past, the plain was characterized by very active fluvial dynamics, including different stream captures and drainage diversions, as testified by the numerous paleoriverbeds in the area (Cianflone et al. 2018).

Fig. 1
figure 1

a General framework of the study area with the indication of the Crati Basin (purple circle indicates the Torano Scalo rain gauge; cyan triangles indicate the Crati, Coscile, and Esaro hydrometric stations). b Position and structures of the Sybaris archaeological site.

Full size image

The main river that crosses the plain is the Crati. Its basin is 2460 km2 wide, with a mean discharge at the mouth of ~ 36 m3/s and a torrential regime with alternating abrupt, high-discharge flood periods and dry periods in the summer season. The Crati is affected by the presence of some artificial lakes along its course (i.e. Tarsia and Esaro dams), while artificial embankments along its lower course limit flood potential. Due to the wide extension of the Crati Basin, it contains many tributaries, the most important of which are the Coscile and Esaro Rivers, with mean discharge levels of ~ 14 m3/s and 8.9 m3/s, respectively.

The Sybaris archaeological site is in the ancient Greek town of Sybaris, founded in 720 B.C. Sybaris is considered among the most important archaeological finds of the Archaic and Classical ages in the Mediterranean area. After the decline of Sybaris in approximately 500 B.C., the Hellenist town of Thurii (approximately 440 B.C.) and the Roman town of Copia (approximately 190 B.C.) were successively developed, overlapping the remnants of the older towns (Greco et al. 1999). The most recent of the three is better represented within the archaeological site due to its position near the ground surface. Historic outcropped settlements are diffused over the four major sectors of the site: the Parco del Cavallo, the Prolungamento Strada, the Casa Bianca, and the Stombi (Fig. 1b).

Moreover, due to its low topographic position (between 2 and 6 m below the ground surface), the ongoing subsidence (Cafaro et al. 2013; Bianchini and Moretti 2015; Cianflone et al. 2015), and the proximity to the coastline, the site is susceptible to groundwater flooding events, which are prevented by a set of well points that continuously work to maintain the water table below the level of the archaeological site.

3 Materials and methods

As reported in Fig. 2, the proposed procedure of multiscenario flood hazard assessment is based first on the probabilistic analysis of rainfall time series to reconstruct IDF curves. Second, this assessment is based on the curve number (CN) estimated on the basis of land use, hydrologic soil group (HSG), and antecedent moisture conditions (AMC) of soils. Considering multiple contributing areas, the procedure consists of runoff hydrographs calculation by the SCS-CN method and the hydrodynamic modelling by the HEC-RAS model, providing maps of depth and extension of flow. Finally, the approach is based on FH maps assembled by overlaying the flood inundation maps.

Fig. 2
figure 2

Flowchart of methodologies applied for FH mapping

Full size image

4 FH estimation

4.1 Probability analysis of hydrological data

As a first step for FH estimation, the annual maximum rainfall data were fit by the type I generalized extreme value (GEV) distribution (Gumbel; e.g. Jenkinson 1955) for a specific return period. The probability function F(h) was expressed as follows:

$$F (h) = \exp \{ - \exp [ - (\alpha (h - u))]\}$$
(1)

where h is the maximum rainfall height and the parameters α and u are estimated considering the mean value of the series (μ) and its standard deviation (σ), respectively. They were calculated through the following equations:

$$\alpha = \frac{1.283}{\sigma }$$
(2)
$$u = \mu - 0.45 \sigma .$$
(3)

F (h) represents the probability of not exceeding the value of the rainfall height h by a random variable. Fixed a return period T > 1, the inverse of the cumulative distribution function (CDF; i.e. a survivor distribution function) provides the values of h on the basis of the return period T and duration d:

$$h\left( d \right) = u\left( d \right) - \frac{1}{\alpha \left( d \right)}\ln \left( { - \ln \left( {1 - \frac{1}{T}} \right)} \right) .$$
(4)

Equation (4) allows us to define the height–duration–frequency (HDF) curves of rainfall, generally described as follows:

$$h\left( {d, T} \right) = a \left( T \right)d^{n}$$
(5)

where a and n are parameters estimated through a graphical and probabilistic approach, respectively.

IDF curves were derived directly by HDF curves for the selected return periods of (1) 10 years, (2) 50 years, (3) 100 years, and (4) 300 years.

Conversely, the annual discharge data were fit by a gamma distribution function (e.g. Yue 2001; Guerriero et al. 2018) to estimate the probability of exceedance (p) and the related return period (T). The related return period had the following form:

$$F \left( x \right) = \frac{{\beta^{\alpha } x^{\alpha - 1} e^{ - \beta x} }}{\Gamma \left( \alpha \right)}$$
(6)

where x is the maximum discharge, α is the shape, β is the scale parameter of distribution, and Γ(α) is the gamma function, which was calculated as follows:

$$\Gamma \left( \alpha \right) = \mathop \int \limits_{0}^{x} x^{\alpha - 1} e^{ - x} {\text{d}}x$$
(7)

As indicated by Guerriero et al. (2020b), the generalized extreme value distribution function was suitable for describing the statistical behaviours of annual maxima; however, it had the disadvantage of underestimating the intensities of very high-return period events. For this reason, the gamma function could be useful to overcome this drawback, especially for time series presenting problems with continuity and length.

4.1.1 Rainfall–runoff model description for hydrograph estimation

The SCS-CN method (SCS 1956, 1964, 1972, 1993) is adopted for constructing synthetic unit hydrographs based on a dimensionless hydrograph, which relates ratios of time to ratios of flow (for more details, see Ramírez 2000). The method is based on the following equation:

$$\frac{{P - I_{{\text{a}}} - Q}}{S} = \frac{Q}{{P - I_{{\text{a}}} }}$$
(8)

where P is the cumulative rainfall (mm), Ia is the initial abstraction (mm), Q is the runoff and S is the potential maximum soil moisture retention after runoff initiation (mm). With respect to runoff, Eq. (8) takes the following form:

$$Q = \frac{{\left( {P - I_{{\text{a}}} } \right)^{2} }}{{\left( {P - I_{{\text{a}}} } \right) + S}}$$
(9)

Further simplifications of the SCS-CN method are represented by the average value for Ia, which is selected as 0.2S, and by the introduction of a parameter called curve number (CN):

$$S = \left( {\frac{1000}{{CN}} - 10} \right).$$
(10)

CN varies between 30 for low runoff potential (and high infiltration) and 100 for high runoff potential, depending on land use, HSG, and AMC. Four HSGs to account for the infiltration abilities of different soils are considered: Group A (low runoff potential); Group B (moderate infiltration rates); Group C (slow infiltration rates); and Group D (high runoff potential). Once the group is defined, CN is evaluated based on land use. Moreover, due to the variability in the antecedent rainfall conditions and the associated soil moisture amount, the CN value exhibits variability in the recognition of the AMC. Specifically, three levels of AMC could be considered: AMC I (dry), AMC II (normal), and AMC III (wet). From a statistical point of view, AMC I, AMC II, and AMC III correspond to 90, 50, and 10% cumulative probability of exceedance of runoff depth for a given rainfall, respectively (Hjemfelt et al. 1982).

The synthetic unit hydrograph is expressed as the ratio of runoff Q and peak of discharge Qp, and the ratio of time t (duration of rainfall) and time to peak tp (Ramírez 2000). The ratio of Q to Qp increases over time and reaches a maximum value (equal to 1) when Q equals Qp; then, it decreases until it reaches a value near zero. The time required to reach the peak and return to zero defines the duration of the hydrograph.

The peak discharge, modulated by the basin area (a) and runoff (Q), is obtained through the equation:

$$Q_{{\text{p}}} = \frac{0.208 \times a \times Q}{{t_{{\text{p}}} }} .$$
(11)

The time to peak tp is calculated as follows:

$$t_{{\text{p}}} = \frac{{t_{{\text{c}}} + 0.133 t_{{\text{c}}} }}{1.7}$$
(12)

where tc is the time of concentration (min).

In the technical literature, there are several empirical equations for estimating the time of concentration with good applicability levels because of the limited information they require.

The values of both the peak discharge and time to peak are applied to the dimensionless hydrograph ratios by the SCS-CN method, allowing us to obtain a unit hydrograph (e.g. Raghunath 2006). Subsequently, the established unit hydrograph is used to develop runoff hydrographs considering cumulative rainfall derived by IDF curves for selected return periods.

Moreover, due to the uncertainty linked to the spatial concentrations of the triggering meteorological events, three variable contributing areas are considered. For this purpose, a topographic analysis is conducted in a GIS environment to delineate a defined watershed area. The algorithm used is an eight-direction pour point algorithm that starts from the digital elevation model (DEM) of the area and the drainage network.

4.1.2 2D hydrodynamic modelling

Runoff hydrographs developed considering IDF curves and contributing area features are used for multiscenario hydrodynamic simulations conducted through the HEC-RAS model (Brunner 2021). The software simulates the movement of water over a topographic surface by solving the Saint–Venant equations:

$$\frac{\delta A}{{\delta t}} + \frac{\delta Q}{{\delta x}} = q_{l}$$
(13)
$$\frac{\delta Q}{{\delta t}} + \frac{{\delta \left( {\frac{{Q^{2} }}{A}} \right)}}{\delta x} + gA \frac{\delta H}{{\delta x}} + gA\left( {S_{0} - S_{f} } \right) = 0$$
(14)

where A is the cross-sectional area normal to the flow, Q is the discharge, g is the acceleration due to gravity, H is the water surface elevation, t is the temporal coordinate, x is the longitudinal coordinate, ql is lateral inflow from tributaries (in our case equals 0), S0 is the channel bottom slope, and Sf is the friction slope. Equations (13) and (14) are solved using the four-point implicit box finite difference scheme.

4.2 Application of the proposed procedure to the Sybaris archaeological site

The proposed procedure of FH analysis was tested through an application at the archaeological site of Sybaris. To this end, the data collected for the study area consisted of: long rainfall time series, three discharge time series, a 5 × 5 m DEM, the Crati Basin drainage network, and a land use map.

Rainfall data (Fig. 3a), in the form of annual maxima for specific durations (1, 3, 6, 12, and 24 h), were derived by the Torano Scalo rain gauge (see Fig. 1a for position) measurements, recorded between 1934 and 2012 (https://www.cfd.calabria.it/). There were 60 units of data, and discontinuities were present; therefore, only the data between 1962 and 2012 (41 units) were considered. Discharge data (Fig. 3b–d) in the form of annual maxima were derived by the Crati, Coscile, and Esaro stations (see Fig. 1a for position), recorded between 1927 and 1966 (31 units) for the Crati station, 1928 and 1978 (29 units) for the Coscile station, and 1929 and 1969 (25 units) for the Esaro station. The 5 × 5 m DEM was provided by the Centro Cartografico Regione Calabria and the drainage network by the Agenzia Regionale per la Protezione dell'Ambiente (ARPACAL) website (https://www.cfd.calabria.it/); the land use map was derived by Corine Land Cover-2018 (EEA 2020).

Fig. 3
figure 3

Annual maximum rainfall for given durations of 1 h, 3 h, 6 h, 8 h, 12 h and 24 h starting from 1962 to 2012, derived by the Torano Scalo rain gauge (a); discharge time series recorded by the hydrometric stations at the Crati River from 1927 to 1966 (b), Coscile River from 1928 to 1978 (c), and Esaro River from 1929 to 1969 (d)

Full size image

Probabilistic analyses were conducted on the annual maxima of rainfall time series by using the inverse Gumbel distribution function, with the aim of estimating the IDF curves, while the peak discharge values for specific return periods were estimated by conducting probabilistic analyses of the annual maxima of discharge time series by using the gamma distribution function. With the first analysis, the runoff–rainfall model considered (SCS-CN method) was able to generate runoff hydrographs on the basis of the parameters estimated by Eqs. 8–12: the runoff Q; the peak of discharge Qp; the time t intended as the duration of rainfall and derived by IDF curves; and the time to peak tp. To estimate the latter, it was necessary to calculate the concentration time tc, which, due to the large size of the analysed basin, required the application of multiple equations to derive a mean value. The relationships used were SCS-CN (SCS 1956), Ventura (1905), Giandotti (1934), Pezzoli (Indelicato 1988), Puglisi (Puglisi and Zanframundo 1978), and Tournon (Merlo 1973).

The contributing areas considered were obtained on the basis of specific thresholds of the drainage areas. We considered three contributing areas corresponding to (A) the entire extension of the Crati Basin; (B) an estimated contributing area of 1016 km2 based on a threshold of the drainage area of 1000 km2; and (C) an estimated contributing area of 286 km2 based on a threshold of the drainage area of 100 km2. The geomorphological parameters (i.e. extent, average slope, maximum, average and minimum elevation, length, and average slope) of the contributing area selected are reported in Table 1. Table 2 reports the values of tc estimated by different relationships for each contributing area. The calculated mean value was finally used to calculate the tp.

Table 1 Contributing area characteristics
Full size table
Table 2 Concentration time tc in hours for the contributing areas A, B, and C estimated by using different empirical equations
Full size table

Regarding the estimation of CN, considering different extensions of the contributing areas, a unique value was derived for each of these areas on the basis of the land use map (Fig. 4). For each of them, the HSG was the B group, and AMC was assumed to be class II due to the lack of information for the study site. Table 3 shows the CN values attributed to different land use areas in the Crati River Basin. The weighted average value of CN used for the analysis was estimated to be ~ 73 for all three contributing areas.

Fig. 4
figure 4

Land use maps for the contributing areas: a A, b B and c C

Full size image
Table 3 CN values of contributing areas A, B and C
Full size table

Finally, for modelling purposes, the 2D hydrodynamic simulations were set in the HEC-RAS model taking into account a computational cell size of 20 × 20 m considering the computational power required for the modelling. The boundary conditions were defined considering one inflow cell located within the Crati River on the upper boundary of the computational area and outflow cells located along the coastline and the upper and lower boundary of the computational area (Fig. 5).

Fig. 5
figure 5

Computational scheme of the hydrodynamic modelling. The figure shows the boundary of the computational area (red line), inflow cell and outflow cells (blue dot and blue line), and Sybaris archaeological area (black square)

Full size image

The channel’s Manning’s roughness coefficient n was set at 0.03 m−1/3 × s according to the type of channel simulated (natural channel; Chow 1959). The hydrodynamic simulation was run considering the developed 47-h (area A)-, 40-h (area B)-, and 30-h (area C)-long hydrographs over a total simulation time of 120 h. The outputs generated inundation maps in terms of the flow depth and extent for each return period considered. FH maps were obtained by combining the inundation maps for the poorest inundation condition of each contributing area, from the lowest layer corresponding to a 300-year event to the highest layer corresponding to a 10-year event. Since the analysis was completed through a reduced complexity approach, the flow velocity was not considered in the FH estimation.

5 Results and discussion

Figure 6 shows IDF curves obtained for the selected return periods, indicating the maxima intensities of 59.3 mm/h, 52 mm/h, 47.4 mm/h, and 36.4 mm/h for return periods of 300 years, 100 years, 50 years, and 10 years, respectively, with a critical duration of 1 h. Conversely, the minimum values, which were identified for a critical duration of 24 h, were 6.9 mm/h, 6.1 mm/h, 5.6 mm/h, and 4.3 mm/h for return periods of 300 years, 100 years, 50 years, and 10 years, respectively.

Fig. 6
figure 6

IDF curves for return periods of 10, 50, 100, and 300 years

Full size image

Considering the estimated runoff, rainfall implemented in the SCS-CN method was identified as the amount corresponding to critical duration, which in this case was equal to the time of concentration tc for each contributing area. For the mean tc estimated at 9.3 h, 7.2 h, and 5.9 h for contributing areas A, B, and C, respectively (Table 2), Table 4 provides the amount of rainfall calculated on the basis of these critical durations (i.e. the tc). This amount represented the term P in the Eqs. 8 and 9.

Table 4 Cumulative rainfall height (in mm) estimated for different tc (h) and return periods (year) related to contributing areas A, B, and C
Full size table

By estimating Qp and tp from Eqs. 11 and 12, respectively, synthetic unit hydrographs and runoff hydrographs were reconstructed for each contributing area and selected return periods. Figure 7 shows the runoff hydrographs for contributing areas A (Fig. 7a), B (Fig. 7b), and C (Fig. 7c). Peaks of discharge appeared to shift towards the left (from Fig. 7a to Fig. 7c) because of the differences in the extent of the contributing areas. This finding implied the durations of the hydrographs (47 h, 40 h, and 30 h for areas A, B, and C, respectively).

Fig. 7
figure 7

Runoff hydrographs estimated for different return periods (10, 50, 100 and 300 years) and for contributing areas A (a), B (b), and C (c)

Full size image

Considering contributing area A, peaks of discharge were reached after 12 h and corresponded to 907.7 m3/s, 1524.7 m3/s, 1808.4 m3/s, and 2278.9 m3/s for 10-year, 50-year, 100-year, and 300-year events, respectively. For contributing area B, peak discharges were reached after 8 h and assumed values of 453.5 m3/s, 779.6 m3/s, 930.8 m3/s, and 1182.8 m3/s for 10-year, 50-year, 100-year, and 300-year events, respectively. Finally, peak discharges were reached after 6 h for contributing area C, with values corresponding to 145 m3/s, 254.3 m3/s, 305.3 m3/s, and 390.7 m3/s for 10-year, 50-year, 100-year, and 300-year events, respectively.

The results in terms of inundation maps depicting the extent and depth of water flow are shown in Fig. 8. As expected, the worst inundation situation occurred, often shortly, after the time of the peak of the hydrographs (≥ 2 h later) and, in particular, after 14 h, 10 h, and 7 h for contributing areas A, B, and C, respectively.

Fig. 8
figure 8

Inundation maps. a–d for contributing area A, e–h for contributing area B and i–l for contributing area C for a 10-year, 50-year, 100-year, and 300-year events, respectively

Full size image

Figure 8a–d refers to contributing area A. The simulations result in the largest flooded area with no significant differences in the sizes of the inundation areas for very low to high-return periods. The Sybaris archaeological area is always affected by river overflow with water depths between a few centimetres and several metres. The Parco del Cavallo site has a flow depth of approximately 4 m for low-return periods, while Stombi, despite being located approximately 1.5 km away from the river, has a flow depth between 2 (10-year; Fig. 8a) and 4 m (300-year; Fig. 8d). Inundation maps for contributing area B (Fig. 8e–h) indicate a smaller flow extent than the previous case for low-return periods. The depth of flow reached 3 m at some parts of the Parco del Cavallo and Prolungamento Strada sites (Fig. 8e–h), while Casa Bianca and Stombi were less involved by the flow. In particular, the Stombi was not flooded for a 10-year event (Fig. 8e), while the maximum depth (almost 3 m) was observed for a 500-year event (Fig. 8h). Finally, Fig. 8i–l indicates that when considering contributing area C, the flow extent and depth in the archaeological sectors greatly reduces. For both 10-year and 50-year events (Fig. 8i, j), no significant overflows were observed; for a 100-year event (Fig. 8k), a limited magnitude flood (a few centimetres in depth) involved the site nearest to the river; and for a 300-year event (Fig. 8l), the flow depth reached approximately 3 m at the Parco del Cavallo and Prolungamento Strada and approximately 2 m at Casa Bianca. Otherwise, in this case, the Stombi site was never involved in flooding. The total inundation area showed a distinct decreasing trend as the return periods of the flow hydrographs decreased from contributing areas A to C.

Overlaying the inundation maps from the 300-year return period event, corresponding to the reference extreme event, to the 10-year return period event, the FH maps were assembled and classified on the basis of the corresponding return period (Fig. 9). Such a classification identify areas prone to flooding by events of specific return periods and represent a useful tool for supporting the future adoption of flood mitigation measures and the development of flood insurance plans (e.g. Linnenluecke and Griffiths 2010). In detail, the results show for contributing area A (Fig. 9a) the wider extension of the floodable area that involve the entire archaeological area, even for low-return periods. Conversely, the FH map for contributing area C (Fig. 9c) appears limited to some parts of the study area, including Parco del Cavallo, Prolungamento Strada, and Casa Bianca, for events with high-return periods (100-year and 300-year event). Even if their involvement in flooding processes depends on both the extent of the contributing area and the as-expected return period, the modelling results did not consider the long-term processes in the floodplain, such as subsidence. Subsidence is responsible for a lowering of ~ 20 mm/yr in the plain due to the neotectonic and glacio-eustatic changes and, above all, the compression of sediments (Cianflone et al. 2015). This and other external factors, such as the effects of climate change and channel adjustment, that are not considered in the FH assessment, could influence the reliability of the analysis (Whitfield 2012; Scorpio and Rosskopf 2016; Nicholls et al. 2021) and must be considered, especially in simulations aimed at investigating long-term effects.

Fig. 9
figure 9

FH maps for the three contributing areas considered: a A, b B, c C; classified on the basis of the selected return period (10, 50, 100 and 300 years)

Full size image

The main advantage related to the proposed procedure is linked to the multiscenario approach. Indeed, considering a well-defined contributing area, it was possible to reproduce the effects of localized rainfall and evaluate the potential for events with defined return periods considering only the features of the contributing area. In contrast, the common approaches consider a procedure that transformed rainfall into flood discharge (e.g. rainfall–runoff model) for which the critical rainfall is generally interpolated to cover the entire analysed basin. These methods aim to spread rainfall over the basin considering a uniformly distributed or concentrated event (e.g. areal reduction factor; Bell 1976) and only subsequently account for basin characteristics. Consequently, in the case of a large basin, these approaches could overestimate the flood discharge and consequently hypothetical infrastructures or lead to an underestimation of them related to the hydrological and morphological characteristics of the basin.

However, there are fundamental challenges in the adoption of a multiscenario approach. One of these concerns the variability of the extent of contributing area. Indeed, in this study case, the contributing area variability was estimated by a flow-direction model based on the morphological characteristics of the basin, but it could differ if derived based on other criteria, such as a stochastic simulation or radar observations of rainfall coverage (if available).

6 Validation of estimated hydrographs

To validate the first part of the analysis results, a comparison between the runoff hydrographs estimated from the SCS-CN method and the values of the peak discharge estimated through a probabilistic analysis from the available discharge time series was conducted. This choice was related to the absence of reliable data (i.e. historical data, geomorphological field observations, or remote sensing data) indicating the area flooded in the past. Although some recent events (see 2013 and 2018 events) caused damage to the Sybaris archaeological site, only incomplete and/or large-scale reports were available (e.g. Trigila et al. 2021). Since three annual discharge time series were available, the probability of exceedance of each specific discharge was fit by the gamma distribution function (Eqs. 6, 7). By such analyses, we calculated the theoretical return period as the inverse of the probability of exceedance. Then, to describe the uncertainty in the estimation of the probability of exceedance of discharge data linked to the features of time series (i.e. sample size and continuity), 95% confidence intervals were derived on the basis of the estimated standard errors for each data series (e.g. Guerriero et al. 2020a).

Figure 10 shows gamma probability models set based on discharge data (bars) measured at the Crati (Fig. 10a), Coscile (Fig. 10b), and Esaro (Fig. 10c) hydrometric stations. Grey dots described the probability of exceedance of discharge data series; continuous lines indicated the probability model for a discharge interval defined on the basis of minimum and maximum of discharge series; overlapped triangles identified the probabilities of exceedance corresponding to return periods of 10 years, 50 years, 100 years, and 300 years; and 95% confidence intervals are denoted by the dashed lines.

Fig. 10
figure 10

Gamma probability models set on the basis of discharge data measured at the Crati (a), Coscile (b), and Esaro (c) hydrometric stations. Solid lines indicate the probability model; triangles indicate the probabilities of exceedance corresponding to return periods of 10 years, 50 years, 100 years, and 300 years; and dashed lines indicate 95% CIs and provide an overview of the uncertainty associated with estimated discharge

Full size image

Discharges measured for each theoretical return period are shown in Table 5. The peak values of runoff hydrographs estimated for contributing areas A, B, and C are also reported. Even if considerable uncertainty affects the estimated discharge data, especially for high-return period events, as indicated by the 95% confidence intervals (i.e. dashed lines of Fig. 10), a qualitative comparison was attempted on the basis of the return period.

Table 5 Comparison between the discharge data derived by probability analyses of time series for the Crati, Coscile and Esaro rivers and peak discharge values for contributing areas A, B, and C derived by the runoff hydrograph for fixed return time events
Full size table

The results showed that the discharge of the Crati River is comparable to the peak value modelled for the hydrographs of contributing area A, especially for low-return periods; due to the extension of contributing area B, the peak discharge could be associated with the sum of the discharges of the Coscile and Esaro Rivers, while the discharge of the Coscile River matches the peak value modelled for the hydrographs of contributing area C.

7 Conclusions

In this paper, a procedure is proposed for estimating multiscenario flood hazards in large river basins, accounting for the potential variability levels of the contributing areas related to the occurrence of meteorological concentrated events or extending over a fraction of the river basin. This approach was tested at the archaeological site of Sybaris in southern Italy by first using the SCS-CN method to derive the runoff from the rainfall on the basis of land use and soil characteristics and subsequently modelling the runoff hydrographs over the flood plain through the 2D hydrologic engineering center river analysis system (HEC-RAS) model. The runoff hydrographs were derived by intensity–duration–frequency curves (IDF) obtained by a probabilistic analysis of rainfall data. Multiple contributing areas were considered to generate runoff hydrographs for different extents of contributing areas due to the unavailability of spatial rainstorm coverage for the basin. In the presence of such data, each scenario could correspond to a predetermined extent.

The applicability and expected results of the proposed procedure are related to (1) the characteristics of the rainfall time series (i.e. number and continuity), representing the basis for a reliable fitting of the probabilistic model; (2) the features of the contributing basin regarding the extension and river characteristics and the land use and hydrological soil conditions, on which the amount of rainfall contributing to flood generation depends (derived by the rainfall–runoff model) and related data, such as land use maps and digital elevation models; and (3) uncertainties linked to the hydrodynamic model derived by the accuracy of topographic data and the impossibility of measuring or directly estimating some parameters required by the model and the presence of validation data. However, the latter is a common issue in flood hazard analysis that can be reduced by adequate-quality hydrological input data. By considering the common availability of rainfall time history data, land cover maps, and digital elevation models, a wide applicability of the procedure is expected. A limitation of the procedure is related to the impossibility of quantitatively validating the inundation maps by historical data or geomorphological field observations. However, the same approach could be easily exported and tested in other fluvial contexts to analyse limited extension areas with available model data and historical inundation maps.