Flood Risk Evaluation in Ungauged Coastal Areas: The Case Study of Ippocampo (Southern Italy)

Abstract

The growing concentration of population and the related increase in human activities in coastal areas require numerical simulations to analyze the effects of flooding events that might occur in susceptible coastal areas in order to determine effective coastal management practices and safety measures to safeguard the inhabited coastal areas. The reliability of the analysis is dependent on the correct evaluation of key inputs such as return period of flooding events, vulnerability of exposed assets, and other risk factors (e.g., spatial distribution of elements at risk, their economic value, etc.). This paper defines a methodology to assess the effects of flooding events associated with basin run-off and storm surge in coastal areas. The assessment aims at quantifying in economic terms (e.g., loss of assets) the risk of coastal areas subject to flooding events. The methodology proposed in this paper was implemented to determine the areas subject to inundation on a coastal area in Southern Italy prone to hydrogeological instability and coastal inundation. A two-dimensional hydraulic model was adopted to simulate storm surges generated by severe sea storms coupled with intense rainfalls in order to determine the areas subject to inundation in the low-land area along the Adriatic coast object of this study. In conclusion, the economic risk corresponding to four different flooding scenarios was assessed by correlating the exceedance probability of each flooding scenario with the potential economic losses that might be realized in the inundated areas. The results of the assessment can inform decision-makers responsible for the deployment of risk mitigation measures.

1. Introduction

Urban development that occurred in the last decades is considered one of the key reasons determining an increase in occurrence and intensity of flood events with consequent social and economic damage in the affected areas [1]. Increase in flood occurrence has been recorded especially in small catchments where an increase in land use has had an impact on flooding areas [2,3,4]. In general, floodplains are areas characterized by the presence of stream channels developed by the combined effects of floods of variable scales and geomorphologic processes [5]. These areas are often mentioned as riparian regions or buffers [6,7] and are easily recognizable from adjacent areas due to variances in hydrological storage [8,9,10] and other environmental aspects [11]; therefore, an accurate flood mapping is particularly important to identify such areas. In fact, flood mapping is a crucial element of flood risk management. Furthermore, in small ungauged catchments, the lack of observed discharge data makes flood mapping very difficult as these data are required to calibrate hydrological and hydraulic models [12,13,14].

The increase in human activities due to global population growth requires the employment of numerical simulations to determine the effects of flooding events with the objective of informing risk assessments and management practices to be adopted in the affected coastal areas [15]. Coastal risk is strictly connected to erosion phenomena and marine inundation [16,17,18]. Several studies have been conducted to investigate coastal flooding areas in recent years [19,20].

With regards to the Adriatic Sea, for example, several coastal inundation maps have been produced [21,22,23]; however, most studies have only focused on storm surge contribution.

Nevertheless, even though flood risk has been studied for a long time with the development of numerous numerical approaches [24], several aspects have not yet been completely investigated such as the influence of contemporary event of river floods and marine inundation, which all require further analyses.

Other studies [25,26] focused on floodplains located in hurricane-prone areas subject to heavy rainfall and storm surge have been carried out. These employ coastal hydrodynamic models coupled with topography-based hydrologic methods. These studies adopt a high-resolution modeling approach which requires high spatial/temporal resolution data sensitive to natural hazard characteristics such as storm intensity, track, and topography. This approach results computationally demanding, thus can be hardly applied to large areas.

Ray (2011) [27] and Torres (2015) [28] analyzed the combined effects of storm surge and inland rainfall using HEC-RAS, a widely used flow model, and generated floodplain maps under hurricane scenarios.

Heavy rain and sea storms, in the context of global sea level rise, increase the probability of coastal flooding in areas such as the Mediterranean Sea [29,30]. Recently, both Adriatic and Ionian sides of Apulia were flooded for thousands of square meters due to sea level surges generated by severe storms in connection with intense rainfalls [31].

As a consequence of the recorded sea level rise and the increased probability of extreme events, in recent years, there has been a growing interest in determining the economic consequences floods and erosions risk along coastal areas [32,33,34].

In line with the latest development made in recent years in this area of study, the purpose of this paper is to provide a methodology to assess the economic risk associated with four flooding hazard scenarios generated in a low-lying area by storm surges coupled with intense rainfalls.

Existing flooding maps in the test area have been generated using only rainfall as input to hydrodynamic models without taking into account the contribution of storm surges [35]. However, an economic assessment of inundation risk has never been carried out.

FLO-2D model has already been used to simulate storm surge inundation [36] and marine ingression coupled with intense rainfall [37]. Unfortunately, due to the lack of validation data, these studies have only conducted a qualitative analysis. In contract with previous attempts, a past storm surge event was simulated using the FLO-2D model at a small spatial scale and validation was carried out with the aim of assessing the model accuracy in flooding prediction.

Hence, 1-, 30-, 50-, and 100-year return period predicted floodplains were used to perform an economic assessment. The economic consequences were determined by multiplying the economic value of elements at risk involved in each flooding scenario by their vulnerability, i.e., the percentage of damage associated with each typology of the element at risk [38,39]. The vulnerability was obtained through flood damage functions which are dependent on the damage of each element at risk with flooding water levels and flow velocity [40,41].

