Assessing Po River Deltaic Vulnerability Using Earth Observation and a Bayesian Belief Network Model

: Deltaic systems are broadly recognized as vulnerable hot spots at the interface between land and sea and are highly exposed to harmful natural and manmade threats. The vulnerability to these threats and the interactions of the biological, physical, and anthropogenic processes in low-lying coastal plains, such as river deltas, requires a better understanding in terms of vulnerable systems and to support sustainable management and spatial planning actions in the context of climate change. This study analyses the potential of Bayesian belief network (BBN) models to represent conditional dependencies in vulnerability assessment for future sea level rise (SLR) scenarios considering ecological, morphological and social factors using Earth observation (EO) time series dataset. The BBN model, applied in the Po Delta region in the northern Adriatic coast of Italy, deﬁnes relationships between twelve selected variables classiﬁed as driver factors (DF), land cover factors (LCF), and land use factors (LUF) chosen as critical for the deﬁnition of vulnerability hot spots, future coastal adaptation, and spatial planning actions to be taken. The key results identify the spatial distribution of the vulnerability along the costal delta and highlight where the probability of vulnerable areas is expected to increase in terms of SLR pressure, which occurs especially in the central and southern delta portion.


Introduction
Coastal areas, and especially deltas, are highly dynamic systems of great importance due to high productivity, giving access to a significant amount of food and raw materials, and being rich biodiversity, hosting several ecosystems and related services [1]. In terms of spatial and temporal scales, deltas are an expression of the complex interplay between natural and human forcing, being a reflection of past and current global change phenomena along coasts as well as the societal dimensions associated with these changes [2]. Nonetheless, deltas are fragile geomorphic features that change dramatically due to human development and climate change putting them perilously out of dynamic equilibrium. These densely populated low-lying coastal plains areas are facing increasing threats and degradation due to a wide range of human [3] and natural pressures, transforming them into hazardous vulnerability index (CVI) is one of the most commonly used and simple methods to asses coastal erosion and/or inundation vulnerability [29][30][31][32][33][34]. Relevant examples of indicator-based approach at the European level include the EUROSION project that identified 13 indicators to support the assessment of coastal erosion risk throughout Europe [35]. An example of a GIS-based decision support system (DSS) is DESYCO (i.e., GIS-based decision support system for coastal climate change impact assessment). It is a tool for the assessment and management of multiple climate change impacts on coastal areas and related ecosystems based on multi-criteria decision analysis (MCDA) to identify and prioritize areas and targets at risk in the investigated region. Available methods based on dynamic computer modelling address specific coastal systems, although they attempt to deal with various coastal processes (e.g., risk assessment of coastal erosion-RACE) or the integrated assessment of coastal vulnerability (e.g., dynamic interactive vulnerability assessment-DIVA; Delft3D) [28,36].
In more recent years, efforts to assess the deltaic system vulnerability attempted to develop a more comprehensive approach, trying to move away from methods that are strictly hazard/physical phenomena based to methods approaching vulnerability in a spatially based manner and integrating relevant biophysical data with socio-economic aspects from different sources, integrating different systems replicable at various scales and in several context [30]. This includes ready-to-use coastal vulnerability approaches, such as coastal vulnerability model promoting ecosystem based management approach (i.e., integrated valuation of ecosystem services and tradeoffs-InVEST), or methods that model complex behaviors as a collection of simple "if-then" rules based on expert knowledge able to easily manage a great amount of variables from different sources (i.e., the fuzzy logic approach) [37]. However, those kinds of approaches can be integrated in defining mutual interactions and correlations among considered variables with respect to varying environmental, land use, and management drivers. The complex interrelatedness between climate, disturbance pressures, and management of the territory increases the need for reliable and consistent information for decision support to prioritize management interventions. However, the selection of an assessment method to be applied in a particular context is also strongly linked to relevant data availability.
Bayesian belief networks (BBNs) are increasingly being used to model complex coastal processes due to their ability to integrate non-linear systems, to deal with data scarcity, to combine quantitative and qualitative data, to correlate multiple variables [38][39][40], and to infer the relationships among them tackling model uncertainties with low computational cost. Accordingly, BBNs are considered as probabilistic network-graphical models that use probability and graph theory in their implementation. The advantage of probabilistic networks is that they provide explicit representations of dependencies or independencies between variables without scientific numeric or functional details, and they are a powerful tool to picture complex socio-ecological systems. Therefore, BBNs can assimilate important factors contributing to coastal change in response to diverse driving forces and can make quantitative, probabilistic predictions under uncertainty [41] that can hold up coastal monitoring and management decisions [42]. For this purpose, BBNs provide a valid alternative approach to the proven empirical and physics-based modelling [43], thus offering important advantages as decision aides to estimate the probability of vulnerability. This is proven by their use in several costal settings and applications to predict episodic coastal cliff erosion [44], to reproduce wave-height evolution in the surf zone [45], to forecast barrier islands response to storms [46,47] and dune retreat caused by coastal storms [48,49], to model hurricane damage [50] and shoreline changes [51], to assess coastal vulnerability to sea level rise and coastal risk analysis [16,49], and to predict the future behavior of coastal systems [52]. It stands to reason, then, that the BBNs approach is ideal to address numerous issues associated with management in the coastal zone, especially in the context of SLR and climate change scenarios [44].
Today, the rapid expansion in the number and availability of data from widely differing sources (especially from Earth observation) has resulted in advances in methodological approaches and computational tools and has created challenges for their effective and efficient processing and vulnerability or risk model implementation. This creates a dilemma as to how best to combine the data sets for maximum utility. This study gave the opportunity to update and expand coastal data sets for a specific study area, mainly based on multiple EO and ancillary data sets pertaining to vulnerability and risk assessment. The proposed approach seeks to increase the efficacy and efficiency of the remotely sensed data uptake for information extraction to support planners and managers in the design and application of more robust and adaptive measures. To deal with risks and uncertainties arising from incomplete knowledge about deltaic environment and its diverse ecosystems, the descriptive/predictive models enable us to consider both natural and human capital together, providing more comprehensive predictions of cumulative impacts. This is done by adopting a BBN-based model to integrate and synthesize a variety of data to derive vulnerability probability in a deltaic system and use it to test relationships between interacting natural and human variables. This testing should improve the assessment of vulnerability in an integrated model including natural and human sub-system interactions, thus providing an effective descriptive/predictive network. Taking into account the data availability and reliability for the identified study area, we developed models based on multisource EO data [53][54][55] to assess parameters influencing vulnerability. The model informs the BBNs with driver, land cover, and land use factors to produce probability maps for a quick scan and spatial visualization of coastal vulnerability and aims at enforcing coastal regulations and supporting management decisions to recognize high-priority intervention hot spots in natural and social subsystems.
The article is organized as follows: after introducing the current methods to assess coastal vulnerability and applications, Section 2 describes the case study, vulnerability assessment (conceptual approach and procedure), the input data used and the steps followed for the definition of key variables, the BBN-based modelling approach, the testing scenarios, and the sensitivity analysis. Section 3 shows vulnerability probability maps for current and future scenarios of SLR, including the identification of hotspot areas and exploring the relative influence of considered variables in determining predicted vulnerability. Section 4 discusses the results of the analysis, data gaps, and future challenges, especially regarding the potential of EO and BBNs to assess vulnerability. Finally, Section 5 summarizes the derived conclusions.

