A Multi-Risk Methodology for the Assessment of Climate Change Impacts in Coastal Zones

Abstract

Climate change threatens coastal areas, posing significant risks to natural and human systems, including coastal erosion and inundation. This paper presents a multi-risk approach integrating multiple climate-related hazards and exposure and vulnerability factors across different spatial units and temporal scales. The multi-hazard assessment employs an influence matrix to analyze the relationships among hazards (sea-level rise, coastal erosion, and storm surge) and their disjoint probability. The multi-vulnerability considers the susceptibility of the exposed receptors (wetlands, beaches, and urban areas) to different hazards based on multiple indicators (dunes, shoreline evolution, and urbanization rate). The methodology was applied in the North Adriatic coast, producing a ranking of multi-hazard risks by means of GIS maps and statistics. The results highlight that the higher multi-hazard score (meaning presence of all investigated hazards) is near the coastline while multi-vulnerability is relatively high in the whole case study, especially for beaches, wetlands, protected areas, and river mouths. The overall multi-risk score presents a trend similar to multi-hazard and shows that beaches is the receptor most affected by multiple risks (60% of surface in the higher multi-risk classes). Risk statistics were developed for coastal municipalities and local stakeholders to support the setting of adaptation priorities and coastal zone management plans.

1. Introduction

Nowadays, international and European policies for disaster risk reduction, sustainable development, and climate adaptation are posing increasing emphasis on systemic risks [1]. The need of holistic risk management was first mentioned in the Hyogo Framework for Action [2] and then confirmed in the UN’s Sendai Framework for Disaster Risk Reduction [3]. Furthermore, multi-hazard risk management was recognized by the Post-2015 Development Agenda as an important way to achieve the Sustainable Development Goals (SDG) and the Millennium Development Goals (MDG) [4]. Similarly, the European Civil Protection Mechanism emphasized the need to address all natural and man-made hazards affecting the environment, society, and economy at multiple spatial and temporal scales [5]. More recently, the IPCC special report on the impacts of global warming of 1.5 °C above pre-industrial levels points out that some regions of the world could be increasingly affected by concomitant changes in several types of hazards [6]. Therefore, the multi-risk paradigm is gaining importance also within the climate adaptation policies (e.g., the European Adaptation Strategy [7]) and the ambitious goals set by the European Green Deal Green Deal [8].

Even if it is recognized that the assessment of multi-risk can make a great difference in building resilience and sustainable development pathways, the achievement of multi-risk governance in the current disaster management and climate adaptation practices is still a big challenge [9].

Multi-risk assessments are particularly relevant in coastal zones that are vulnerable to climate change and several interacting hazards and impacts (e.g., sea-level rise, coastal erosion, and storm surge) [10,11]. Although several methods and tools (e.g., vulnerability index and indicators, GIS-based decision support systems) have been already developed to assess and manage coastal vulnerability to climate change and natural hazards [12], they cannot be easily applied for an integrated assessment of multiple hazards and risks affecting different systems and sectors. Moreover, some of them (e.g., dynamic computer models) are data-intensive and require specialized knowledge and the use of complex programs [13,14].

Even if in the last decade the interest about the multi-risk assessment has increased (e.g., [11,15]), especially in applications and initiatives aimed at the assessment of risks derived from different natural and man-made hazardous events, there are still few examples of applications in the climate change community [16]. A detailed description of the methodologies available in the scientific literature, as well as their assumptions and limitations, can be found in [17,18]. The more practical approach, representing a screening of the investigated hazards/risks, is the qualitative multi-risk method, allowing a first estimate of risks based on end-user and expert judgment (e.g., [19,20,21,22]). The semi-quantitative methods (e.g., weighted sum and cause–effect matrices) provide a deeper evaluation of the relationships between agents and processes [22] and the exposures of elements at risk [23]. Finally, quantitative multi-risk methods (e.g., probabilistic approaches, Bayesian networks, and agent-based models) allow a more robust analysis of the multi-risk dynamics [19,20,21,22,24,25] based on more sophisticated computational models [26,27,28].

Scientific challenges related to the application of multi-risk methodologies to climate-driven impacts are related to understanding the combination of multiple hazards contributing to societal or environmental risk (i.e., the compound events) and of their probability [29]. Moreover, current issues are related to: the identification of key variables and thresholds for the assessment of the impacts of both gradual climate changes and extreme events; the evaluation of multiple elements potentially at risk and of their interconnections; and the assessment of dynamic vulnerability to multiple hazards that may affect, jointly or consecutively, the same target [16]. Finally, a crucial point is to develop multi-risk information that is usable and easily understandable by stakeholders and decision-makers, in order to define sustainable adaptation pathways [9].

Going beyond the traditional risk-by-risk approach commonly used to evaluate the effect of individual hazards in coastal areas (e.g., [30,31,32,33,34,35,36]), this paper presents a pioneering Multi-Risk Assessment (MRA) methodology to measure the effect of climate change on multiple (interacting) hazards on exposed vulnerable sectors. The methodology allows:

• identifying a set of common metrics and operative steps for the assessment of multi-hazard risks in coastal areas, under current and future climate conditions;
• classifying the vulnerability of natural and human elements at risk to different types of hazards; and
• providing a relative ranking of hotspot multi-risk areas at the regional (sub-national) scale, employing Multi-Criteria Decision Analysis (MCDA).

It has been conceived as a screening risk tool to be used in coastal areas without the need of excessive data collection or computer processing capacity, with the final goal of mainstreaming a multi-hazard risk perspective in decision-making practices. The implementation of the MRA methodology in a test case makes it possible to quantify, integrate, and compare multiple risk pathways for a selected set of hazards, showing how cumulative risk analysis can be useful to improve cross-sectorial decision-making, climate proofing, and adaptation planning.

After presenting the pilot case of the North Adriatic coast and the available input data (Section 2), the paper introduces the methodological steps of the multi-risk methodology (Section 3), the discussion of results (Section 4 and Section 5), and the concluding remarks (Section 6).

2. The North Adriatic Coastal Area

The study area considered for the application of the MRA methodology (Figure 1) comprises the Italian coast of the North Adriatic region from the national border between Italy and Slovenia to the mouth of the southern tributary of the Po Delta system (i.e., Po di Goro). It includes the coast of the Veneto and Friuli Venezia Giulia Regions, for an overall length of about 286 km along the Adriatic Sea [37]. It covers the territory below the 2 m elevation contour line, for a total surface of 3531 km². From a morphological point of view, the coast includes a variety of coastal landforms spanning from the low-lying sedimentary coasts of the Po River Plain to the rocky and higher shorelines of Trieste, interrupted by several river mouths and lagoons. Moreover, it hosts several environmental and socioeconomic systems highly vulnerable to climate change (e.g., the Po River Delta, the historical city of Venice, wetlands, and habitats of international importance). Some examples of climate-related impacts affecting the region are the intensification of extreme flood events in the city of Venice [38,39], which have been particularly damaging in autumn 2018 and 2019 (www.comune.venezia.it). Moreover, the heavy rainfall events and the related pluvial floods in urban areas [40] as well as the risk of water quality deterioration and loss of ecosystem services in coastal [41] and marine waters [42] are an increasing concern for the region. Despite the widespread presence of artificial protections (e.g., seawalls, breakwaters, submerged defenses, and river banks in the Po River Delta), the possible effects of coastal flooding and erosion remain a big issue for coastal zone management [43,44,45] and are expected to increase in the future due to climate change [46,47]. For these reasons, the Italian Ministry of the Environment recognized the North Adriatic coast as a fragile region, highly vulnerable to climate change and sea-level rise impacts [48]. Even if several measures such as the mobile barriers designed to protect Venice from high tides (i.e., the MOSE project that should be completed by the end of 2021: www.mosevenezia.eu) and other types of natural and semi-natural coastal protections (e.g., beach and sandbanks nourishment, dune restoration, and creation of new dunes with the use of natural vegetation) [49,50,51,52,53,54] are already in place to protect the study area, sustainable and integrated coastal planning should be reinforced to consider the impacts that could amplify one another across the land-sea interface [55].