In the present paper, a “combined hydraulic modeling” is presented, whereas, in most studies, the risk assessment in a coastal area takes into account only storm surge contribution, even though, in low lying areas with the presence of water channels, the classical approach could provide an underestimation of economic risk (as shown in the “Results and Discussion” section). Moreover, the FLO-2D model output was validated by employing real data, proving the model to be an effective tool in identifying coastal areas exposed to storm surge inundation. In addition, it is very important to highlight that the geographical object of this study, despite being a flood-prone area, has not been the subject of similar studies in the past and therefore this study improves the knowledge around the identification of the most vulnerable areas at a local scale under the guidelines provided by the 2007/60 Flood Directive.

The paper is structured as follows. Section 2 illustrates the features of the case study. Section 3 provides details on the numerical model adopted. Section 4 provides a summary of the results obtained. Section 5 provides the conclusions of the case study.

2. The Study Area

The examined area is located within the Gulf of Manfredonia in the Adriatic Sea (Apulia Region, Southern Italy) and the coastline is oriented north–south with a narrow sandy beach backed by salt marshes and cultivated crops.

The study area is crossed by an important hydrographic network and is bounded on the north by the Cervaro River, on the south by the Peluso canal and the Carapelle River, and on the east by the Adriatic Sea (as shown in Figure 1).

The rainfall regime is Mediterranean with an average rainfall of 442.1 mm per year and an average of 63.4 rainy days per year. The maximum rainfall occurs in the autumn–winter period, with a strong dryness during the summer months [42].

Elevation is lower than 15 m above mean sea level and the area also includes a large part of the ancient Salso Lake. This area has historically been mostly uninhabited and wild, being naturally rich in marshes and swamps, and it is classified as a SIC (Site of Community Importance for the European Commission Habitats Directive—92/43/EEC) [43] (SIC IT9110005 “Capitanata wetlands”).

Starting from the 1970s, some local landowners gave impetus to an extensive reclamation of the area, which had already begun since the early twentieth century, transforming a highly natural landscape into intensively cultivated agricultural land and built-up areas. A large network of reclamation channels, currently in a poor state of conservation, was specifically built. The expansion of the cultivated land, to produce vegetables close to the coastline, destroyed the coastal dune cordon reducing it to a narrow embankment to protect crops. In that period, the Villaggio Ippocampo tourist complex was also built and further developed in the following years and nowadays in the summer period, there are more than 15,000 residents.

Since the mid-20th century, the long sandy beaches of the Gulf of Manfredonia have suffered a significative shoreline retreat due to a massive anthropic action [44]. At present, the retreating coastline extends from Margherita di Savoia to the Ippocampo tourist village and the erosive process is progressing towards the port of Manfredonia [31]. Due to the hundreds of protection structures built in the past 40 years, the examined area has especially been observed in the last decades, becoming a “case study” [18].

In addition to coastal erosion issues, the study area is also aggravated by flooding hazard [35] because, due to the presence of low-lying and depressed areas, coastal flooding occurs also during ordinary storms. Furthermore, sea storms are often associated with significant meteoric events with consequent river overflow. However, it should be noted that, despite the twentieth-century anthropization, half of the territory to be examined can be classified as “natural area”. About 50% of the natural areas are listed as habitats indicated in the Habitats Directive 92/43/EEC [43].

Trying to analyze the land-use, it is evident that the area, as already mentioned, is characterized by a strong anthropization, which totally transformed, in the last decades, the natural environments. Specifically, the construction of the “Ippocampo” tourist complex took place on brackish moist soils characterized by halophilous vegetation, which was recognized by the European Community as a habitat of community importance. The phenomenon of anthropization has also affected the beaches and the entire dune system, both for the construction of the shores and the cultivation of the land. Together with the total obliteration of the dune cordon to make way for the bathing structures, various degrees of transformations of coastal dunes are found, ranging from the construction of an embankment to protect the internal spaces with the planting of exotic species (e.g., Carpobrotus acinaciformis), to the reduction of the original dune and its typical dune vegetation (e.g., Ammophila sp.). These data highlight the importance of the examined area from the point of view of environmental conservation.

3. Materials and Methods

3.1. Numerical Model Description

To assess the coastal inundation risk zones in the case study area along the Adriatic coast (Figure 1), two-dimensional hydraulic modeling software (FLO-2D software) has been used to simulate a storm surge approaching a low sandy beach added to river flooding.

The FLO-2D model [45] is generally used for different types of hydraulic tests such as the propagation of floods and debris flow events. It allows us to simulate flooding on a complex topography and a given roughness based on the conservation of volume. O’Brien [36] also simulated ocean storm surge scenarios occurring in Hawaii by mapping the inundation areas. To simulate inland flooding, FLO-2D requires the topography of the area, which can be represented by Digital Terrain Model (DTM) points, roughness coefficients, and input hydrographs. The FLO-2D model produces time step flow depth and velocity maps and the maximum flow depth map outlining the maximum area of inundation.

The entire calculation domain extends 65.5 km² and consists of a square mesh grid with a resolution of 15 m. In the model, the rainfall input has been entered in the form of hydrograph (discharge versus time), whereas the storm surge has been treated as a boundary element with a stage–time relationship.

3.2. Rough Data Collection

The fundamental characteristics of the climate were derived by analyzing and processing the data for the fifty years 1951–2000 collected by the Hydrological Annals (Hydrographic Service of the Ministry of Public Works) [42], referring to the “Manfredonia” thermopluviometric station (2 m a.s.l.) a few kilometers from the study area.