Study Area: Po Delta
The study focuses on the Po River Delta, the largest delta in Italy, which, like other major world deltas, is highly threatened by SLR [56][57][58][59] and natural subsidence (i.e., tectonic, sediment loading and sediment compaction) enhanced by anthropic activities (i.e., groundwater and methane extraction, drainage, and intensive farming practices) [60], making this coastal environment greatly vulnerable to inundation and erosion [61]. One of the main features of this area in terms of geomorphological, economic, and cultural value is that the Po Delta (Figure 1a) [62] underwent many changes of shape [63] and uses throughout history but still represents one of the most productive and biodiversity-rich areas. The system of the Po Delta is a relatively young (~4000-5000 years) alluvial fan that progrades from the Pianura Padana to the Adriatic Sea. The Po Delta covers an area of about 400 km 2 that extends seaward for about 25 km and is divided into five active branches. It is characterized by many dune systems and barrier islands that are sometimes connected to spits that delimiting lagoon areas along its coastline. Taking into account its natural variability and complexity, it is considered to be one of the most important European natural areas. Accordingly, it encompasses areas of great natural beauty and monuments of historical interest and was included on the UNESCO World Heritage list in 1995. The actual landform of the Po Delta is the result of river sedimentation processes but also of manmade actions that, over the centuries, have regulated the waters and reclaimed the land. It is an asymmetric delta dominated by shore-parallel, southward advection of finer-grained sediments. This asymmetric setting is also due to the loss of hydraulic efficiency of the northern fluvial branches in favor to the three southern branches. Until the middle of the 20th century, progradation of the delta was noticeable due to the abundant sediment supply [63]. Moreover, in the past few decades, Water 2020, 12, 2830 5 of 30 the natural rates of subsidence (i.e., tectonic, sediment loading, compaction of Holocene sediments) estimated to be 2-4 mm/year were enhanced by human activities (i.e., groundwater and methane extraction, sediment supplies, drainage, and intensive farming practices) and experienced as much as 3 m of land subsidence from the 1930s to the 1970s [60]. Elevated subsidence rates, decreasing sediment supply due to dam and reservoir interception, and sediment mining in the lower river course caused major delta erosion after the 1950s. Indeed, over the last 50 years, coastal erosion affects both the emerged and the submerged (prodelta) systems [64], and a marked reduction of sediment supply led to a fast degradation and shoreline retreat [65] accelerated by human activity [42]. Most significant changes in terms of land use and land cover evolution is represented by a decrease in saltmarshes and a broadening of agricultural areas [66]. Additionally, over the last century, an artificial surface significantly spread out, and natural swamp and salt marsh areas decreased by 85% [66]. These and other interventions have deeply modified the physical and ecological characteristics of the Po Delta, determining a drastic change in existing ecosystems, together with a general decrease in the standards of environmental quality. Furthermore, the delta evidences multiple threats that will probably be exacerbated in the following decades by the effects of climate change (acting on wave dynamics and sea level rise), with potential environmental and socioeconomic impacts. Natural gas withdrawal led to high subsidence rates in the Po Delta [67], and large areas of the delta are now more than two meters below sea level, which is expected to have significant impacts related to increasingly SLR throughout the 21st century [13,33].
Water 2020, 12, x FOR PEER REVIEW 5 of 31 Indeed, over the last 50 years, coastal erosion affects both the emerged and the submerged (prodelta) systems [64], and a marked reduction of sediment supply led to a fast degradation and shoreline retreat [65] accelerated by human activity [42]. Most significant changes in terms of land use and land cover evolution is represented by a decrease in saltmarshes and a broadening of agricultural areas [66]. Additionally, over the last century, an artificial surface significantly spread out, and natural swamp and salt marsh areas decreased by 85% [66]. These and other interventions have deeply modified the physical and ecological characteristics of the Po Delta, determining a drastic change in existing ecosystems, together with a general decrease in the standards of environmental quality. Furthermore, the delta evidences multiple threats that will probably be exacerbated in the following decades by the effects of climate change (acting on wave dynamics and sea level rise), with potential environmental and socioeconomic impacts. Natural gas withdrawal led to high subsidence rates in the Po Delta [67], and large areas of the delta are now more than two meters below sea level, which is expected to have significant impacts related to increasingly SLR throughout the 21st century [13,33]. Observed SLR trends from tide gauge data between 2010 and 2019 also showed a general rise in the Adriatic Sea level, which ranges from 6.7 cm to 7.2 cm [68][69][70][71]. This situation is particularly alarming for the coastal restoration, mostly in the areas around the Po River Delta [9,33,66]. Indeed, the SLR [68,69] expected over the next 40-50 years will accentuate the transgressive phenomenon.
Together with the aforementioned aspects making this region particularly vulnerable to coastal erosion and flooding, anthropic pressures (e.g., tourism, industry, and fishery) and the effect of climate change (acting on wave dynamics and sea-level rise) are expected to affect the mobilization of sediments along the shoreline in the next few decades, with potential environmental and socioeconomic impacts [13]. At the delta-plain scale, many studies have already evaluated the Observed SLR trends from tide gauge data between 2010 and 2019 also showed a general rise in the Adriatic Sea level, which ranges from 6.7 cm to 7.2 cm [68][69][70][71]. This situation is particularly alarming for the coastal restoration, mostly in the areas around the Po River Delta [9,33,66]. Indeed, the SLR [68,69] expected over the next 40-50 years will accentuate the transgressive phenomenon.
Together with the aforementioned aspects making this region particularly vulnerable to coastal erosion and flooding, anthropic pressures (e.g., tourism, industry, and fishery) and the effect of climate change (acting on wave dynamics and sea-level rise) are expected to affect the mobilization of sediments along the shoreline in the next few decades, with potential environmental and socioeconomic impacts [13]. At the delta-plain scale, many studies have already evaluated the potential impacts and vulnerability aspects, including the consideration of most relevant deltaic processes [19,72], among them natural and anthropogenic subsidence rates [33,60,66,73] and related relative sea level rise (RSLR) [33,59,73,74], changes in the land use [66], and coastal erosion processes [9,13]. Despite the fact that it has been demonstrated that the Po River has resumed delta progradation since 2010, this has occurred especially in its northern portion [65], highlighting the fact that the analysis of spatial and temporal vulnerability variations is required, since processes are spatially non-uniform.

Vulnerability Assessment
Coastal vulnerability assessment is a capable method for accounting for the interaction of the multiple environmental threats that determine a delta's exposure [75,76]. BBN models have been chosen as a proper method for scenario simulations because they can assimilate different data distributions, including categorical, discrete, and continuous data.
Indeed, once BBN is trained, it is able to answer the following question: "What is the probability of target class (i.e., vulnerability), given the input distribution for the considered parameters under diverse sea level rise scenarios?" Once the probability that answers this question is calculated for every spatial units (cells) in the entire study area, a GIS-based probable vulnerability map is generated to produce belief maps based on the exploratory analysis of BBNs.
To characterize the vulnerability conditions of the study area by merging different kinds of variables from multiple data sources, a procedure made of several methodological phases has been developed ( Figure 2). These phases include (1) data input, what includes definition, categorization, and homogenization of the variables; (2) BBN model design, which includes the organization of conditional probabilities between variables; and (3) vulnerability assessment, which includes the testing scenarios and the sensitive analysis.
Water 2020, 12, x FOR PEER REVIEW 6 of 31 relative sea level rise (RSLR) [33,59,73,74], changes in the land use [66], and coastal erosion processes [9,13]. Despite the fact that it has been demonstrated that the Po River has resumed delta progradation since 2010, this has occurred especially in its northern portion [65], highlighting the fact that the analysis of spatial and temporal vulnerability variations is required, since processes are spatially non-uniform.