Differently from recent studies focusing on single-hazard risk assessments in the North Adriatic region (e.g., [30,32,37,43,55,56,57]), this study is a first attempt to adopt a multi-hazard risk perspective to analyze the cumulative impacts of climate-related impacts in the North Adriatic coast, using publicly available datasets at the regional (sub-national scale).

2.1. North Adriatic Dataset for Multi-Risk Assessment

Several data in graphic format or database were requested and retrieved from previous research projects and from local public authorities (mainly Veneto and Friuli-Venezia Giulia region) to characterize climate related hazards. The elements potentially at risk (i.e., receptors) and their vulnerability single hazard scenarios necessary to apply the multi-hazard approach in the North Adriatic coastal zone were selected from recent risk assessment studies developed for the North Adriatic coast [34,35,36,56,58], as summarized in

Table 1. Summary and technical description of single hazard scenarios employed in multi-risk assessment for the North Adriatic coast.

Single Hazard Map	Description	Technical Features
Sea-level rise (source: [35])	Maps of inundated areas (cm) according to projected sea-level rise scenarios for the year 2100 (spatial resolution 25 m).	- Emission scenario: A1B SRES; - Climate-forcing model: EBU-POM (spatial resolution 28 km). - Impact Model: SHYFEM high resolution ocean and sea circulation model (Spatial resolution: 2.5 km–50 m); - Reference period: 1960–1990 - Time scenario: 2070–2100; - Hazard metric: projected sea-level at 2100 (low 17 cm, high 42 cm).
Storm surge (source: [58])	Maps of inundated areas (cm) according to different storm surge return periods and sea-level rise scenarios (spatial resolution 25 m).	- Method: Joint Probability Method (JPM, Pugh and Vassie 1979); - Storm surge return period: 20, 50, 100, 200, 500 years; - Spatial resolution: 28 tide gauge stations; - Hazard metrics: projected sea-level at 2100 (low 17 cm, high 42 cm, from the impact model SHYFEM), mean sea-level, astronomical, and meteorological tide (from the JPM).
Coastal erosion (source: [36])	Maps of seasonal coastal erosion hazard for the thirty-year period 2070–2100 (spatial domain: 1 km from the shoreline, spatial resolution 25 m).	- Emission scenario: A1B SRES; - Climate-forcing model: COSMO-CLM Regional Climate Model (spatial resolution 14 km); - Impact model: coupled wave-ocean ROMS-SWAN model (spatial resolution: 5–2 km); - Reference period: 1960–1990 - Time scenario: 2070–2100; - Hazard metrics: wave height, bottom stress.

.

The single hazard maps described in Table 1 were obtained following the single-hazard Regional Risk Assessment approach implemented by the DEcision support SYstem for COastal climate change impact assessment (DESYCO) [57]. These maps integrated individual hazard metrics produced by impact models (e.g., projected sea-level rise scenarios and wave height) with pathway (e.g., elevation of the territory) and attenuation factors (e.g., artificial protections) evaluating the potential areas affected by sea-level rise inundation, storm surge flooding, and coastal erosion. Specific information about the impact models used in the assessment can be found in the works of Umgiesser et al. [59], Carniel et al. [60], and Torresan et al. [35].

The sea-level rise hazard maps produced for the North Adriatic area allowed the identification of areas potentially inundated by sea-level rise at 2100, ranking the hazard intensity based on the projected water height over each grid cell [34]. The purpose of these maps is to provide a preliminary look at the extent of areas prone to coastal inundation, based on the highest accuracy elevation data publicly available at the regional scale. Digital Elevation Models (DEMs) used in the assessment were derived from regional technical offices of Veneto and Friuli-Venezia Giulia with a horizontal resolution of 5 and 10 m, respectively. Both DEMs were elaborated from technical regional maps (contour lines and ground control points). To homogenize the spatial resolution of all the raster maps used in the risk assessment, both DEMs were converted into 25-m grid layers that were already considered as a suitable resolution for meso-scale regional risk assessments [61,62]. Given the flat geology of the coastal plain, spatial interpolation was performed without significant loss of accuracy compared to the original data. However, to have more precise multi-risk results, it is recommended to refine the analysis using inundation maps based on a higher resolution DEM (e.g., Light Detection and Ranging (LiDAR) data) or quantitative hydrodynamic models.

Moreover, even if subsidence is an important factor for the assessment of future inundation risk, especially for the low-lying coastal plains of the Venice lagoon and the Po River Delta (with values of 1–2 mm/year) [63,64], this factor could not be considered in the present analysis due to the lack of homogeneous data mapping the phenomenon for the overall coast of Veneto and Friuli-Venezia Giulia Regions. Moreover, as highlighted by Rizzi et al. [58], extrapolating subsidence values for the future based on past trends could lead to big errors, especially when the main causes of subsidence are anthropogenic, such as on the North Adriatic coast.

The storm surge maps estimated the potential flooded areas due to extreme events of different return periods (e.g., 20, 50, 100, 200, and 500 years). In addition to the effect of astronomical and meteorological tides, the assessment considered the potential effect of sea-level rise in 2100, the elevation of the territory, and the distance from the coastline [56].

Finally, as described by Gallina et al. [36], the seasonal coastal erosion maps (i.e., spring, summer, autumn, and winter) represent coastal areas potentially exposed to higher rates of coastal erosion due to climate change in the future scenario 2070–2100, considering the action of waves and bottom stress, as well as the presence of artificial protections and the distance from the sea [34,35].

All single hazard maps consider the A1B emission scenario as climate forcing, representing an intermediate case compared to the more intense A2 and the weaker B1 storyline families [65]. The A1B projections are equally plausible and comparable with the more recent simulations of the Representative Concentration Pathway (RCP) 6 [66], an intermediate forcing pathway between RCP4.5 and RCP8.5. In fact, the overall performance of the Coupled Model Intercomparison Project Phase (CMIP) 5 ensemble [66] shows similar values to that of the older CMIP3 ensemble [65] considering the mean global temperature projections, which range 1.7–4.4 °C for the A1B and 1.8–3.7 °C for RCPs [67,68]. In addition, the projected sea-level rise used in this application (i.e., 42 cm) renders a medium value compared to the Special Report on Emissions Scenarios (ranging 18–52 cm, from [65]) and to the more recent RCPs (i.e., 26–98 cm from [10]).