The sea level analysis was carried out exploiting tidal data acquired from 2007 to 2015 by the tide station located in the port of Manfredonia belonging to Apulia Region Meteomarine Network (also known as SIMOP) [46]. That regional meteomarine network has been operating since 2006, collecting wave, wind, and tidal data along the Apulian coastline [47].

The Manfredonia tide station, in particular, was deployed because there was a lack of data exactly in an area at high risk for coastal flooding [21]. The tidal station is equipped with an ultrasonic gauge incorporating an internal temperature sensor to compensate for temperature changes. Identification of anomalous data, spikes, and timing errors has been carried out following quality control checks [48], and spurious records have been removed from the time series. The amplitude of tidal oscillations recorded in the Gulf of Manfredonia during the observation period is of the order of tens of centimeters and is compatible with sea level variations in the Mediterranean Sea.

The analysis of sea level data showed an increase in the mean sea level and an increase in tidal peaks starting from 2009, with a maximum recorded on 6 February 2015 equal to +0.70 m a.s.l. [31]. A sudden increase in mean sea level has been reported from winter 2008–2009, with a variation between 2007 and 2009 of about 10 cm, followed by a decrease in 2011, while no significant changes have been found in meteorological tides [49].

For the examined site, offshore observed wave data are not available, hence we referred to a wave time series collected about 130 km south of Ippocampo by a directional wave buoy moored offshore Monopoli in the Adriatic Sea belonging to the Italian National Buoy Network.

The observed wave dataset covers an 18-year period, from 1990–2007, and consists of significant wave height (H_m0), peak wave period (T_p), and mean wave direction (D_m).

Offshore depth information was taken from the nautical maps produced by the Italian Navy Hydrographic Office (IIM) and a specific field campaign was carried out using a single-beam echo sounder associated with a GPS/RTK positioning system. Collected data were used to create a detailed bathymetric map up to a depth of about −20 m (Figure 2).

A high-resolution Digital Terrain Model (DTM) generated from LiDAR data covers the area and it has been further improved by additional topographic data acquired during a field survey. LiDAR data were used to create a structured computational grid with a cell size of 15 m. The entire calculation domain extends 65.5 km² and the grid resolution was chosen to reduce the computational cost.

Flow resistance or “roughness” is represented by Manning’s friction coefficient “n” (s/m^1/3) and, for floodplains, it depends on the substrate and vegetation density [50,51]. The surface roughness, in each simulation, was assumed equal to the default value in the model (0.04 s/m^1/3) for all cells inside the computational domain.

3.3. Hydrologic Modeling

Following a consolidated technical practice in the regional area, the VA.PI. Apulia method [52] provided the estimation of the precipitation input.

To this purpose, the intensity–duration–frequency curves of precipitation (IDF) were defined. Each IDF curve is associated with both a scenario and a hydrographic sub-catchment. The study considered four return periods (T_R = 1 year, T_R = 30 years, T_R = 50 years, and T_R = 100 years) corresponding to four different scenarios. To suitably define the floodplain, since the whole hydrographic network included in the study area is composed of several river branches, it was necessary to identify a hydrological input for each of them. For this purpose, the overall domain was divided into sub-catchments, each pertaining to a stretch of the hydrographic network.

A monomial exponential expression derived from two-parameter law described the rainfall IDF curves:

$$i_m(t, T_R)=a{t}^{n-1}$$

(1)

where i_m is the average intensity of the precipitation, a function of the duration of the rainfall event (t) and the return period (T_R); a is a parameter depending on the return period (T_R) and considered area; and n is a parameter depending on the considered area.

A regional analysis of the annual maximum rainfall intensity, performed with a probabilistic model based on the use of a Two Component Extreme Value (TCEV) Distribution [53], maximum likelihood estimator, and hierarchical estimation of regional model parameters [54], allowed the estimation of the rainfall IDF curves.

The regional analysis was based on the VA.PI. method, provided by the National Group for the Defense from Hydrogeological Disasters [52], through the probabilistic analysis of annual maximum rainfall of assigned duration (1, 3, 6, 12, and 24 h) [42].

The VA.PI. method is a national report, articulated through a series of “regional reports”, that presents a synthesis of the studies carried out in different areas of the national territory. Each regional report explains how to use the entire procedure developed at the local level.

The regional report for the Apulia Region defines an effective procedure [52]. Exactly, through a “regionalization” process, homogeneous zones for rainfall were identified within the territory of the Apulia Region. An equation was associated with each zone to evaluate the relative IDF curve, as a function of the rainfall depth (expressed in mm).

The whole area examined corresponds to the relationship:

$$h (t, z)=22.23t^{0.247}$$

(2)

where h represents the rainfall depth; t is the time of the precipitation storm; and z is the mean elevation of the basin. In this case, Zone 2 is not dependent on the parameter z.

Studies by Eagleson [55] and Penta [56] made it possible to apply the procedure to the basins of Southern Italy.

Moreover, the growth factor K_T allowed pondering the IDF curve dependence from the return period (T_R). Statistical studies suggest the growth factor is constant for the whole northern Apulia Region [52], individuated within the previously mentioned VA.PI. Apulian Regional Report.

In particular, for these zones, the VA.PI. report defines the following equation to rate the growth factor value:

$$K_T=0.5648+0.415\ln T_R$$

(3)