Vulnerability Assessment
Coastal vulnerability assessment is a capable method for accounting for the interaction of the multiple environmental threats that determine a delta's exposure [75,76]. BBN models have been chosen as a proper method for scenario simulations because they can assimilate different data distributions, including categorical, discrete, and continuous data.
Indeed, once BBN is trained, it is able to answer the following question: "What is the probability of target class (i.e., vulnerability), given the input distribution for the considered parameters under diverse sea level rise scenarios?" Once the probability that answers this question is calculated for every spatial units (cells) in the entire study area, a GIS-based probable vulnerability map is generated to produce belief maps based on the exploratory analysis of BBNs.
To characterize the vulnerability conditions of the study area by merging different kinds of variables from multiple data sources, a procedure made of several methodological phases has been developed ( Figure 2). These phases include (1) data input, what includes definition, categorization, and homogenization of the variables; (2) BBN model design, which includes the organization of conditional probabilities between variables; and (3) vulnerability assessment, which includes the testing scenarios and the sensitive analysis. First, all layers were clipped on the Po Delta study area, considering only natural and social subsystems, for removing data outside the purpose of the investigation (i.e., artificial surfaces and agricultural areas). The dimension of the spatial units were selected at the beginning of the assessment, based on the purpose of the analysis and on the spatial resolution of the available data. Our aim was to generate assessment units that will respond uniformly in terms of probable  First, all layers were clipped on the Po Delta study area, considering only natural and social subsystems, for removing data outside the purpose of the investigation (i.e., artificial surfaces and agricultural areas). The dimension of the spatial units were selected at the beginning of the assessment, based on the purpose of the analysis and on the spatial resolution of the available data. Our aim was to generate assessment units that will respond uniformly in terms of probable vulnerability and can be treated as single units for the purpose of adaptation or intervention planning in the future. We created a set of spatial coastal assessment units considering the radius of influence of coastal erosion (RICE area) as the landward limit [9,13,30]. According to the EUROSION definition, RICE is represented by a buffer of 500 m from the coastline extended to areas lying below −5 m [35,57], and it is considered as a proxy of the terrestrial areas, which may potentially be subject to coastal erosion or flooding in the coming period of 100 years. The threshold representing the RICE was set equal to 1 km inland on the Po Delta, as proposed in Gallina et al. [13] for the northern Adriatic area, covering a total delta surface of 64 km 2 and 62 km of the coastline. Likewise, the delta is split into three macro-sectors ( Figure 3a) considering diverse coast orientation based on main wave climate [77] as the distribution of wave height, period, and direction averaged over a period of time (i.e., 2007-2017). Each macro-sector was divided into homogenous 500-m-long segments along the coastline (Figure 3b,c). Overall, there are 115 cells with an average length of 500 m and a width of 1 km, which have been overlaid with vulnerability assessment input data. To test the network, we explored several cases where we constrained the model with specific SLR scenarios, evaluating the probability of vulnerability, given a specified 100% probability that the rate of SLR refers exactly to each of the four diverse scenarios as described in Section 2.5 ( Table 2).
Water 2020, 12, x FOR PEER REVIEW 7 of 31 RICE is represented by a buffer of 500 m from the coastline extended to areas lying below −5 m [35,57], and it is considered as a proxy of the terrestrial areas, which may potentially be subject to coastal erosion or flooding in the coming period of 100 years. The threshold representing the RICE was set equal to 1 km inland on the Po Delta, as proposed in Gallina et al. [13] for the northern Adriatic area, covering a total delta surface of 64 km 2 and 62 km of the coastline. Likewise, the delta is split into three macro-sectors ( Figure 3a) considering diverse coast orientation based on main wave climate [77] as the distribution of wave height, period, and direction averaged over a period of time (i.e., 2007-2017). Each macro-sector was divided into homogenous 500-m-long segments along the coastline (Figure 3b,c). Overall, there are 115 cells with an average length of 500 m and a width of 1 km, which have been overlaid with vulnerability assessment input data. To test the network, we explored several cases where we constrained the model with specific SLR scenarios, evaluating the probability of vulnerability, given a specified 100% probability that the rate of SLR refers exactly to each of the four diverse scenarios as described in Section 2.5 ( Table 2).

Data Input
The identification of datasets counts on a broad availability of EO data processed by the authors and of already processed products from the European Copernicus Land Monitoring Service (CLMS) and the Copernicus Marine Environment Monitoring Service (CMEMS). This study takes into account other data sources including a variety of physical, ecological, and social variables to characterize the spatial pattern and the frequency distribution of the coastal vulnerability model.
Twelve variables are identified as relevant for simulating deltaic vulnerability scenarios, bearing in mind that delta systems are rapidly developing areas with changing landforms, sediment deposition rates, and vegetation cover. The physical, ecological, and social variables have been grouped in three broad classes representing the sources (driver factors-DFs), the pathway (land cover factors-LCFs), and the receptors (land use factors-LUFs) of the vulnerability [78,79] (respectively colored in red, green, and blue in Table 1).
It was assumed that the susceptibility to erosion decreases if the soil has high vegetative cover that contributes to higher wave energy dissipation. In response to rising sea levels and subsidence and diminishing fluvial sediment discharge [4], most deltas would naturally reduce their size under wave forces and migrate to shallow parts of the basin by switching or inundation. The deltaic fringe

Data Input
The identification of datasets counts on a broad availability of EO data processed by the authors and of already processed products from the European Copernicus Land Monitoring Service (CLMS) and the Copernicus Marine Environment Monitoring Service (CMEMS). This study takes into account other data sources including a variety of physical, ecological, and social variables to characterize the spatial pattern and the frequency distribution of the coastal vulnerability model.
Twelve variables are identified as relevant for simulating deltaic vulnerability scenarios, bearing in mind that delta systems are rapidly developing areas with changing landforms, sediment deposition rates, and vegetation cover. The physical, ecological, and social variables have been grouped in three broad classes representing the sources (driver factors-DFs), the pathway (land cover factors-LCFs), and the receptors (land use factors-LUFs) of the vulnerability [78,79] (respectively colored in red, green, and blue in Table 1).
It was assumed that the susceptibility to erosion decreases if the soil has high vegetative cover that contributes to higher wave energy dissipation. In response to rising sea levels and subsidence and diminishing fluvial sediment discharge [4], most deltas would naturally reduce their size under wave forces and migrate to shallow parts of the basin by switching or inundation. The deltaic fringe (i.e., subaerial and subaqueous parts of the deltaic coast interacting with and being modified by waves, tides, and currents) responds by barrier and dune buildup and sediment redistribution. Areas with less vegetation coverage were then assumed to be more susceptible (less resilient) to coastal erosion than wider vegetated areas [13].
In Table 1, the 12 variables are listed and described, highlighting data source types and formats, spatial resolutions, and temporal coverage. The use of multitemporal data acquired allowed for exploring diverse surface physical processes involved in the Po Delta dynamics. The use of an EO time series dataset from seven to 30-year timespans ( Table 1) depends on the observed phenomena and on the objective of the data analysis. The case of vegetation and sediment dynamics led the use of a 30-year timespan for the analysis and classification of images considering the annual and seasonal change detection goal. Specifically, it allows us to identify land cover transitions and confirms the assumption of a high rate of changing conditions. All the continuous, raster variables are harmonized by a spatial resampling to obtain 100-m-cell grid resolutions and, in addition to the vector data, have been clipped on the reference study area buffer (RICE area).
Among the 12 variables, the vertical velocity was obtained using an EO radar dataset coming from the Envisat ASAR satellite sensor (2003-2010) with the small baseline subset (SBAS) technique and kriging interpolation. The variable is determined by the combination of ascending and descending line of sight (LOS) velocities. The vertical displacements velocity map is derived by interpolation (kriging) of points for a regular 100-m-cell grid.
The other large processing was the EO optical dataset based on Landsat and Sentinel 2 satellite datasets (1990-2019) with the linear spectral mixture analyses (LSMA) technique [80,81] (Tables S1 and S2). A 30-year stack of images derived from the two satellite sensors, geometrically and atmospherically corrected, was first masked to keep only the coastal land and the very shallow waters and discard sea and lagoon areas that surround the barrier islands. The LSMA processing to retrieve from the mixed reflectance signal of each pixel the fractional coverage abundance of sediments, vegetation, and water (dark) was based on a three multitemporal endmembers and fully constrained spectral mixture model on a regular 30-m-cell grid, as done in Filipponi et al. [81]. The derived variables called coverage factor and protective distance stand for the vegetation and sediment fractional abundance and are used as cover percentage maps to provide proxies for wave energy dissipation factor and coastal erosion protection [25]. In particular, the protective distance variable is obtained by calculating the empirical orthogonal function (EOF) of the 30-year layer stack, as done in Filipponi et al. [82], to obtain the temporal signal decomposition into dominant variability and to extract the most common variations sediment spatial patterns and temporal trends. The EOF technique was used to discriminate seasonal changes from the inter annuals and the alongshore averaged value of the maximum sediment variation in the last 30 years was considered as the potential protective distance acting as a defense buffer against erosion processes [83].
In order to have a complete source-pathway-receptor-consequence chain [54], the action of offshore waves were propagated down to each of the three coastal orientation and 115 cells (Figure 3c) to generate a point matrix with the frequency distribution of wave height, period, and direction averaged over the period of the analysis (2007-2017).
The waves were used in conjunction with the DEM to obtain the impact scale classification as in Sallenger et al. [14] to define four regimes representing different levels of impact and thresholds defining where processes and magnitudes of impacts change dramatically. The regime is calculated using the general expression Equation (1) for beach run up parameterization [15] applied to the Water 2020, 12, 2830 9 of 30 115 transects orthogonal to the shoreline spaced of 500 m, representing topographic variability along the coastline for each defined segment: where H 0 = wave height; L 0 = deep-water wave length; β f = beach steepness; and T p = wave period calculated using the empirical formula Equation (2) [84][85][86]: where H s = significant wave height.
In accordance with Sallenget et al. [14], susceptibility to erosion (and consequently vulnerability) is higher for "inundation" regime where surge is sufficient to completely and continuously submerge the barrier islands. In the case of "overwash" regime, wave runup overtops the berm or, if present, the foredune ridge, and the associated net landward sand transport contributes to net migration of the barrier landward. Runup exceeds the threshold of the base of the foredune ridge in the "collision" regime, while in the "swash" regime runup is confined to the foreshore.
The adaptive capacity of the social systems depends on space and functioning organization, social and economic development, government and governance structure, and the presence of institutions that collect and store knowledge and experiences [87]. Here, adaptive capacity of the Po Delta social subsystem is posed by the resilience index [88], involving socio-economic factors that affect community resilience at the municipality scale, this being the primary administrative unit in case of emergencies and strategic in the urban planning activities. Specifically, it defines the metrics for community resilience drivers customized for the Po Delta from stakeholders' interviews and states the variables able to measure social, economic, community capacity, infrastructure, and environment resilience subcomponents, proposing a list of factors indicated by the literature as appropriate for the resilience concept clustered into relevant dimensions. Table 1. List of input variables used in the Bayesian belief network (BBN) and description.