As far as sea-level rise is concerned, there is a high uncertainty on future projections due both to the greenhouse gas emissions considered and to the assumptions done about the dynamics of terrestrial ice melting [45]. The present application employs sea-level rise projections produced by a dedicated model-chain composed of climate and high-resolution ocean/sea circulation models for the North Adriatic Sea [34,35], which allow properly evaluating the consequences associated to sea-level rise, storm surge and coastal erosion.

Various physical and environmental data are needed to characterize the spatial pattern and distribution of targets (e.g., beaches, wetlands, and agricultural areas) and define appropriate indicators of multi-vulnerability (e.g., vegetation cover, percentage of urbanization, and presence of protected areas) to cumulative climate change impacts in the studied area. The available dataset for the application of the exposure and multi-vulnerability assessment to the case study area of the North Adriatic coast is summarized in

Table 2. Available dataset for the characterization of exposure and vulnerability in the North Adriatic coast (VE, Veneto Region; FVG, Friuli Venezia Giulia Region).

Dataset	Spatial Domain	Reference
Land cover Map (1:10,000)	VE	VE Region, 2009
Monitoring Land Use/Cover Dynamics (MOLAND) (1:25,000)	FVG	FVG Region, 2000
Digital Elevation Model (DEM) (5 m)	VE	VE Region, 2007
Digital Elevation Model (DEM) (10 m)	FVG	FVG Region, 2006a
Natural reserves, Regional Parks, Sites of Community Importance (SIC)/ Special Protection Areas (ZPS) (1:150,000)	VE	VE Region, 2005
	VE	VE Region, 2008
	VE	VE Region, 2006
	FVG	FVG Region, 2007
	FVG	FVG Region, 2008
Soil type, Geologic map (1:100,000)	VE	VE Region, 2009
Soil type, Geologic map (1:150,000)	FVG	FVG Region, 2006b
Population and Housing Census	VE	ISTAT, 2001
	FVG	ISTAT, 2001

.

All case study data were pre-processed using GIS spatial operators to provide a homogeneous raster dataset for the application of the environmental multi-risk methodology with the final resolution of 25 m × 25 m, representing the highest feasible detail, according to scale of the analysis and the available dataset.

3. Multi-Risk Methodology for the Assessment of Natural and Climate-Related Impacts

The main aim of the MRA methodology is to estimate the cumulative impacts of climate-related hazards, providing a rapid assessment tool able to scan areas and targets at risk from different types of hazards, on a multiplicity of elements at risk (e.g., beaches, wetlands, population, urban and agricultural areas, etc.). The final goal of such methodology is to mainstream the multi-hazard risk perspective in decision-making practices, as claimed by the most recent international policies on disaster risk reduction and climate adaptation [3,4,5,7].

Specifically, the multi-risk methodology pursues the following objectives:

• to identify multi-hazard risk and vulnerability indicators suitable for the assessment of current and future climate change impacts on vulnerable natural and human targets;
• to provide a relative ranking of areas and targets more heavily affected by multiple risks related to climate change at the regional (sub-national) scale, employing expert-based Multi-Criteria Decision Analysis (MCDA) tools; and
• to assist local communities and stakeholders (e.g., water, soil, and coastal management authorities) in the definition of sustainable climate adaptation and disaster risk reduction strategies.

In line with the recent definitions of the IPCC [10,11], risk is considered as the potential adverse consequences for natural and human systems resulting from the interactions of climate-related hazards with vulnerabilities of the exposed systems (i.e., elements at risk). Moreover, going beyond the traditional hazard-by-hazard approach usually proposed by the coastal risk assessment community (e.g., [36,56,69,70,71]), the risk framework of this paper adopts a multi-risk perspective, considering that different hazards can be interrelated, potentially affecting a variety of receptors (i.e., elements at risk) featured by different (hazard-dependent) vulnerabilities. Specifically, according to the recent glossary proposed by the UNISDR [72], multi-hazard refers to the hazardous (chronic or acute) events that may occur simultaneously, cascadingly, or cumulatively over a certain region and timeframe, taking into account the potential interrelated effects, while multi-vulnerability considers how the sensitivity of receptors changes in relation to different hazards [16]. Figure 2 summarizes the four main steps of the MRA methodology: multi-hazards, exposure, multi-vulnerability, and multi-risk assessment.

As described in the next subsections, each step of the MRA requires the management of heterogeneous input data with different units of measures (e.g., climate projections, hydrodynamic simulations, publicly available land use, and geo-biophysical data) that are first geo-referenced and converted into homogeneous grid maps (raster data), and then normalized to ensure comparability across different spatial layers and their aggregation with MCDA. According to an acknowledged practice in natural hazard and environmental risk assessment [73], our MRA integrates the use of expert knowledge when a quantitative estimate of hazards or vulnerabilities is not achievable due to lack of resources or data. The experts involved in the MRA testing in the North Adriatic Coast are mainly the authors of the paper, representing a multi-disciplinary team of environmental scientists, geographers, and natural hazard and climate change risk experts, as well as experts with extensive experience in the application of Multi-Criteria Decision Analysis (MCDA). The experts’ opinions were first based on the survey of scientific literature and historical data on coastal erosion and flood hazard risk assessment on the case study area, and then complemented by discussions with other experts in relevant fields such as hydrodynamic and climate modelers, as well as local stakeholders and territorial experts (mainly civil engineers) from regional/local coastal zone management authorities (Table S1) who joined several workshops in the framework of national and European projects where the intermediate results of the analysis were presented and discussed (e.g., FP7 CLIM-RUN and PEGASO projects [74,75]). The output of the analysis (including multi-risk maps and statistics presented in Section 4) could support the definition of comprehensive climate adaptation and disaster risk management strategies.

3.1. Multi-Hazard Assessment

The multi-hazard assessment evaluates the potential influences among different types of hazards, taking into account possible interactions that may amplify the overall risk [76].

It is performed through the following procedure: (1) select the hazards and timeframe of analysis; (2) assess the hazard interactions; (3) analyze hazard probabilities; and 4) aggregate and normalize the multi-hazard score.

The hazards and timeframe to be considered in the MRA depend on the scope of the study and the data availability. Single hazard metrics can be derived from climate models (e.g., temperature or precipitation projections) and physical impact models (e.g., sea-level rise, currents velocity, and bottom stress) [34,35], or by reviewing single-hazard scenarios available from previous projects in the considered region. The selected metrics should represent the intensity of the hazard (e.g., projected water level for inundation and wave energy for coastal erosion) and capture the influence of climate change (e.g., anomalies between baseline and future scenarios).

In our analysis for the North Adriatic coast, the hazard scenarios for the future timeframe 2070-2100 were retrieved from previous single-hazard assessment (Table 1) and combined in the MRA to represent the worst conditions, following a precautionary approach. More specifically, the higher sea-level rise projection of 42 cm at 2100 was considered in the multi-hazard assessment, associated with the more extreme storm surge event (500-year return period) and the winter coastal erosion scenario (representing the worst pattern for the bottom stress and wave height hazard metrics). Heterogeneous hazard scenarios expressed in different units of measure (e.g., sea-level rise, wave height, and bottom stress) were first normalized in the same 0–1 scale over the study area, and then integrated in the assessment of hazard interactions.