For each duration of rainfall, values of corrected rainfall depth (h_c) are estimated:

$$h_c=h (t,z) K_T(T_R)$$

(4)

The following relationship couples the corrected rainfall depth (h_c) with the mean intensity of rainfall:

$$h_c=i_{m}t=a{t}^n$$

(5)

Thus, the IDF curve can also be expressed, in terms of rainfall depth, as:

$$i_m(t, T_R)=a{t}^{(n-1)}$$

(6)

The IDF curves led to the determination of the peak flow rates of the basins hydrographs, which is an input of the flooding modeling through the runoff curve number method proposed by the USDA Soil Conservation Service Agency (SCS-CN) [57].

This method is based on the use of the empirical parameter called Curve Number (CN), widely used in hydrology to estimate the surface runoff from the mass rainfall. The CN parameter depends on the soil type, land use, treatment, and moisture condition [58], and it can vary between 0 and 100. Following the indications of the local basin authority contained in the PAI Puglia [35], in this case, all values of CN were referred to a wet condition of the soil (AMC-III), for the benefit of security.

For each sub-basin considered in the computational domain, floods with a semi-distributed model were estimated, based on the average height, slope, surface, length of mainstream, and CNIII values of the relative basin [59].

Table 1. CN-III and peak flow rates summary.

Sub-Catchment	CN(III)	Q Peak (m3/s)
		TR = 1 year	TR = 30 years	TR = 50 years	TR = 100 years
TRIBUTARY 1	87	11.08	29.64	34.75	41.70
TRIBUTARY 2	86	1.05	13.26	15.48	18.56
TRIBUTARY 3	88	0.83	9.39	10.92	13.04
CHANNEL 1	88	0.59	7.51	8.76	10.51
CHANNEL 2	90	0.53	4.66	5.36	6.33
CHANNEL 3	90	0.66	5.33	6.12	7.19
CHANNEL 4 UPSTREAM	90	1.02	12.19	14.20	16.99
CHANNEL 4 DOWNSTREAM	87	0.32	2.97	3.43	4.06
PELUSO 1	89	20.70	55.35	64.90	77.87
PELUSO 2	92	0.13	1.24	1.43	1.69
PELUSO 3	92	0.09	0.95	1.11	1.32
PELUSO 4	89	0.02	0.29	0.33	0.40

reports the values of CN-III and the peak flow rate in m³/s for each sub-catchment.

The shape adopted for the hydrograph is triangular, with the discharge peak at the time of concentration (t_c) of each sub-catchment.

The time of concentration (t_c) parameter was estimated through the relationship:

$$t_c=\frac{t_l}{0.6}$$

and the Mockus formula, provided in [60], to evaluate the Lag Time (t_l)

$$t_l=0.342\frac{L^{0.8}}{s^{0.5}}{(\frac{1000}{\mathrm{CN}}-9)}^{0.7}$$

(7)

where L is the maximum flow length, expressed in kilometers; s is the mean slope; and CN is the value of the Curve Number parameter.

Fiorentino [61] conceived a method for estimating flood volumes based on the use of flood reduction curves, which led to the assessment of the hydrograph’s overall duration.

The so-called flood reduction curve expresses the relationships between the flood peaks and the volumes passing through. It represents the ratio between the index flow rate on the variable duration and the peak flow rate for a given return period (T_R), described by the following relationship:

$${\epsilon }_{t,T_R}=\frac{Q_{t,T_R}}{Q_{T_R}}$$

(8)

where ${\epsilon }_{t,T_R}$ is the coefficient of flood reduction, depending on the duration and the return period (T_R); $Q_{t,T_R}$ is the index flow rate, depending on the duration and the return period (T_R); and $Q_{T_R}$ is the peak flow rate, depending on the return period (T_R).

Bacchi [62] explained that the dependence on the return period, under suitable hypotheses on the frequency distribution of the flows, can be neglected. Thus, the reduction curve can also be defined as the ratio between the average of the extremes of a given duration $(\mathsf{\mu }(Q_t))$ and the index peak flow rate $(\mathsf{\mu }(Q))$.

$${\epsilon }_t=\frac{\mathsf{\mu }(Q_t)}{\mathsf{\mu }(Q)}$$

(9)

Conceptual approaches have been developed to allow an application even in ungauged basins. Several authors have proposed formulas in which the various parameters of the reduction curve are directly related to the geomorphoclimatic characteristics of the basin, such as the basin lag time (t_l) and the coefficient n, exponent of the two parameter IDF curve.

In 1985, Fiorentino [61] proposed a formulation of the reduction curve, based on the runoff model’s conceptual hypothesis of a linear reservoir:

$${\epsilon }_t=(\frac{k}{t})(1-e^{-\frac{t}{K}})$$

(10)

It is a monoparameter formulation, in which the parameter k [63] is dependent on the lag time (t_l) and the IDF exponent n:

$$k=1.027 t_l e^{2.277 n}$$

(11)

Once the k value for each sub-catchment is established, the assessment of the flood volume ($V_{t,T_R})$ is then possible using the following expression:

$$V_{t,T_R}=Q_{t,T_R}\times t$$

(12)

Replacing the expression of ${\epsilon }_{t,T_R}$ in Equation (12):

$$V_{t,T_R}={\epsilon }_{t,T_R}\times Q_{T_R}\times t$$

(13)