Variables (and Periods) Description Data Type Processing Method and Data Source
Elevation (2007) The elevation of the area is partially below mean sea level (msl) due to subsidence induced by hydraulic engineering works, extraction of gas-bearing water, and reduction of sediment supply and embanking.  The mean value of the maximum sediment spatial variations in the last 30 years is considered as the potential protective distance acting as a defense buffer against erosion processes. It is related to the dune presence and to the deltaic small-scale coastline changes.
Raster of the alongshore averaged values of the coastline sediment variation in the last 30 years, continuous, terrestrial, 30-m-cell grid.
Vectorial shore classified as geomorpho-sedimentological types, discrete, terrestrial-intertidal. The spatial and temporal variability of global surface water presence and its long-term changes (i.e., 32 years) calculated on 3 million Landsat images [90].
Protected natural area The network of core natural habitat types represents the natural capital of the area and corresponds to high exposure levels.
Vectorial polygons classified as protection designation, discrete, terrestrial.
It is the Natura 2000 network based on the Habitat and Birds Directives reporting and designation. Data source: CLMS.

Resilience Index (2017)
The resilience index is referred to the adaptive capacity of the Po Delta social subsystem and it involves socio-economic factors that affect community resilience at municipality scale.
Vectorial points at municipality scale, classified as degree of resilience, discrete, terrestrial.
Metrics obtained from stakeholders' interviews, to score infrastructure and environment resilience's subcomponents. Data source: [88].
* Processed by the authors.