The analysis of the relationships among the hazards (e.g., how a hazard could be negatively affected by another operating in the same area and in the same temporal window) is performed by means of an hazard influence matrix [22] (Table S2). This matrix evaluates how much the different hazards (h_i) influence each other (i.e., the synergies among the hazards) by means of weights (w_a,b) For example, h₁ influences h₂ with a magnitude represented by w_1,2 and vice versa. If the synergy among h₁ and h₂ is higher (i.e., very strong), w_1,2 should assume the maximum value of 1, while a missing interaction among the investigated hazards will be translated in the lowest value (i.e., 0). Where possible, a quantitative estimate of hazard interactions should be addressed in a quantitative way (e.g., by using physical modeling of hazard (co-)occurrence). However, in situations where resources and data are limited, scores and weights are usually assigned based on scientific literature, professional intuition, and expert judgement [72].

In this application, a team of environmental risk experts represented by the authors of this paper performed the assignation of the influence weights (Table S5). The weighting procedure was performed considering the recent scientific evidences about storm surge flooding, sea level rise, and coastal erosion risks in the North Adriatic region [35,36,58]. Qualitative and individual evaluations involved in the process can limit the validity and reproducibility of the results. It is worth noting that the common scale adopted for the assignation of weights to hazard interactions removes a certain amount of subjectivity, adding credibility to the multi-hazard analysis (Table S3).

Table S5 highlights that sea-level rise is not affected by the influence of storm surge and coastal erosion (i.e., influence weight equal to 0), as, considering the hydrodynamic model SHYFEM used for this application (Table 1), it is mainly forced by regional climate variables related to atmospheric pressure and winds [35]. Therefore, the hazard influences considered in the assessment include sea-level rise affecting storm surge and coastal erosion, and storm surge affecting coastal erosion. Specifically, an influence weight equal to 0.8 was assigned to the influence of the sea-level rise to storm surge events due to the strong effects that sea-level rise could have on tidal range and the evolution of storm surges, exacerbating the intensity of flood events [11,58,77].

Sea-level rise will also affect coastal erosion; however, in the literature, there are still controversial studies for what concerns the nature of this relationship (e.g., [10,36,78,79]). Therefore, an influence weight equal to 0.5 was assigned to the relationship between sea-level rise and coastal erosion, providing a certain but intermediate importance to the synergy that the two investigated hazards could have. Finally, looking at the relation between storm surge and coastal erosion, an influence weight equal to 0.8 was provided. In fact, the effects that storm surges could have on coastal erosion are quite certain in the scientific literature: they will increase the probability and the rate of the beach removal (e.g., [2]).

The assessment of hazard interactions acting on the same spatial units (i.e., the raster grid cells) is then performed according to Function 1 in Table S4. If there are no relationships among the hazards, the cell maintains its initial score. Otherwise, the multi-hazard weighted score (h_i^’) of each cell is increased considering the potential synergy of overlapping hazards, represented by hazard influence weights. To provide a normalized multi-hazard score (h), all values (h_i^’) are aggregated into a weighted sum (Function 2 in Table S4).

Afterwards, the expert should define the probabilities of single hazards happening in the same timeframe (not necessary simultaneously), based on the available hazard data. The best option is to derive the probabilities from numerical modeling, but, unfortunately, there is still a significant lack of information about probabilities associated to climate change projections and related hazards [16]. The second option is to consider the analysis of past events happened with the intensity of concern in a given timeframe: the more numerous are the events, the higher is the probability of the related hazard. In addition, the inverse of the return period could be used to represent the probability of a hazard event as already suggested by several studies and directives concerning natural hazards risk assessment (e.g., [73,80]). However, if no probabilities can be defined for each hazard, the most precautionary and conservative approach is to consider a probability score equal to 1, assuming the investigated hazard as certain event in the timeframe. For the testing application of the MRA in the North Adriatic coast, a value equal to 1 was assigned to the sea-level rise and coastal erosion scenarios, as it was not possible to reconstruct the probabilities from historical data or numerical models for these phenomena. Instead, for the storm surge extreme event, the inverse of the return period (i.e., 1/500 = 0.002) was defined as probability value in the multi-hazard assessment.

Once probabilities of single hazards are defined, the disjoint probability is applied to each cell obtaining a total multi-hazard probability value p. Probabilities of each individual hazard is generally expressed as a number between 0 and 1, where 0 indicates that the hazard is not affecting the investigated cell (impossible outcome), and 1 indicates that the hazard will certainly happen.

Function 3 (Table S4) considers that the hazards could happen both individually (i.e., a single hazard h₁ or h₂ will affect the cell) or simultaneously (e.g., h₁ happens together with h₂ in the same cell). Accordingly, the higher is the number of the hazards affecting the cell, the higher is the total probability.

The last step of the multi-hazard assessment is aimed at the aggregation of the hazard scores and probabilities. Therefore, the final score is comprehensive of the normalized multi-hazard weighted score (h) and of the probability of occurrence of the hazards affecting the cell, according to Equation (1):

h^p = h * p

(1)

where h^p is the multi-hazard score with the associated probability; h is the multi-hazard score of the investigated cell weighed, normalized in [0,1] (Function 2); and p is the probability of the investigated cell according to the disjunctive probability function (Function 3).

The assignation of the influence weight allows the implementation of Functions 1 and 2 (Table S4), which provide a weighted and normalized multi-hazard score for each cell of the study area. The score, visualized in the multi-hazard map (Section 4.1), ranges from 0 to 1, where 0 means no hazards affecting the cell and 1 represents the cells where the number, weights, intensities, and probabilities of hazards are relatively higher.

3.2. Exposure Assessment

The main aim of the exposure assessment is to identify all the elements at risk (i.e., the receptors), representing valued natural or human attributes that can be threatened by the hazards [76]. The process considers all the relevant features of the territory such as environmental resources, physical infrastructure, socioeconomic activities, and human health.

The exposure score (e) is evaluated as follows:

$$e=\{ \begin{array}{c} 0 if no receptors is present in the investigated cell\\ 1 if there is the presence of one or more receptors\end{array}$$

(2)

According to the main characteristics and vulnerabilities of the territory, the scale of analysis and the available dataset (Section 2), the exposure assessment in the North Adriatic coast was focused on environmental receptors including beaches, river mouths, wetlands, protected areas, natural/semi-natural systems, and agricultural area. Relevant socioeconomic elements such as infrastructure and population were also considered in the analysis, selected as potential risk hotspots in urban areas.

The result of this step is a map showing the localization and the geographic extent of all the investigated receptors, including areas where multiple receptors are located and therefore more socioecological values are potentially at risk.

3.3. Multi-Vulnerability Assessment

The main aim of this step is the assessment of receptors’ vulnerability to multiple hazards in the case study area. It considers the presence of different elements at risk and their physical-environmental predisposition (sensitivity) to be adversely affected by multiple natural and climate-related hazards, which means that different receptors are differently sensible to the same risk [36].

The assessment is based on the selection and aggregation of multiple physical and environmental vulnerability factors (vf) representing the sensitivity of receptors to the different investigated hazards. According to Torresan et al. [37], this step is carried out through the development of a vulnerability matrix in which the main column is an inventory of the selected natural and climate-related hazards affecting the analyzed area and the main row lists the major receptors at risk. For each receptor, a set of vulnerability factors (vf, grey cells in the intersection between a hazard and a receptor) is defined representing physical, ecological, and socioeconomic parameters (e.g., slope, vegetation cover, wetland extent, and percentage of urbanization) relevant for measuring vulnerability to different hazards. The same relative scoring system and guideline is used to normalize quantitative and qualitative vulnerability factors in a common scale where 0 indicates no vulnerability and 1 indicates the maximum vulnerability in the analyzed region and case study (Table S3). Then, the vulnerability assessment continues with the aggregation of normalized vf for each cell of the case study. If a cell is characterized by one receptor (e.g., r₁) that is affected by one hazard (e.g., H₁), only the biophysical and environmental vulnerability factors related to H₁ and r₁ of the matrix are considered. If a cell is characterized by two or more receptors and two or more hazards, the aggregation process considers all the vf identified for those receptors and hazards taken only once.

Moreover, individual vulnerability factors can be weighted to represent their relative importance in the estimation of vulnerability associated with selected natural and climate-related hazards. In addition, in this case, experts are asked to provide a hierarchy of importance to vulnerability factors according to a common 0–1 scale (Table S3). This means that a specific vulnerability factor (vf₁) can be highly relevant for one hazard (e.g., weight 1) and strongly less important for another hazard (e.g., 0.4).

Therefore, the final weight for a vulnerability factor vf, possibly part of several receptors and affected by several hazards, is defined as the maximum of the weights of obtained by each vf in relation to each hazard H_i.

The aggregation of the weighted vulnerability factors is then performed by applying the “probabilistic or” operator [81], defined as follows:

$$\mathit{v}\mathit{f}=\underset{\left|\mathit{V}{\mathit{W}}_{\mathit{c}}\right|}{⨂}\left(\mathit{V}{\mathit{W}}_{\mathit{c}}\right)$$

(3)

where:

• vf = physical and environmental multi-vulnerability score of the investigated cell to the investigated impact(s).
• $\mathit{V}{\mathit{W}}_{\mathit{c}}$ = set of physical and vulnerability factors weighted by the corresponding ${\mathit{w}}_{\mathit{f},\mathit{c}}$. They are related to all impact(s) and receptors present in the investigated cell c, i.e. $\mathit{V}{\mathit{W}}_{\mathit{c}}=\left\{\mathit{v}{\mathit{f}}_{\mathit{c}}\cdot {\mathit{w}}_{\mathit{v}\mathit{f},\mathit{c}}|\mathit{v}{\mathit{f}}_{\mathit{c}} is a vulnerability factor present in cell c\right\}$.

Shown below is an example of the “or” operator at work (Equation (4)):

$${⨂}_{\mathit{i}=\mathbf{1}}^4\left[{\mathit{f}}_{\mathit{i}}\right]={\mathit{f}}_{\mathbf{1}} ⨂{\mathit{f}}_{\mathbf{2}} ⨂{\mathit{f}}_{\mathbf{3}} ⨂{\mathit{f}}_{\mathbf{4}}$$

(4)

where

• f_i = ith generic factor f

$$f_1 ⨂f_2=f_1+f_2-f_1 f_2=F_1$$

$$F_1 ⨂f_3=F_1+f_3-F_1 f_3=F_2$$

$$F_2 ⨂f_4=F_2+f_4-F_2 f_4={⨂}_{i=1}^4\left[f_i\right]$$

Due to the associative and commutative proprieties, the process can be repeated for all the vf that are present in each cell of the case study.

The fuzzy “or” operator is used to apply a conservative approach in the vulnerability assessment process. If just a vulnerability factor (vf) assumes the maximum value 1 (e.g., presence of dunes, shrubbery, and eroding coast), the higher score 1 will be assigned to the overall multi-vulnerability index. On the other side, the more numerous are the vulnerability factors that are present in a cell, even with low vulnerability scores, the greater is the final multi-vulnerability. For instance, if a cell has 3 vf with a score of 0.2 (e.g., advancing large wetlands in delta areas), they result in a final multi-vulnerability of 0.6.

Table 3. Multi-vulnerability matrix applied to the case study of the North Adriatic coast.

Multi-Vulnerability Matrix
	Beaches	River Mouths	Wetlands	Protected Areas	Natural and Semi-Natural Systems	Agricultural Areas	Urban Areas
Storm surge	Slope angle	Slope angle	Slope angle	Slope angle	Slope angle	Slope angle	Slope angle
	Vegetation cover	Vegetation cover	Wetland typology	Vegetation cover	Vegetation cover	Agricultural use
	Coastal typology		Wetland extent
			Vegetation cover
Coastal Erosion	Slope angle	Slope angle	Shoreline evolution	Slope angle	Slope angle	Slope angle	Slope angle
	Coastal typology	Shoreline evolution	Wetland extent	Shoreline evolution	Vegetation cover	Agricultural use	% of urbanization
	Shoreline evolution	Mouth typology	Vegetation cover	Vegetation cover
	Dunes	Vegetation cover
	Vegetation cover

shows the multi-vulnerability matrix developed for the North Adriatic coast.

Vulnerability factors were identified considering only storm surge and coastal erosion hazards, since a sea-level rise inundation event was assumed to affect all the receptors in the same way, causing a permanent loss of all submerged receptors, based only on the elevation of the cells that is considered in the hazard map. Therefore, each cell of the territory was considered to have the same vulnerability (i.e., equal to 1) to sea-level rise.

To apply Equation (3), the classification and scoring of vulnerability factors was provided by the authors considering previous studies [36,37] and the scientific literature included in Table S6. The classification and assumptions employed in the rating of vulnerability factors is summarized in the Supplementary Materials (Table S6).

Moreover, to define the relative importance of each vulnerability factor in the multi-vulnerability of coastal areas, weights were assigned by the MCDA expert team, according to the linguistic evaluations provide in Table S3. This process tends to remove subjectivity and add credibility to the assessment, improving the communication of results among the research community and towards the end-users of information (e.g., local authorities and decision-makers). The weighting process was implemented considering also the output of a roundtable organized with several local stakeholders and experts involved in the development of the DEcision support SYstem for COastal climate change impact assessment (DESYCO [74]). The local experts joining the table (including civil engineers and coastal geologists) were aware of some preliminary results of the MRA, which have influenced their opinion during the assigning of weights.

Overall, the local stakeholders and experts have joined a series of workshops organized by the authors of the paper within the FP7 research projects CLIM-RUN and PEGASO. As specified in the related publications [74,82], the local stakeholders and experts engaged in the MRA included public institutions with a specific mandate for Integrated Coastal Zone Management (ICZM) in the Veneto and Friuli-Venezia Giulia Regions (e.g., environmental protection agencies, managers of protected areas, and representatives of geological surveys and regional meteorological services). The participatory process engaged 37 people involved by means of surveys, local workshops, and focus groups (Table S1). The early and continuous stakeholders’ involvement accompanied the development and refinement of the MRA framework, and turned out to be very useful to frame the right research questions in terms of future time scenarios, geographical scale and resolution, choice of hazards, receptors, vulnerability factors, and risk thresholds and weights to be considered in the analysis. Moreover, it was useful to identify relevant research products (risk maps and statistics) responding to real end-users’ needs.

The assigned weights are reported in Table S7. Coastal typology, slope angle, and shoreline evolution gained the higher weights (i.e., 0.8) in the present assessment as they represent geo-physical characteristics very important for the assessment of receptors vulnerability, compared to other vulnerability factors; vegetation cover and dunes were considered to have medium-high weight (i.e., 0.6); and medium and low weights (i.e., 0.5 and 0.4) were assigned to the vulnerability factors most related to ecological characteristics (i.e., wetland extent, river mouth, agricultural use, and wetland typology) and the percentage of urbanization, respectively.

The result of the multi-vulnerability assessment is a map (Section 4.3) highlighting the areas more vulnerable to multiple natural and climate-related hazards according to the intrinsic characteristics of the investigated territory.

3.4. Multi-Risk Assessment

The main aim of the multi-risk assessment is a relative evaluation of the cumulative risks affecting the study area, taking into account correlation and probabilities of the analyzed hazards, the presence and localization of multiple elements at risk, and their vulnerability to the different types of hazards [76].

According to IPCC [10], the aggregation of hazard, exposure, and vulnerability scores allows the evaluation of the multi-risk in the case study, as shown in Equation (5):

$$\mathit{r}={\mathit{h}}^{\mathit{p}}\cdot \mathit{e}\cdot \mathit{v}\mathit{f}$$

(5)

where:

• r = multi-risk score;
• h^p = multi-hazard score with the associated probability (Equation (1));
• e = exposure score related to the union of the geographical area of the investigated receptors (Equation (2)); and
• vf = physical and environmental multi-vulnerability score of the investigated cell to the investigated impact(s) (Equation (3)).

According to the proposed multi-risk equation, one or more interacting hazard(s) could affect one or more receptor(s) and the vulnerability of the receptors can change according to different hazards. This approach is more useful than a single-hazard risk assessment since it considers, in a systemic way, the multiple hazards affecting a specific region, allowing the definition of more comprehensive adaptation strategies [11].

The multi-risk equation aggregates the results of each step (multi-hazard, exposure, and multi-vulnerability) and provides an overall view of the causes and consequences of the different types of hazards in a single multi-risk map (Section 4.4).

4. Results

This section presents the results of the implementation of the MRA to the pilot case of the North Adriatic coast, describing the obtained multi-risk maps and statistics as well as the potential policy implications for climate change adaptation and disaster risk reduction. As clarified in Section 2.1, the multi-risk maps derive from the combination of existing single-hazard scenarios available for the North Adriatic region. Specific details concerning the single-hazard maps employed in the assessment for coastal erosion, storm surge flooding, and sea-level rise inundation hazards are available and can be found in [35,36,58].

4.1. Multi-Hazard Maps

The final result of the multi-hazard assessment (Figure 3) highlights that circa 75% of the study area could be affected by more than one climate-related hazard. The multi-hazard map was classified using the equal interval method allowing the division of the scores into 5 equal sized classes (i.e., very low, low, medium, high, and very high). The equal interval is best applied to data ranges in order to emphasize the amount of an attribute value relative to other values. Moreover, this method is widely used for simplifying maps understanding and comparability amongst different scenarios [83].

From the maps in Figure 3I, it is possible to see that the higher multi-hazard classes are located near the coastline where there is an overlapping of the three considered hazards (i.e., sea-level rise, storm surge, and coastal erosion). Generally, the multi-hazard decreases with the distance from the coastline, going progressively from the highest to the lowest scores. In Figure 3I(a,b), the presence of artificial protections attenuates the multi-hazard score (streaked effect) near the shoreline. Figure 3I refers to the northern part of the Venice lagoon (Figure 3I(b) and the fragile area of the Po Delta (Figure 3I(c). In these areas, the high and very high multi-hazard classes denote the joint presence of coastal erosion (close to the Adriatic Sea), sea-level rise, and storm surge drivers.

Overall, the multi-hazard assessment allows the identification of coastal areas that could be affected by multiple hazards. As summarized in Figure 3II, the study area will be affected by different hazards for almost 75% of the territory (i.e., 2635 km²). Specifically, the 72% will be affected by one or two hazards (i.e., lower multi-hazard scores) and the 3% (i.e., higher multi-hazard scores) by all the three hazards.

Therefore, the multi-hazard map is useful for: getting a general picture of the investigated coastal hazards; supporting broader planning decisions in the long term; and identifying multi-hazard hotspots. Moreover, by exploring the single-hazard layers, it could be used to support sectorial planning decisions at the regional/sub-regional level (e.g., flood and coastal erosion mapping). Specifically, by integrating the effect of future climate change scenarios, the multi-hazard map is relevant for forward-looking analysis and for long-term planning decisions related to infrastructure development, human settlements and ecosystem conservation. Finally, this information could increase the awareness of coastal managers about multiple hazards affecting the same region, supporting the achievement of higher level of protections and resilience against disasters and climate change, as explicitly requested by the European Adaptation Strategy [7] and the Post 2015 Development Agendas [3,4].

4.2. Exposure Map

The exposure map (Figure 4) allows the spatial visualization of the investigated elements at risk in the North Adriatic study area.

Most of the territory is characterized by agricultural areas (ca. 70% in light green), which, in some areas, overlaps with protected areas and river mouths (e.g., the Natural reserve of the Po River Delta) highlighting complex and fragile ecosystems that should be preserved. Moreover, 11% of the North Adriatic coastal area is represented by urban areas (e.g., Venice municipality and around Grado and Marano lagoon), which includes infrastructures, buildings, and inhabitants and could be affected by multiple-risks causing damages to people, economic activities, and the environment. Other small receptors relevant at the local level are analyzed (i.e., beaches, natural and semi-natural systems, and wetlands) as they represent important features for the North Adriatic coastal area, increasing the socioecological and economic importance of this region. Overall, the assessment of exposure is instrumental to identify the elements potentially at risk (Section 4.4), such as the number of people, the total surface of buildings, and croplands located in potentially flooded areas, as required by the Sendai Framework for Disaster risk reduction (Targets A–C) [3].

4.3. Multi-Vulnerability Map

The multi-vulnerability map (Figure 5) was classified using the equal interval method allowing the division of the scores into five equal sized classes (i.e., very low, low, medium, high, and very high) [78].

Figure 5I highlights that generally multi-vulnerability is high or very over the whole territory. This is due to the widespread presence of low-sloped areas in the study area, and the predominant presence of arable lands, which are considered the more vulnerable class for the vulnerability factor agricultural use (Table S6). As shown in Figure 5II, the ranking of vulnerable receptors (including the percentage of vulnerable areas for the elements at risk identified in Section 4.2) indicates that beaches, wetlands, protected areas and river mouths have a very high multi-vulnerability score ranging from 95% for beaches (ca. 13 km²) to 52% for river mouths (ca. 422 km²), while agricultural areas are featured by 54% of surface within the highest multi-vulnerability class. In addition, the vulnerability factor “percentage of urbanization” (Table S6) plays an important role in increasing the multi-vulnerability score, especially where the rating is higher than 10%, e.g. in the hinterland area of the Venice lagoon.

It is worth noting that the low and very low multi-vulnerability classes are scarce in the North Adriatic area, with a percentage ranging from less than 1% (for protected and urban areas) to ca. 20% (e.g., for the urban areas near Trieste that are featured by steeper slopes and a lower percentage of urbanization).

Finally, it is possible to see that overall the multi-vulnerability hotspots are located in the Po Delta and in the Venice lagoon where there is an overlay of exposed sensitive receptors, such as protected areas, river mouths, and wetlands. Protected areas are highly vulnerable due to their ecological and biodiversity relevance for different European directives (e.g., Nature 2000 network and Habitat and Bird Directives).

We recognize that land subsidence can accelerate marine ingression in low-lying coastal plains. For this reason, data on subsidence rates should have been considered in the multi-vulnerability assessment. The Italian coast of the North Adriatic Sea is characterized by diffuse and uneven vertical land movements, with some local sinking hotspots [46]. In the Veneto Region, the city of Venice and its lagoon are particularly vulnerable to land subsidence [84], with large local variability [85]. Even though the historical city of Venice and the central lagoon are relatively stable (with an average land subsidence of 1–2 mm/year), the northeastern coastland bordering the lagoon is characterized by higher sinking rates (4–6 mm/year). Moreover, some hotspots (up to 30 mm/year) can be observed at the lagoon inlets, where the new structures of the MOSE works (i.e., the mobile barriers designed for the protection of the city from high tides) are under construction. The salt marshes (which naturally protect the lagoon of Venice from storm surges and erosion), are also featured by variable subsidence rates reaching 20 mm/year. Other significant—but localized—sinking movements have been observed in correspondence of newly built-up areas (e.g., the tourist villages near the city of Jesolo to the north of the Venice lagoon), highlighting the need to monitor the anthropogenic land subsidence due to the establishment of new settlements, structures, and infrastructures. South of the Venice lagoon, the city of Chioggia and its surroundings are more stable (subsidence rate of about 1 mm/year), whereas higher sinking rates can be observed near the mouths of the Brenta and Adige Rivers. Subsidence up to 10–15 mm/year have been measured in the Po River Delta [86], which is a well-known vulnerable hotspot to relative sea-level rise [87].

Going to the Friuli Venezia Giulia Region, a recent analysis of land subsidence around the coastal area of the Grado and Marano lagoon [88] showed a regional southwest gradient (going from the hinterland toward the sea) with subsidence rates reaching 4 mm/year in the area adjacent to the mouth of the Tagliamento River. Even if the average subsidence of the Grado-Marano lagoon is about 2 mm/year, the northern coastal area bordering the lagoon ranges from a stable condition to a subsidence of about 8 mm/year. Finally, the littoral zone of the Grado-Marano lagoon shows uneven subsidence rates due to the combination of natural and local anthropogenic drivers (e.g., groundwater exploitation, loads induced by newly built-up areas for the tourism sector).

The present multi-vulnerability assessment did not integrate subsidence information due to the lack of accessible quantitative data with homogeneous spatial coverage, at the time of our analysis. However, this lack would not have changed the overall multi-vulnerability score of the main hotspots shown in Figure 5I (i.e., the red areas around the Po River Delta and those bordering the central and northern part of the Venice lagoon). On the contrary, the consideration of the recent subsidence studies in the multi-vulnerability map could have increased the vulnerability score of some local sites in the norther littoral zone of Venice (Jesolo municipality) and along the shoreline of the Grado and Marano lagoon, where recent anthropic drivers of subsidence (e.g., newly urbanized areas) are contributing to local subsidence phenomena. For this reason, we recommended monitoring this parameter in future site-specific analysis, in order to refine the investigation of subsiding areas that could increase present and future vulnerability to marine ingression.

Combined with the exposure assessment, the multi-vulnerability map can be useful for local and regional authorities to identify vulnerable hotspots where the resilience to multiple hazards (i.e., permanent/temporary inundation and coastal erosion) should be increased. For instance, it can help decision-makers in a first screening of infrastructure or buildings susceptible to coastal flooding (e.g., roads and houses located in gentle sloping and low-lying areas) and in the identification of urban areas where the drainage systems could be improved. Relevant examples of measures aimed at increasing the resilience (or reducing the vulnerability) of coastal communities are represented by nature-based solutions providing multiple benefits for human and environmental well-being [89]. These measures can be set-up for different receptors to mitigate multiple hazards and include all the actions aimed at conserving, restoring and reconstructing ecosystems (e.g., reinforcement of arboreal communities in areas susceptible to flooding; restoration of wetlands to mitigate the effects of floods; and stabilization of sandy dunes with vegetation to reduce coastal erosion).

4.4. Multi-Risk Map

The multi-risk map (Figure 6I) is the final result of the MRA. Similar to the previous maps, it is classified using the equal interval method, allowing the division of the scores into five equal-sized classes (i.e., very low, low, medium, high, and very high) [71].

Generally, multi-risk scores are lower than the multi-hazard ones as they are multiplied by multi-vulnerability (i.e., scores ranging 0–1). However, in Figure 6I, the same trend of the multi-hazard map is still visible. The hotspot risk zones shown in Figure 6I(a) are beaches, which is the receptor most affected by multiple risks, with a percentage of about 60% in the very high and high multi-risk classes (Figure 6II). The medium multi-risk areas around the Venice and Marano/Grado lagoons, instead, are mainly associated to the receptors wetlands and protected areas (Figure 6I(a,b)). Finally, around the Po River Delta (Figure 6I(c)), the higher multi-risk classes are located near the shoreline while decreasing risk scores can be observed going inland, following the trend of the multi-hazard map.

Figure 6II shows the ranking of receptors affected by multiple risks and, together with the statistics included in Figure S1, supports the identification of the more affected municipalities in the study area (e.g., Staranzano, Porto Tolle, Grado, Cavallino-Treporti, Caorle, and Venice).

Overall, the multi-risk assessment can be useful to highlight hotspots for future foreshore protection (e.g., priority areas for beach erosion and sea-level defenses) within Integrated Coastal Zone Management plans. Moreover, such maps can be useful to support ecosystems management in compliance with international environmental restrictions (e.g., Nature 2000 sites) and national and regional laws for the defense of the littoral zone [88]. The multi-risk maps are particularly helpful since they provide an integrated vision of the elements at risk from multiple hazards, considering the potential effects of climate change. This systemic view can be very helpful to reduce the likelihood that risk reduction efforts targeting one type of impact increase the exposure or vulnerability to other impacts [11,88].

Additional information useful for coastal managers can be obtained from the statistics presented in Figure 7, concerning the infrastructure and population at risk in urban areas, for the coastal municipalities of the study area.

Even if all the municipalities show relatively higher percentages of infrastructure (km) and population (number of citizens) in the lower multi-risk classes, it is possible to note that highly urbanized municipalities (e.g., Jesolo and Caorle) are partly affected by higher multi-risk scores. This information can help regional authorities in the identification of the municipalities more threatened by multi-risk, supporting an efficient allocation of risk management funds.

Overall, multi-risk maps and statistics can be used as a quick scan tool to rank coastal areas and municipalities requiring higher attention for the implementation of disaster risk reduction and climate adaptation strategies (e.g., coastal sediment management plans, flood proofing, or relocation of existing developments). The results can be helpful particularly helpful for the implementation of National Adaptation Plans, as required by the European Adaptation Strategy [73].

5. Discussion

The proposed MRA methodology is a first attempt to consider how multi-hazard interactions could trigger cumulative risks in flood-prone coastal areas. This is an important advancement compared to conventional single-hazard risk assessments, since it provides a systemic and forward-looking analysis able to consider: (1) how compound hazards could trigger risk amplifications on the same spatial units and timeframe; and (2) how the vulnerability of multiple elements at risk could vary according to different types of hazards.

Our MRA is conceived as a decision support tool aiding multi-risk screening at the regional/sub-national scale. One of the strengths of the approach is that it can be applied with publicly available datasets, to offer an improved basis for understanding and communicating multi-risk patterns to local authorities, stakeholders, and decision-makers. Moreover, by varying the parameters of the MRA model, it is possible to understand the relative importance of the main factors influencing multi-risk, as well as the effect of different climate and risk management scenarios. The final goal of the assessment is to provide an overall picture of the multiple threats affecting low-lying coastal areas in view of climate change, supporting the definition of holistic adaptation strategies avoiding maladaptation (e.g., due to fragmented single-hazard or single-sector measures). Overall, the MRA employs a wide range of physical, environmental and socioeconomic indicators widely used to help decision makers in the definition of disaster risk reduction and climate adaptation policies [72,90,91,92].

However, differently from usual applications of Multi-Criteria Decision Analysis focusing on individual elements at risk, ecological functions, or services [93,94,95,96,97,98], our methodology sets up a conceptual framework considering multiple cause–effect relationships and allows progressively systemizing different sources of empirical and modeling data to achieve an holistic estimate of coastal risks affecting both natural and human systems. Moreover, compared to recent studies assessing coastal vulnerability in the North Adriatic region [35,36,58] and Mediterranean region [99], the operative steps of the MRA allow an improved understanding of multi-hazard interactions in the coastal system, and of their cumulative impact in the final categorization of risks.

Even if the proposed MRA goes beyond the traditional hazard assessment by considering the potential relations among selected types of climate-related hazards, the methodology could be improved by integrating information related to other types of natural or man-made hazards (e.g., river floods, subsidence, and contamination of water resources). Particularly, land subsidence is an important factor that should be considered in the multi-hazard analysis, since it may amplify both sea level rise and storm surges flooding and the related risk scenarios in the North Adriatic coastal plains [30,58]. However, even if specific studies have been conducted to assess the rate of subsidence over the Venice lagoon [64,85,87,88], we could not consider this phenomenon in our MRA analysis since comprehensive and homogeneous spatial data for the whole study area (i.e., the interregional/sub-national scale) were not available when we developed the study. Nevertheless, MRA is an iterative process that should be regularly updated when new datasets related to coastal hazards or vulnerability indicators are available. Another important constraint was that there is still large uncertainty in the extrapolation of future subsidence values based on past trend observations, especially when the causes of subsidence are both natural and man-made, such as for the North Adriatic region [100]. We recognize that a more detailed analysis focused on specific regions (e.g., the city of Venice or the Po River Delta) should consider subsidence as a relevant triggering factor in the multi-hazard assessment.

The assessment and management of multiple natural hazards and climate-related risks requires an extension of commonly applied risk assessment and management frameworks toward the identification and quantification of multiple stressors interactions [16]. The MCDA based approach adopted in this paper is suitable for a preliminary consideration of risk interactions (which are still not well understood) and provides a relative (not absolute) quantification of multi-risk hotspots. However, most variables described in the MCDA rely on expert knowledge and depend on the individual perception of the system, failing to provide a strong quantitative assessment and validation of results [28].

More detailed knowledge on multi-risk patterns, therefore, may be obtained by progressing to a quantitative level of analysis at a site-specific spatial scale, applying more sophisticated complex system analysis tools such as Bayesian networks, agent-based models, or system dynamics [27]. Such methods are increasingly being applied for cumulative impact assessments of climate and anthropogenic hazards in different context (mountain regions, freshwater resources management, and marine areas) since they can provide a stochastic assessment of multi-risk values, combined with an improved characterization of uncertainties [101,102].

6. Conclusions

Decision makers and local authorities need to have a clear vision of systemic risks affecting human health and safety, the environment, the economy, and society at large to develop robust adaptation pathways [103,104]. However, such assessment is notoriously difficult due to the complex and uncertain combinations of hazardous events that might happen, the multiplicity of physical processes involved, and the continuous changes and interactions of vulnerable socioecological systems. Moreover, multi-risk dynamics are characterized by complex domino effects that could amplify the overall impact on exposed elements. Those include chains of hazardous events and increasing vulnerability of receptors, together with other types of correlations among the different risk components [9]. Understanding how different hazards and vulnerabilities sources interact to determine emerging risk states is one of the open issues for the scientific community, which is attempting to shift from the single-risk paradigm toward the more challenging systemic risk concept [29].

The proposed methodology should be considered as a quick scan tool, allowing the semi-quantitative evaluation of the relationships among different hazards (and vulnerabilities) in the same area. A more detailed and quantitative assessment should consider the joint probabilities of single and compound hazards and how vulnerability changes according to different hazard thresholds [83]. Since it does not require a large amount of data, the proposed MRA is flexible and scalable in a wide range of conditions and geographical regions. Moreover, it has been conceived with a modular approach, allowing an easy upgrade to additional natural hazards, climate change impacts, or vulnerable sectors. For instance, the multi-hazard assessment that has been focused thus far on climate-driven hazards affecting the shoreline at the regional scale (e.g., storm surges, sea-level rise, and erosion) could be in principle expanded to the analysis of additional hazards relevant at the local scale (e.g., local subsidence, heavy rainfall events, and urban heat islands), given the availability of more detailed data on this phenomena. Specifically, a comprehensive analysis of vertical land movements should be considered to achieve a better understanding of the potential correlations between climate and non-climate drivers contributing to future marine ingression. Additionally, the multi-vulnerability assessment, now focused on the physical and environmental dimensions of vulnerability, can be improved by integrating adaptive and coping capacity indicators (e.g., income level, education, and safety network) and other social factors related to individual and collective behaviors and risk perceptions. Moreover, parameters and algorithms for the economic assessment could be in principle included in the multi-risk modeling, allowing the estimate of consequences in terms of direct and indirect costs.

Finally, we recognize that the exposure and vulnerability dynamics are very important for analyzing multi-risk patterns, especially in coastal areas [105]. Therefore, future improvements of our MRA methodology should provide a more realistic comprehension of the distributions of the elements at risk across space and time (e.g., land use, population, and infrastructure dynamics) and of their vulnerability.

To properly support decision-making processes and assess the efficacy of different adaptation pathways, future testing of the MRA should consider how different management options can contribute to the reduction of the overall multi-risk score.

Finally, it is important to underline that the multi-risk methodology allows a relative ranking of the cumulative investigated impacts and risks and should be re-applied each time improved input information and data (e.g., new models and more accurate site-specific indicators) become available. To this respect, a future ambitious improvement of the MRA includes the exploitation of big data analytics and Information and Communication Technologies (ICT) to run a huge amount of risk scenarios and cascading uncertainties in space and time, in order to provide robust science-based outputs for decision-making in a multi-risk perspective.