The overall duration of the input hydrograph has been established to obtain a flood volume ($V_{ input T_R}$), for each return period considered, higher than that ($V_{T_R}$) calculated by maximizing the Equation (13).

3.4. Coastal Modeling

To investigate coastal flooding, the specific input data required for the FLO-2D model are water levels at the shoreline, as a function of time, calculated as the sum of storm surge, astronomical tide, wave set-up, and wave run-up.

Storm surge is defined as the sum of the effects that barometric pressure, wind stress, Coriolis force, and wave breaking have on sea level. The sea surface elevation recorded by tide gauges is the sum of astronomical and meteorological tide induced by wind, storms, and atmospheric pressure disturbances, which can be split up performing harmonic analysis [64,65]. The meteorological tide can be predicted by implementing the extremal analysis of tidal residuals to obtain events with a given return period.

The 1-, 30-, 50-, and 100-year return period storm surges were estimated in the Gulf of Manfredonia exploiting a hindcast tide dataset calculated using a numerical model Regional Ocean Modeling System (ROMS) [31].

To investigate wave-induced sea level set-up, the available wave dataset observed offshore Monopoli was transferred offshore Ippocampo by using the geographical transposition of wave gauge data [66] based on the empirical Sverdrup–Munk–Bretschneider (SMB) model [67]. Estimated wave data were analyzed and the distribution of percentage frequency of wave direction and height was analyzed. Wave climate offshore in the Gulf of Manfredonia is moderate and the highest sea storms approach from ENE (Figure 3).

In addition, the extremal wave analysis was carried out fitting all the waves coming from homogeneous sectors [68] above a threshold value (peak over threshold method [69]) to a Weibull extreme-value distribution (

Table 2. The 1-, 30-, 50-, and 100-year time return wave height and water level.

Sea Level	TR = 1 year	TR = 30 years	TR = 50 years	TR = 100 years
Hm0 (m)	3.4	5.63	6.03	6.04
R2% (m)	0.68	1.12	1.20	1.27
Tidal residual (m)	0.28	0.36	0.38	0.4
Astronomical tide (m)	0.2	0.2	0.2	0.2
Total sea level (m)	1.16	1.68	1.78	1.87

).

Extremal waves were propagated towards the coast using the Simulating WAves Nearshore (SWAN) model [70,71] using a nesting approach to improve the resolution in the coastal area. A high-resolution 10-m grid, extending from shoreline to −10 m contour, was nested in a coarse grid (at 200-m resolution) covering a large area in the Gulf of Manfredonia.

The SWAN model was run to estimate wave heights and breaker depths using the values of 1-, 30-, 50-, and 100-year return period waves.

Given the breaking wave characteristics, it is possible to evaluate wave run-up, defined as the maximum vertical extent of wave uprush on a beach or structure above the still water level (SWL).

Several formulas for predicting wave run-up have been developed from experimental studies [72,73,74,75]; in this work, the Stockdon formula was used [76]:

$$R_{u2\%}=1.1 \left(0.35 {\beta }_f {(H_0{L}_0)}^{0.5}+{\left[\frac{H_0{L}_0(0.563 {\beta }_f^2+0.004)}{2}\right]}^{0.5}\right)$$

(14)

where R_u2% is the run-up level exceeded by two per cent of the incoming waves, H₀ is offshore wave height, L₀ is the offshore wavelength, and bf is the foreshore slope. The beach slope was evaluated by considering the mean slope of the vertical region whose height is two times the breaking depth [16]. As shown in Figure 2, a very slight slope (about 0.3%) can be detected along the sea bottom in the examined area.

The Stockdon equation already includes wave setup and requires offshore wave characteristics, thus the use of this formula is very convenient for cases of application to large areas and where bathymetric data are not available.

The sum of tidal levels, storm surge, and wave run-up provides the total water level, as reported in Table 2. For each selected scenario, the water level time variation was represented as a triangular shape with the peak value equal to the estimated total sea level [77].

3.5. Risk Assessment

Generally, for natural calamitous events (landslides, floods, etc.), the risk can be quantified as the product of hazard, vulnerability, and exposure [78]. Hazard represents the probability of occurrence of a specific hazard scenario with a given return period. Vulnerability is the degree of damage of the elements at risk resulting from the occurrence of a natural hazard of a given intensity [79]. Finally, exposure defines the system at risk (its components, elements, functions, activities, etc.) and its value (social, economic, functional, etc.).

The quantitative assessment of risk, in monetary terms, was carried out by computing the economic losses consequent to potential flood events for different return periods (1, 30, 50 and 100 years). Firstly, all elements at risk were recognized and typologically classified in the study area. In total, 20 types of exposed assets were identified: electric cabin, aqueduct cabin, swimming pool, sports area, cistern, shed, agricultural shed, building, tree planted area, arable crop, unqualified garden, nude soil, shrub pasture, Mediterranean shrub area, olive grove, orchard, vineyard, greenhouse, urban road, and extra-urban road. These elements at risk were subsequently grouped into five categories, namely buildings, other structures, infrastructures, specialized land use, and unspecialized land use, which are characterized by different damage functions.

Among the elements at risk, the population was not considered for different reasons. First, it represents an intangible loss [80]. Moreover, assessing the temporal and spatial distribution of people in an area is an objectively difficult goal and quantifying the economic value associated with injuries or deaths is a complex ethical task. Population is a mobile asset [81] and, for this reason, the evaluation of exposure of persons would require the calculation of the conditional probability of persons being present in buildings or on the roads at the occurrence time of natural hazards [78].