BBN Model Design
BBN models based on probabilistic structure are best suited to perform scenario analysis such as quantitative vulnerability and risk analysis, error propagation, and decision making, while traditional deterministic modelling cannot provide the sophisticated spatial analytics required to model uncertainty within the natural environment. There are four main benefits to using BBNs compare to other methods. The ability to include uncertain or unavailable input variables coupled with expert knowledge makes BBNs a very powerful tool in decision support systems. Second, they describe relationships of causes and effects using a graphical framework (Phase 2) that provides for quantification of probabilities and clear communication of the results among scientists, managers, and the broader community, improving the integration of probabilistic modelling into adaptive management practices. Third, BNNs incorporate an explicit causal structure having more variables that interact directly and indirectly (Phases 3 and 4). Finally, BBNs perform computationally fast uncertainty assessments in complex systems and therefore are ideally suited for a spatial vulnerability trends distribution modelling, improving models' predictive capability, where results are consistent if the aim is the definition of the maximum and minimum conditions of vulnerability (Phase 5).
It is important to be aware of the system complexity when building a Bayesian belief network (BBN); thus, constructing a model network first involves defining the purpose of the modelling and the variables that best describe the process to be studied in order to identify the relationship between the variables (which variables are impacted by which variables) and illustrating the relevance of some factors over others. The model network was designed with the 12 selected variables to address the most probable vulnerable areas following difference phases (Figure 2).  Figure 4). Based on the discrete quantities, a score from 1 to 4 contributing in increasing the final probable vulnerability has been assigned.
Water 2020, 12, x FOR PEER REVIEW 12 of 31 on the discrete quantities, a score from 1 to 4 contributing in increasing the final probable vulnerability has been assigned. • Phase 2. A generic conceptual model ( Figure 5) was used to develop an integrated network able to represent the relationships between input variables and the interconnection with sourcepathway-receptor-consequence (SPRC) model applied to the vulnerability assessment. BBN is a directed acyclic graph (DAG) where the selected variables are represented by nodes and the arrows indicate the direct cause-effect relationships between the nodes [91] ( Figure 5). The network is solved when nodes have been updated using Bayes' Rule Equation (3): where P(A) is the prior distribution of parameter A; P(A|B) is the posterior distribution, the probability of A given new data B; and P(B|A) the likelihood function, the probability of B given existing data A. Our hypothesis is that the functional relationship can be learned from available data set and cast in terms of a probabilistic estimate in the form Equation (4): The conditional probability distribution of the Vulnerability (V), that is the "child node," to a particular set of inputs variables can be obtained from the joint probability of driver factors (DFs), land cover factors (LCFs), and land use factors (LUFs) that are the "parents nodes". • Phase 2. A generic conceptual model ( Figure 5) was used to develop an integrated network able to represent the relationships between input variables and the interconnection with source-pathway-receptor-consequence (SPRC) model applied to the vulnerability assessment. BBN is a directed acyclic graph (DAG) where the selected variables are represented by nodes and the arrows indicate the direct cause-effect relationships between the nodes [91] ( Figure 5). The network is solved when nodes have been updated using Bayes' Rule Equation (3): where P(A) is the prior distribution of parameter A; P(A|B) is the posterior distribution, the probability of A given new data B; and P(B|A) the likelihood function, the probability of B given existing data A. Our hypothesis is that the functional relationship can be learned from available data set and cast in terms of a probabilistic estimate in the form Equation (4): where P(A) is the prior distribution of parameter A; P(A|B) is the posterior distribution, the probability of A given new data B; and P(B|A) the likelihood function, the probability of B given existing data A. Our hypothesis is that the functional relationship can be learned from available data set and cast in terms of a probabilistic estimate in the form Equation (4): The conditional probability distribution of the Vulnerability (V), that is the "child node," to a particular set of inputs variables can be obtained from the joint probability of driver factors (DFs), land cover factors (LCFs), and land use factors (LUFs) that are the "parents nodes".  [78,79]. BBN conceptual model Figure 5. Bayesian belief network (BBN) conceptual model and associated source-pathway-receptorconsequence (SPRC_chain applied to the vulnerability assessment [78,79]. BBN conceptual model indicating that nodes driver factors (DFs), land cover factors (LCFs), and land use factors (LUFs) have a common effect on node V. This implies a conditional dependence between nodes DF, LCF, and LCF that are in effect competing to explain V [92]. Solid line circle indicates parent nodes, and the colored solid arcs indicate that all variables of the parent category influence all variables of the child category; thick line circle represents target node.
The conditional probability distribution of the Vulnerability (V), that is the "child node," to a particular set of inputs variables can be obtained from the joint probability of driver factors (DFs), land cover factors (LCFs), and land use factors (LUFs) that are the "parents nodes".
• Phase 3. To construct the BBN, the GeNie Academic Version 2.4 software package (BayesFusion, LCC, Pittsburgh, PA, US) was chosen. Classes corresponding to assigned score categories as described above resolve each parent or child node. The connection of each node in the BBN is determined by the conditional probability table (CPT) and probability table (PT), where the degree of belief (probability) that a variable will be in a particular state given the state of the parent variable [93] is defined. • Phase 4. The entries in the CPTs and PTs can be discovered using a range and combination of methods, including directly observed data and probabilistic or empirical equations, results from model simulations and elicitation from expert knowledge. Some of the parameter learning functions are also based on expectation maximization (EM) algorithm (i.e., it requires only the model structure to be known beforehand and iteratively calculates maximum likelihood estimates for the parameters given the data and the model structure [38]) available in GeNie [94], implemented to provide prior probability for root nodes, while target node conditional probabilities is reported as a specific matrix representing the frequency of vulnerability scores associated with each input variables from one to four per cell. Vulnerability is selected as target node, the output variable that directly represents the phenomenon [92] ( Figure 5). The prior probability is provided for regime, vertical velocity, dune, protective distance, coverage factor, global water surface dynamic, protected natural areas, and resilience index. • Phase 5. Accordingly, the output of running the network contains the posterior probability distribution of the vulnerability node implementing extracted and reclassified values from resampling a GIS-based data set at 100 × 100 m resolution.

Vulnerability Assessment
The optimized model was used to predict vulnerability probability for the present SLR scenario and then for the three future SLR scenarios [56]. Once parameters have been learnt, a prediction can be made by updating the BBN with "evidence," whereby the states of SLR node in the systems assume certain probabilities and the resultant conditional probabilities of the target variable are estimated. Vulnerability worst-case projection of global mean SLR to 100 years [56] (the highly pessimistic) has been considered (i.e., the more precautionary conditions). Three SLR future scenarios have been finally implemented in the BBN as evidence (Table 2): conservative scenario, pessimistic scenario, and highly pessimistic scenario. In the testing phase, many factors such as changes in wave height are characterized by high uncertainty in the northern Adriatic [57,58]. Except for the value of SLR, the values of all the other variables are assumed to be unchanged in hundred years, even average vertical velocity and surface water dynamics are considered in accordance with the present trends. Starting from posterior vulnerability probability obtained from BBNs models testing scenarios, the hotspot areas of the study coast are identified.

Sensitive Analysis
Broadly, in a Bayesian network, the sensitivity analysis refers to the analysis of the influence and influence degrees of multiple causes (node states) on the result (target node) [38,95]. The sensitivity analysis function of GeNie 2.4 can be used to identify which variables in the system have the highest influence on the outcome of the model. This is done by investigating the effect of small changes in numerical parameters (i.e., probabilities) on the output parameters (e.g., posterior probabilities). Highly sensitive parameters affect the reasoning results more significantly having the highest influence on the delta vulnerability target node. GeNIe implements an algorithm proposed by Gaag et al. [96] that performs simple sensitivity analysis in Bayesian networks. The algorithm efficiently calculates a complete set of derivatives (i.e., the value of the first derivative of the posterior probability of the selected state of the target node at its current value over the considered parameter) of the posterior probability distributions over the target node over each of the numerical parameters of the Bayesian network. These derivatives give an indication of importance of precision of network numerical parameters for calculating the posterior probabilities of the target. If the derivative is large for a parameter, then a small change in that parameter may lead to a large change in the posteriors of the target variable. If the derivative is small, then even large changes in the parameter make little difference in the posteriors. The tornado plot of the sensitivity test shows the most sensitive parameters for a selected state of the target node sorted from the most to least sensitive.
Strength of influence is calculated from the CPT of the child node and essentially expresses some form of distance between the various conditional probability distributions over the child node conditional on the states of the parent node taking the plain average over distances in the normalized mode.

Data Input Discretization
Input variables have been transformed in categorical scales, and a score, ranging from 1 (low vulnerability) to 4 (high vulnerability), was assigned to each variables' category (Table 3). Table 3. Summary of variables included in the BBN, discretization methods applied and distribution classes, intervals, scores and states associated with the input variables identified in the vulnerability matrix for BBN model in the Po Delta area colored based on the types of variables representing the sources, driver factor (red), the pathway, land cover factor (green), the receptors, land use factor (blue). A score of 1 means a very low contribution to vulnerability, and a value of 4 means a high contribution to vulnerability.

Interval Score State
Expert decision/literature review Elevation above sea level represent the 38% of emerged areas, mostly concentrated in the northern macro-sector. The central and southern macro-sectors emerged areas are mostly below sea level reaching elevation values of −3 m above sea level in some specific sectors.    Table 3. Cont.

Interval Score State
Equal intervals/literature review In the northern macro-sector 9.7 km 2 could be characterized by the "collision" regime, whilst 6.4 km 2 could be characterized by the "overwash" regime. In the central part of the delta "collision" regime is more probable north of Po di Pila mouth, while the "overwash" regime is expected in the sector adjacent and south of Po di Pila mouth. In the southern part of the study area 'collision' regime is more probable, while along the coastline located close to Po di Tolle mouth, the occurrence of the "overwash" regime is expected. The cover map based on the per-pixel vegetation abundance obtained with the LSMA, the average value of the vegetation abundance (AVA) in % of the pixels that each polygon contains has been calculated, obtaining 45% as the maximum AVA. The final cover map differs between vegetated and non-vegetated emerged areas, depending on the AVA being higher or lower than the 50th percentile. Mapped vegetated areas are mostly distributed close to Po di Mistra, Po di Pila, and Po di Tolle and south of Po di Goro.

Input Variables Discretization Method and Main Distribution Results
Interval Score State

Equal intervals
The potential protective distance acting as a defense buffer against erosion processes obtained using the EOF approach, divided in four classes, from less than 40 m up to 1500 m. The highest values are found along the coastline, between Po di Mistra and Po di Pila mouths, and in the southern part, close to Lido di Volano. In the central macro-sector in correspondence of Po di Pila mouth, a loss of permanent water is observed as a consequence of conversion of permanent water into land. Rather, in the southern macro-sector close to the Po di Gnocca mouth land converted into seasonal water, while in correspondence of the sandy spit (the Scanno in Sacca di Goro) located south of Po di Goro mouth a loss of permanent water occurred.

BBN Model
To address the question of how vulnerable the Po Delta sub-systems are, we provided a probabilistic prediction giving the (in) dependence relationships among the variables and between variables and the target node, designing a BBN graphical model using a DAG.
The resulting graphical network (Figure 6), derived from the hypothesis that the vulnerability is a function of 12 input variables having 1537 possible combinations of input variable states. The target node conditional probabilities matrix represents the frequency of vulnerability scores associated with each input variables with values from one to four per cell.
Water 2020, 12, x FOR PEER REVIEW 18 of 31

BBN Model
To address the question of how vulnerable the Po Delta sub-systems are, we provided a probabilistic prediction giving the (in) dependence relationships among the variables and between variables and the target node, designing a BBN graphical model using a DAG.
The resulting graphical network (Figure 6), derived from the hypothesis that the vulnerability is a function of 12 input variables having 1537 possible combinations of input variable states. The target node conditional probabilities matrix represents the frequency of vulnerability scores associated with each input variables with values from one to four per cell. Figure 6. BBNs developed to model vulnerability in the Po River Delta. The nodes are grouped and colored based on the types of variables representing the sources, driver factor (red), the pathway, land cover factor (green), the receptors, land use factor (blue), and the consequence target node (beige). Arrows denote the direction of probabilistic influence implemented in software rather than causal relationships between the variables. Figure S1 reports two of the four testing scenarios applied to the whole delta for present and highly pessimistic scenarios.

Testing Scenarios
Focusing on reference present scenario at macro-sector scale (Figure 7), probability of "high" vulnerability accounts for the 7%, 22%, and over the 7% for the northern, central and southern macrosectors, respectively. The probability of "moderate" vulnerability constitutes the 22% of both central and southern macro-sectors, and the 36% of the northern. Looking at "low" vulnerability, the variability of probability among macro-sectors is smaller, being 24%, 22%, and 22% for the northern, central, and southern areas, respectively. "Very low" vulnerability shows probability values equal to 33% for all the macro-sectors. Focusing on the highly pessimistic scenario at macro-sector scale a Figure 6. BBNs developed to model vulnerability in the Po River Delta. The nodes are grouped and colored based on the types of variables representing the sources, driver factor (red), the pathway, land cover factor (green), the receptors, land use factor (blue), and the consequence target node (beige). Arrows denote the direction of probabilistic influence implemented in software rather than causal relationships between the variables. Figure S1 reports two of the four testing scenarios applied to the whole delta for present and highly pessimistic scenarios.

Testing Scenarios
Focusing on reference present scenario at macro-sector scale (Figure 7), probability of "high" vulnerability accounts for the 7%, 22%, and over the 7% for the northern, central and southern macro-sectors, respectively. The probability of "moderate" vulnerability constitutes the 22% of both central and southern macro-sectors, and the 36% of the northern. Looking at "low" vulnerability, the variability of probability among macro-sectors is smaller, being 24%, 22%, and 22% for the northern, central, and southern areas, respectively. "Very low" vulnerability shows probability values equal to 33% for all the macro-sectors. Focusing on the highly pessimistic scenario at macro-sector scale a growing probability of "high" vulnerability is evident, increasing at 24%, 38%, and 24% for the northern, central, and southern macro-sectors, respectively, to the detriment of the "moderate" vulnerability probability category. growing probability of "high" vulnerability is evident, increasing at 24%, 38%, and 24% for the northern, central, and southern macro-sectors, respectively, to the detriment of the "moderate" vulnerability probability category.
The probability of "moderate" vulnerability noticeably increased in the pessimistic scenario resulting in 52%, 40%, and 54% for the northern, central, and southern macro-sectors, respectively, to the detriment of the "low" vulnerability probability category. The characteristic behaviors of the vulnerabilities for each SLR scenario along each of the 115 cells divided in three macro-sectors is represented in Figure 8. As expected, the results show mostly a constant increase of 16-17% of probability of "high" vulnerability comparing the present and highly pessimistic scenarios. Taking as a reference the present scenario and looking at a conservative scenario projected SLR of 0.52 m, it appears that the overall vulnerability of the delta will not change considerably within 100 years, and "low" vulnerability probability increases at the expenses of "very low" vulnerability probability across all the three macro-sectors. Looking at a pessimistic scenario projected SLR of 0.74 m, the "moderate" probability of vulnerability rises almost 17 percentage points compared to reference scenario reaching 54% in the southern sectors and to almost 16% in the other sectors.  The probability of "moderate" vulnerability noticeably increased in the pessimistic scenario resulting in 52%, 40%, and 54% for the northern, central, and southern macro-sectors, respectively, to the detriment of the "low" vulnerability probability category.
The characteristic behaviors of the vulnerabilities for each SLR scenario along each of the 115 cells divided in three macro-sectors is represented in Figure 8. As expected, the results show mostly a constant increase of 16-17% of probability of "high" vulnerability comparing the present and highly pessimistic scenarios. Taking as a reference the present scenario and looking at a conservative scenario projected SLR of 0.52 m, it appears that the overall vulnerability of the delta will not change considerably within 100 years, and "low" vulnerability probability increases at the expenses of "very low" vulnerability probability across all the three macro-sectors. Looking at a pessimistic scenario projected SLR of 0.74 m, the "moderate" probability of vulnerability rises almost 17 percentage points compared to reference scenario reaching 54% in the southern sectors and to almost 16% in the other sectors. Vulnerability probability maps have been created by transferring the probabilistic results of vulnerability for each SLR scenario into the GIS environment ( Figure 9) and by defining the level of probability in each cell identified in the Po Delta.  Vulnerability probability maps have been created by transferring the probabilistic results of vulnerability for each SLR scenario into the GIS environment ( Figure 9) and by defining the level of probability in each cell identified in the Po Delta.
Water 2020, 12, x FOR PEER REVIEW 21 of 31 Figure 9. Vulnerability probability map at cell spatial scale considering SLR increasing. The color palette is divided into four color groups corresponding to the vulnerability target classes from very low (blue), low (green), moderate (yellow), and high (red). Each vulnerability target class is divided in seven classes of increasing probability score from lighter to darker hues.

Identification of Hotspot Areas
The calculated vulnerability probability has allowed us to identify the hotspot areas. Figure 9. Vulnerability probability map at cell spatial scale considering SLR increasing. The color palette is divided into four color groups corresponding to the vulnerability target classes from very low (blue), low (green), moderate (yellow), and high (red). Each vulnerability target class is divided in seven classes of increasing probability score from lighter to darker hues.

Identification of Hotspot Areas
The calculated vulnerability probability has allowed us to identify the hotspot areas. The very low probability accounts for the northern sector and, as expected, the increasing SLR would reduce the low probability of vulnerability and even maintaining the spatial pattern due to the different exposure of the systems to the waves intensity and frequency.
It is of interest that the social resilience in the northern macro-sector is higher, and this is an influencing factor on the lower exposure of receptors (Figure 9a,e,i,m).
The low vulnerability probability has a heterogeneous pattern moving along the cells and, unless the central delta, the values never reach very high probability (>75%) (Figure 9b,f,j,n).
Moderate vulnerability probability highlights, as hot spots the southern Po Delta, where the sandy spits and barrier are the dominant typology of the coast, as with very relevant dune and wetlands habitat, determining higher receptors exposure due to the natural capital. The dunes and the saltmarshes do not reach elevations enabling a barrier role to the increasing SLR (Figure 9c,g,k,o).
High vulnerability probability is recurrent in the central delta that is the ebb area with a higher rate of changing conditions due the interaction between sediments and seawaters. Here, even the exposure is lower in respect to the northern and southern macro sectors, the elevation is the lower, the overwash regime dominates, and no adaptive factors are able to reduce the probability of vulnerability (Figure 9d,h,l,p).

Sensitivity Analysis
The results of sensitivity tests and the strength of influence for a "high" vulnerability state show that the elevation and the water occurrence dynamic affect the uncertainty of vulnerability the most ( Figure 10). We found that, among variable waves, H 0 has the highest impacts on vulnerability. The other parameters seem to have the same degree of influence on the target variable, with lower values for protective distance variable. To analyze the sensitivity under the highly pessimistic SLR scenario, we used a tornado diagram to investigate the impacts on the "high" vulnerability state (vulnerability = high). The analysis provides the range (10% of its current value up and 10% down) and the direction (negative or positive) of changes in the target state. The most sensitive combination of the variables is regime = low; protective distance = low; protected areas = no; resilience = moderate; SLR = highly pessimistic; H 0 frequency = moderate ( Figure S2). The very low probability accounts for the northern sector and, as expected, the increasing SLR would reduce the low probability of vulnerability and even maintaining the spatial pattern due to the different exposure of the systems to the waves intensity and frequency.
It is of interest that the social resilience in the northern macro-sector is higher, and this is an influencing factor on the lower exposure of receptors (Figure 9a,e,i,m).
The low vulnerability probability has a heterogeneous pattern moving along the cells and, unless the central delta, the values never reach very high probability (>75%) (Figure 9b,f,j,n).
Moderate vulnerability probability highlights, as hot spots the southern Po Delta, where the sandy spits and barrier are the dominant typology of the coast, as with very relevant dune and wetlands habitat, determining higher receptors exposure due to the natural capital. The dunes and the saltmarshes do not reach elevations enabling a barrier role to the increasing SLR (Figure 9c,g,k,o).
High vulnerability probability is recurrent in the central delta that is the ebb area with a higher rate of changing conditions due the interaction between sediments and seawaters. Here, even the exposure is lower in respect to the northern and southern macro sectors, the elevation is the lower, the overwash regime dominates, and no adaptive factors are able to reduce the probability of vulnerability (Figure 9d,h,l,p).

Sensitivity Analysis
The results of sensitivity tests and the strength of influence for a "high" vulnerability state show that the elevation and the water occurrence dynamic affect the uncertainty of vulnerability the most ( Figure 10). We found that, among variable waves, H0 has the highest impacts on vulnerability. The other parameters seem to have the same degree of influence on the target variable, with lower values for protective distance variable. To analyze the sensitivity under the highly pessimistic SLR scenario, we used a tornado diagram to investigate the impacts on the "high" vulnerability state (vulnerability = high). The analysis provides the range (10% of its current value up and 10% down) and the direction (negative or positive) of changes in the target state. The most sensitive combination of the variables is regime = low; protective distance = low; protected areas = no; resilience = moderate; SLR = highly pessimistic; H0 frequency = moderate ( Figure S2). Figure 10. The result of the tests on sensitivity (i.e., from higher sensitivity variables in vivid red to lower sensitivity variables in grey) and strength of influence for a "high vulnerability" state considering the present and highly pessimistic scenarios. The thickness of the arrows indicates the strength of influence between the nodes that they connect (Table S3 in supplementary material). The Figure 10. The result of the tests on sensitivity (i.e., from higher sensitivity variables in vivid red to lower sensitivity variables in grey) and strength of influence for a "high vulnerability" state considering the present and highly pessimistic scenarios. The thickness of the arrows indicates the strength of influence between the nodes that they connect (Table S3 in

Discussion
In light of increasing anthropogenic and climatic pressures and the uncertain future that climate change poses, especially related to SLR, there is an unmet need for information solutions and theories that can enable the effective implementation of risk mitigation measures in low-lying coastal areas [5]. The main objective of the analysis was to assess and map vulnerability with target classes (i.e., very low, low, moderate, or high vulnerability) to help decision makers in future planning and management actions, implementing coastal regulations, and recognizing high-priority intervention hot spots in natural and social subsystems. The peculiarity of the integrated approach described here lies in the extraction of vulnerability values using EO data directly implemented in a BBN model, with the key advantage of generating conditional probabilities on the past and present data and exploring several future SLR scenarios.

Data Input Relevance
With the aim to overcome previous model-driven tools' weaknesses related to the capability of representing vulnerability on a spatial basis and variables mutual interactions, EO-based methods are suitable to include the complex physical and ecological variability of the Po Deltas' patterns and processes [59,60,[65][66][67]. The use of EO from satellites and Copernicus information along with modeling plays a key role in the provision of spatial products representing land surface attributes in a continuous way determining an advantage for a valuable data-driven approach and in coastal processes understanding. First, the combination of satellite data coupled with spatially referenced ancillary measures, including vegetation, topography and bathymetries, climate, and geology, [4,23,24,40,53,54,81,82,95], has proved to be a cost-effective approach to classify and quantify coastal parameters [23,53,54,82,83] influencing the delta's vulnerability (Tables 1 and 3). Second, the use of EO data contributes to a better understanding and regular monitoring of landforms and coastal dynamics. Specifically, the use of SAR-based interferometry for the monitoring of land subsidence reveals the heterogeneity of the subsidence rates along the coastal system, and the spatial pattern of this variable can be used as a source of vulnerability as well as a determinant of uncertainties in the interpretation of the relative sea level rise scenarios effects [33,59]. The use of an EO optical dataset based on Landsat and Sentinel 2 satellite 30-year data time series (Tables S1 and S2) definitely supports expert knowledge to deeply understand surface physical processes related to vegetation and sediments dynamics that highly influence the susceptibility to erosion processes across scales.

BBN Model Development
The use of BBNs allows the possibility to integrate variables acquired from different sources. Phase 1 of the model design accommodates the limited ability of BBNs to deal with continuous parameters derived from EO data input [38,39] and for running the model over the discrete sectors and subsectors, discretized values (i.e., from 2 to 4 intervals) treated as categorical (Figure 4) (Table 3). This allowed us to build a meaningful BBN and to capture complex empirical distributions and non-linear relationships among the variables. The BBN-developed model ( Figure 5) is applied using current conditions, which define the baseline, using future SLR scenarios taking in consideration the Intergovernmental Panel on Climate Change SLR projections [56] to define a plausible future vulnerability probability occurrence. For each identified sector in the three macro-sectors, the values in the input data are used as evidence in the network, and inference is performed to obtain the posterior probability distribution of the target nodes (Figures 7 and 8). Then, the posterior distributions of the target nodes are written into a spatial file making the results a combination of each cell (Figure 9). The simplicity of the links between variables in a BBN means that it is possible to model very complex systems and also to model these kinds of data with a different number of state variables [39]. Generally, in the BBN results, weaknesses emerge in the case of the arbitrary definition of conditional probabilities when expert knowledge is not available or there is a lack of "evidence" in terms of information on the intensity and extent of the past hazardous events and on the consequent damage impacts [63,65].

Vulnerability Assessment and Key Findings
The aim of mapping coastal vulnerability required running the BBN in a spatially explicit way able to include the spatial variability of the main coastal features, and the knowledge on the delta system allowed developing a segmentation of the coastal delta adopting as reference spatial boundary for vulnerability assessment [13,30,31,35,57,58] (Figure 3b,c) the landward limit of the radius of influence of coastal erosion (RICE) area.
The ongoing climatic changes will further intensify the present intensive erosive phenomena related to the reduction of fluvial discharges and subsidence [62,63]. Here, the results confirmed the outcome provided by Tosi et al. [59] using a multi-band SAR methodology integrating COSMO-SkyMed and ALOS-PALSAR images, since the subsidence rates increase eastward toward the central tip of the delta, mainly composed by unconsolidated, highly compressible thick prodelta mud. Consequently, there is a greater probability of high vulnerability hot spots in the central delta (Figure 9d), increasing in view of the highly pessimistic SLR scenario (Figure 9p). Human activities (i.e., drainage and intensive farming practices in the upstream portion of the delta [66], methane groundwater extractions, sediment supplies) have caused a worsening of the situation increasing subsidence rates related to natural consolidation of recent compressible low-permeable deposits [60,63,66] that, in view of the expected SLR, will continue to occur for decades in the portions of the Po Delta formed over the last two to three centuries [60]. Moreover, subsidence in concomitance with the heavy reduction of the sediment load of the Po and enhanced erosive phenomena [60,63], as observable in the reduced protective distance of southern sector, result in the greater percentage of moderate vulnerability probability hot spots in the current scenario ( Figure 9) (Table 3). Here, the greater probability of moderate vulnerability is rather influenced by LUFs, and the low social resilience is particularly affected by low social and institutional resilience [88] and the extensive presence of protected natural areas.
Therefore, observing hot spots vulnerability variation in the whole delta, non-climatic environmental and human-induced (i.e., coverage factor, protective distance and subsidence) and socio-economic (i.e., resilience) variables turned out to be relevant in determining expected higher vulnerability ( Figure 10) as for hot spots located in the central ebb area (Figures 7-9). Indeed, hot spots placed in the central macro sector show higher rate of changing conditions due the interaction between sediments and seawaters [4,65]. Considering that the exposure is lower in respect to the northern and southern macro sectors, the elevation is the lower and no adaptive factors (i.e., here lower economic and institutional resilience components diminish the adaptive capacity) are able to reduce the probability of vulnerability (likewise confirmed by the sensitivity analysis shown in Figure S2), the use of the heterogeneity of the subsidence rates in the BBN analysis has been relevant.
Furthermore, the inclusion of the LCFs (i.e., coverage factor, dunes presence and geomorphology) ( Table 1) contributes to explain the lower probability of high vulnerability for the majority of northern hot spots, with some exceptions for the northward hot spots (Figures 8 and 9), where despite higher adaptive capacity and lower subsidence rates, exposure is significantly higher (i.e., the absence of dunes and sub-merged vegetation acting as proxy for wave energy dissipation, prevalence of overwash regime, lower protective distance and developed beaches < of 1 km) (Figure 4) (Table 3). Highly heterogeneity of vulnerability to SLR, resulting in significantly different probability values of contiguous cells, particularly visible in both "moderate" and "high" vulnerability maps (Figure 9), reflects the combined result of the spatialized variables derived from EO data. Although the sensitivity analysis ( Figure S2) shows that vulnerability is more influenced by low energetic but more frequent events (Table 3), the frequency of storm-related erosive marine processes in the last decade has remained substantially constant, showing only limited deviation within the normal climate variability especially in the northern portion [65] and highlighting the importance to include the interaction between waves and coastal morphology to define potential impacts (i.e., regime). The results could be enriched by integrating the model flow with real cases (i.e., introducing "evidence") and adding significant additional factors and their combined effects that are important on local scales (e.g., coastal structures, river discharge, total suspend matter, wetlands typology distribution, longshore transport, sediment budgets, hydrological modification and deforestation, changes in coastline, and coastal erosion) may significantly improve predictions of sites experiencing coastal vulnerability. In this view, introducing evidences and limitations by informing the CPTs with evidentiary data and outcome states could lead to an improvement of the effective variance reduction that mainly depends on the influence of each input variables to the network.
The 12 variables frequency distribution used to generate categories all across the stretch makes the cells only comparable to each other, for the results of the Po Delta study, and they are not comparable with other deltas of the world because the frequency distribution of values of other areas, even considering the same set of variables, would be significantly different.

Management Implications
We build on the current state of knowledge with an approach capable to support short to medium-term coastal management decisions in identifying key variables and their interactions [40], prioritizing areas that are vulnerable, and determining planning interventions considering natural and social subsystems, thus improving deterministic models, which remain arguable when applied to real cases [38,39]. The consideration of EO time series dataset and the involvement of relevant stakeholders address the re-thinking of existing policy and planning constraints and adaptation strategies. In interpreting model outcomes, BBN results confirm that this is a valuable tool for vulnerability assessment, since it can incorporate different typology of data-both quantitative (e.g., obtained from existing models, monitoring and from site-specific measurements) and qualitative-obtained mostly from expert stakeholder knowledge, providing a probability-based approach. In addition, the probabilistic and analytical structure of the variables network takes into account uncertainties, which are often high in vulnerability scenarios. The graphic representation ( Figure 6 and Figure S1) enables a simplified communication for stakeholders and decision makers about the interaction among variables to the probability of vulnerability because (i) it makes easy deducing the probabilities of different causes given the consequences ( Figure 10) and (ii) it provides fast responses to queries to decision makers, in contrast to simulation-based models that requires time and computational power [38]. Findings of the study put in evidence that the synergistic effects of various pressures, particularly human-induced processes, appear to exacerbate the increasing of high vulnerability in the perspective of a highly pessimistic SLR scenario, at least at the local scale. Consequently, the identification of hotspot areas and the preparation of a correct coastal vulnerability assessment to SLR need to take into consideration also non-climate drivers, in the perspective of a local-level model downscaling. The study also highlights the greatest strength in the use of the social-ecological components to address the natural and human capital exposed to future SLR. The resilience index in particular provided the stakeholders knowledge integration that is considered one of the most relevant inputs when building the BBN. The protected natural areas along with the coverage factor provided the inclusion of ecosystem services components as the variables that account for the naturalistic role of the delta and protection from erosion. Considering the CF as a key indicator, we can state that the socio-ecological resistance and resilience evaluation could be based on the interaction of the physical, biological, and social components.

Conclusions
From the obtained outcomes, we are able to draw the following main conclusions: • EO satellite data keeps proving its high potential for monitoring the Earth, easing spatial and temporal analysis in very complex and dynamic areas.

•
The implementation of EO data into the BBN model provides a statistical framework performing inference on a very large amount of different data that describes what is known about the physical, ecological, and social parameters of the system and, perhaps more importantly, what is known about the causal relationships between them.

•
The BBN approach, based on the definition of conditional probabilities, associated to the use of a GIS tool, provided the trend of distributions of vulnerability along the shore, which is especially useful for management and planning purposes.

•
The overall assessment reveals that almost the 35% of the Po Delta is highly vulnerable to the highly pessimistic SLR scenario. The high heterogeneity of vulnerability to SLR reflects the combined results of the selected variables: 29% of the northern cells, 100% of the central cells, and 23% of the southern cells fall into vulnerability target classes that range from 0.34% to >72% of probability to being highly vulnerable. Our findings highlight potential changes in the sediment budget considering the major vulnerability of the central ebb that is the gate between the river sediment supply and the costal sediment transport.

•
The outcomes highlight the strength of influence of each variable on the target and their evident dependence and interconnection. Wave climate mostly depending on events frequency, land subsidence, and changes in land cover, which are the more sensitive variables playing a key role in vulnerability to SLR.

•
The study mainly highlights the need (i) to update and expand coastal data sets, specifically concerning local scale environmental and human-induced processes, and (ii) to include both bio-geophysical and socio-economic aspects within delta analysis of the distribution and change in risk components, identifying vulnerable areas and communities, and where changes in hot spots may occur in the future. This can inform other analyses, such as the design of surveys and data collection, as well as identify policy needs and indicate where adaptation actions are likely to be required.

•
The findings of this study are also valuable for the further performance of a gap analysis between the available tools and operational user needs. Emerged results provide a contribution in term of completeness of the information for a clearer overview of the expected vulnerability pattern and key variables. These insights could be useful to assess and verify the compliance and the capability of existing products to the user and stakeholder needs and to define a suitable strategy that could fulfil those needs.
The authors consider the study a definition of a spatially explicit tool for vulnerability assessment in deltaic systems and a contribution in terms of knowledge to the local management and planning component. This approach is a starting point for researchers and institutions to develop new down streaming services and evolve more detailed and scalable models to policies, with the aim of enhancing Earth monitoring and the detection of changes in coastal areas.
Supplementary Materials: The following are available online at http://www.mdpi.com/2073-4441/12/10/2830/s1, Table S1. The main characteristics of the satellites and band sets used. Table S2, Acquisition dates of the Landsat images.; Figure S1. Bayesian network developed in GeNie 2.4 to model the vulnerability assessment of the Po Delta (Italy) for reference scenario (a) and future SLR scenarios (b); Figure S2. Tornado diagrams for highly pessimistic scenario in the central macro-sector.; Table S3. Strength of influence.