The procedure used for evaluating quantitatively both the exposure of elements at risk and the economic risk was derived from [38,39]. The exposure assessment was carried out by overlaying of the elements with each flooding hazard map and evaluating the amount of assets in each flooding areas and quantifying the economic values for the various flood hazard scenarios (Figure 4, Figure 5 and Figure 6). Each type of element at risk was quantified in terms of areal extent (square meters or hectares). The spatial distribution of assets was combined individually with the four flooding hazard maps to obtain the number of elements potentially affected by flooding with 1-, 30-, 50-, and 100-year return period.

The spatial distribution of elements at risk was derived from the Regional Technical Map of Apulia at 1:5000 scale, for structures and infrastructures, and from Land Use Map of Apulia at 1:5000 scale, for land use categories, both available at geo-portal of Apulia Region [82].

Therefore, to evaluate the economic exposure associated with the different flooding hazard levels, unit market values or unit construction costs were assumed for each element at risk category. Starting from buildings, the unit market value (Euros per square meter) relative to 2019 was obtained from the Observatory of Real Estate Market (OMI) instituted by the National Territorial Agency. OMI is a cadastral database and provides maximum and minimum values for different kinds of buildings, such as residential, commercial, agricultural, etc., and for different zones of the inhabited area, e.g., center, suburbs, etc. In the area occupied from the village, residential, mainly villas, and commercial buildings are present. To obtain an overall unit economic value, a weighted mean among the market values related to residential and commercial buildings was calculated. The weights were derived considering the percentage distribution of commercial and residential buildings (about 30% and 70%, respectively). For sports areas, electric cabin, aqueduct cabin, agricultural sheds, cisterns, and swimming pools, the unit construction costs (Euros per square meter or cubic meter) were calculated from values defined in 1988–1989 by the National Territorial Agency. These values were updated, by calculating the inflation rate using data of ISTAT (National Institute of Statistics). In addition, the unit agricultural values (Euros per hectare) were obtained from the National Territorial Agency. The most recent data for the study area (included in the Foggia Province) are related to 2012. Finally, for roads, unit maintenance costs obtained from the Regional Price List of Apulia for 2019 were used. This unit maintenance cost includes mainly cleaning road work, i.e., mud removal from the roadway, drainage ditches, inspection wells, and scarps.

The quantitative evaluation of risk associated with the four flood hazard scenarios, in terms of economic expected losses, was carried out starting from the exposure estimation, according to Figure 7 and Figure 8. In particular, the exposed values associated to each flooding scenario were quantified by computing the amount of elements at risk (square meters, cubic meters, and hectares) involved in each flooding scenario and, subsequently, multiplying the amount of the individual exposed assets with their unit economic value.

For each category of elements at flood risk, the potential damages were obtained by taking into account their vulnerability to flood damages, expressed by flood-damage functions. Velocity and flood depth values, resulting from the hydraulic modeling and varying for each flood hazard scenario, permitted calculating a vulnerability (Ve) value, ranging from 0 to 1 (

Table 3. Flood damage functions.

Categories of Elements at Risk	Vulnerability (Ve)
	h (m); v (m/s)
Buildings	Ve(h) = 0.5∙h
Other structures		0.5∙h	if v < 2
	Ve(h,v) =
		0.35∙h∙(1 + 0.25∙v)	if v ≥ 2
Specialized land use	Ve(h,v) =	h	if v < 0.25
		h∙(1 + v)	if v ≥ 0.25
Unspecialized land use	Ve(h,v) =	0.25∙h	if v < 0.25
		0.25∙h∙(1 + v)	if v ≥ 0.25
Infrastructures	Ve(h)—Damage function from [40,41]

). Therefore, Ve defines the percentage of damages of economic assets in function of physical characteristics of floodwaters. In fact, each category of the five above-mentioned was associated with a different vulnerability flood-damage function. The flood damage functions of buildings, other structures, specialized, and unspecialized land use were referred to [40], whereas that of infrastructures was referred to [41].

The consequences were assessed by multiplying the exposed economic value of assets by the Ve values of each category of elements at risk. For each flooding scenario, the total consequences (or expected losses) were calculated by summing the values related to the five categories of elements at risk. The economic losses values associated with each flooding scenario were multiplied with the corresponding exceedance probability of hazard scenarios to obtain the related risk in monetary terms.

4. Results and Discussion

The first set of simulations was run to check the reliability of the inundated areas and the accuracy of the simulated flooded areas had been checked by comparing the model outputs with the effects of an observed event.

In November 2009, a severe meteo-marine event struck the examined area, causing flooding and considerable damages. In the hours following the disaster, a quick field survey was carried out in the area to detect inundated regions and numerous photographs had been taken (Figure 5). Initially, following the common practices for the evaluation of the risk on the coastal areas [21,22,23], a FLO-2D model based only on storm surge input was defined exploiting the water levels and wave data collected in the area. However, it proved unable to fully return the observed floodplain, because only the areas closest to the coast were flooded.

Therefore, a combined model that included both storm surge and flood input was run and the accuracy of the simulated inundation areas was tested comparing the model results with the observed flooded areas. Quite good accordance was found, as shown in Figure 6, even if some differences, especially in the agricultural area, can be explained considering that observed inundation boundaries have been collected during an expeditious survey along the damaged roads in the days following the event. In the urban zone, where the post-event flooding area was detected with extreme detail, the numerical results perfectly overlap.

After the model results validation, the FLO-2D [45] model processed for each selected scenario (T_R = 1-, 30-, 50-, and 100-year return period) two flooding simulations, one considering only the storm surge contribution and one using the storm surge combined with the rainfall.

Storm surge and rainfall depth were modeled considering the same return period, since in the case of a small catchment, mostly located in the coastal areas, the same synoptic weather system may induce extreme sea level and heavy rainfall [83]. Moreover, Bevacqua et al. [84] and Moftakhari et al. [85] confirmed a substantial decrease in return periods if the joint occurrence probability of sea water level and river discharge is considered. Hence, in these cases, in favor of safety, it could be appropriate to consider the probability distribution of the extreme events of the two random variables substantially coincident, because the meteorological forcing is the same. Therefore, the presented combined model represents the worst-case scenario and could be a useful tool for coastal area management.

The simulation provided results in the form of flood maps in terms of flow depth and velocity. These maps presented zones affected by values of both flow depth and velocity that were negligible for the realistic evaluation of the connected hazard distribution. Technical literature called these areas “marginals”. The delimitation of the “marginal inundation areas” was carried out following the methodology suggested by the River Tevere Basin Authority [86].

The applied methodology conservatively considers a water depth above 0.2 m and a velocity flow greater than 0.3 m/s as the limits of danger, whereas the flooded areas with smaller values are defined as marginal areas. In particular, all water depth values >0.2 m were considered to identify flooding areas, regardless of velocity values. At the same time, all velocity values >0.3 m/s were considered, regardless of depth values.

Following this approach, all flooded areas resulting from the hydraulic simulations were mapped (see Figure 7 for the storm surge model and Figure 8 for the integrated model). To develop hydraulic modeling, FLO-2D is preferred, for which the numerical stability was preliminarily verified. This aspect made possible simultaneously modeling both the storm surge input and the flood hydrographs for the different return periods, in a more detailed calculation domain. The inundation contours simulated with the FLO-2D model (see Figure 8) highlight the presence of large areas at risk of marine ingression in all the investigated scenarios and the flooded areas increase further if the rainfall is included in the model.

The zone most affected by flooding includes built-up areas, wetlands, and agricultural land, constituting a potential hazard for people, facilities, and viability. The presence of channels and rivers heavily affects the result of the analysis. The water channels, in fact, constitute a preferential way of marine ingression directing water inland for hundreds of meters and in the case of combining with heavy rain, due to the increased water depth in the channel, the flooding event involves a greater area. The extension of flooded areas increases as the year return period increases, but it is important to highlight that even a one-year return period event will inundate the urban area. Consequently, this scenario was also considered in the risk assessment of the study area.

The output of flooding simulations, in terms of flooding hazard maps, was combined with the distribution in the study area of elements at risk to carry out the exposure analysis according to Figure 4. This procedure was applied for both flooding simulations, i.e., both considering only the storm surge contribution and using the storm surge combined with the rainfall (combined model). In Figure 9, the elements at risk, subdivided into structures, infrastructures, and land use, involved in flooding areas derived from the four combined models are highlighted.

Subsequently, for each of flooding scenarios and each category of elements at risk, the amount of exposed assets potentially involved in the flooding events and their economic values were computed. Therefore, the consequences were obtained by multiplying the economic values of the exposed assets with vulnerability values, deduced from the damage functions in Table 3. Finally, the total economic losses associated with each flooding scenario were computed by summing the consequence values of the five categories of elements at risk. The results of exposure and risk assessment are synthesized in

Table 4. Amount and economic value of exposed assets and overall consequences, for each element at risk category and flooding scenario, considering only the effect of storm surge.

Elements at Risk	TR	Exceedance Probability	Amount (sqm)	Economic Value	Consequences	Total Economic Losses
Buildings	1	1	163	€891,831	€163,634
Other structures	1	1	365	€35,470	€4661
Road	1	1	36,762	€1,028,209	€207,380
Specialized land use	1	1	1903	€4563	€1783
Unspecialized land use	1	1	585,449	€315,691	€34,665	€412,123
Buildings	30	0.033	1864	€3,432,475	€633,715
Other structures	30	0.033	632	€61,271	€13,060
Road	30	0.033	57,393	€1,604,410	€414,963
Specialized land use	30	0.033	1903	€4563	€4080
Unspecialized land use	30	0.033	728,335	€421,024	€84,052	€1,149,870
Buildings	50	0.02	2234	€4,379,509	€774,118
Other structures	50	0.02	632	€61,271	€13,640
Road	50	0.02	59,299	€1,657,620	€435,358
Specialized land use	50	0.02	1925	€4616	€4212
Unspecialized land use	50	0.02	739,636	€427,151	€88,295	€1,315,623
Buildings	100	0.01	2407	€4,721,521	€880,093
Other structures	100	0.01	632	€61,271	€14,794
Road	100	0.01	64,811	€1,811,336	€474,252
Specialized land use	100	0.01	1,925.01	€4616	€4319
Unspecialized land use	100	0.01	759,036	€441,781	€96,043	€1,469,500

and

Table 5. Amount and economic value of exposed assets and overall consequences, for each element at risk category and flooding scenario, considering the combined effect of storm surge and rainfall.

Elements at Risk	TR	Exceedance Probability	Amount (sqm)	Economic Value	Consequences	Total Economic Losses
Buildings	1	1	917	€1,777,765	€271,608
Other structures	1	1	365	€31,210	€3689
Road	1	1	98,534	€2,703,967	€517,629
Specialized land use	1	1	56,000	€127,087	€49,089
Unspecialized land use	1	1	4,050,000	€6,241,545	€466,578	€1,308,593
Buildings	30	0.033	5967	€10,291,706	€1,867,053
Other structures	30	0.033	633	€61,568	€13,090
Road	30	0.033	270,945	€7,408,776	€1,592,121
Specialized land use	30	0.033	316,000	€559,501	€251,005
Unspecialized land use	30	0.033	12,780,000	€21,062,380	€1,912,155	€5,635,423
Buildings	50	0.02	7030	€11,999,597	€2,233,376
Other structures	50	0.02	633	€61,568	€13,759
Road	50	0.02	302,982	€8,278,693	€1,817,806
Specialized land use	50	0.02	400,000	€683,770	€319,958
Unspecialized land use	50	0.02	14,440,000	€23,938,683	€2,291,826	€6,676,725
Buildings	100	0.01	7931	€13,425,906	€2,678,795
Other structures	100	0.01	633	€61,568	€14,847
Road	100	0.01	350,619	€9,558,971	€2,125,702
Specialized land use	100	0.01	524,000	€861,115	€412,420
Unspecialized land use	100	0.01	16,680,000	€27,819,598	€2,826,192	€8,057,956

. The first one is related to modeling with only storm surge input, while the second one refers to the combined model (storm surge and rainfalls). The final values of risk, or economic losses, corresponding to a different return period of flooding events are also shown in Figure 10 in the form of risk curves.

Table 4 and Table 5 show that economic losses corresponding to the first scenario, flooding at one-year return period, are already relevant if compared with the other three scenarios. Therefore, it represents the risk scenario to be considered to plan mitigation measures. Moreover, the amount, in square meters, of structures (buildings and other structures) undergoes slight variations among the last three scenarios (T_R = 30, 50 and 100 years). This is due to the configuration of flooding areas, which, as shown in Figure 9, are nearly similar for these scenarios. This effect is even more evident for the storm surge modeling (in Table 4), observing the values of total economic losses for 30-, 50-, and 100-year scenarios. Land elevation in coastal areas subject to inundation from storm surge is below mean sea level and, hence, a low return period event is sufficient to completely flood the urban area.

Furthermore, the elements at risk mostly affected by flooding events are agricultural areas and roads, which represent preferential flow paths. Finally, the risk curves in Figure 10 highlight that the combined flood modeling provides more complete results in terms of potential economic losses. This aspect is confirmed in Figure 7, which shows that the storm surge alone could not sufficiently justify the flooding areas.

5. Conclusions

Quantitative risk estimation methods are often limited by the lack of data on calamitous events, information on damages, and related costs, although an increase availability of data is lately being recorded in the technical literature with regard to risk related to natural hazards.

In this paper, a “combined hydraulic modeling” is proposed to identify the hazard map in an ungauged coastal area. In contrast with most studies, in which the risk assessment in a coastal area only takes into account the storm surge contribution in adherence to the classical approach, this study shows that in low lying areas with presence of water channels the classical approach would lead to an underestimation of the economic risk. To improve the estimate of the economic risk, a two-dimensional hydraulic model was adopted to assess the coastal inundation risk zones in the area along the Adriatic coast. The model enabled the simulation of a storm surge approaching a low sandy beach simultaneously affected by river flooding.

Firstly, the output of the model was validated by employing real data to demonstrate its effectiveness in identifying coastal areas exposed to storm surge inundation.

Secondly, the output of flooding simulations was combined with the distribution of the elements at risk in the study area to determine the level of exposure. This procedure was applied for both flooding simulations, i.e., in the presence of only storm surge contribution and in combination with rainfall. Subsequently, the amount of exposed assets potentially involved in the flooding events and the related economic values were computed for each flooding scenarios and for each category of elements at risk.

Finally, the total economic losses associated with each flooding scenario were computed by summing the consequence values of all categories of elements at risk.

The results show that economic losses corresponding to the first scenario analyzed are already relevant when compared with the other three scenarios for floods of one-year return period. These results represent the risk scenario to be considered to plan mitigation measures. The results also show the surface area amount, in square meters, of civil structures subject to slight variations in the last three scenarios considered with return periods of 30, 50, and 100 years. This becomes even more evident for storm surge modeling when observing the values of total economic losses for the 30-, 50-, and 100-year return period scenarios. The particular trend of the risk curves determined can be explained by the peculiar topography of the flooded area are characterized by negative elevation respect to the medium sea level. This is the reason a low return period event can also be sufficient to completely flood the urban area.

The risk curves for the flood combined modeling provide more complete results in terms of potential economic losses. This was confirmed by the necessity of the combined hydraulic model which is proposed in present work.

The approach adopted in this work can deliver a concrete support to decision-makers responsible for the identification of measures aimed at reducing the flood risk and for the implementation of flooding protection policies detailing safety interventions. A further development of this work will account for the indirect and intangible losses at basin scale in accordance with the 2007/60 Flood Directive.