A Bayesian Network methodology for coastal hazard assessments on a regional scale: The BN-CRAF

Hazard assessment is one of the key elements to be included in any coastal risk assessment framework. Characterizing storm-induced erosion and inundation involves the assessment of the coastal response under the forcing of a stochastic source (the storm), acting on a variable morphology (the beach) and inducing some damages. Hazard assessment under any present or future scenario will be affected by uncertainties either associated to the models used, the definition of climate conditions, and the characterization of the coastal morphology. In this context, Bayesian Networks (BN) can effectively address the problem as they allow accounting for these uncertainties while characterizing stochastically the system response and giving insight on the dependencies among involved variables. In this work, a BN-based methodology for storm-induced hazard assessment at regional scale is presented. The methodology is able to account for uncertainties associated with included models and forcing conditions through Monte-Carlo simulations. It produces distributions of erosion and inundation hazards under given scenarios allowing conditioned hazard assessments as a function of storm and morphological variables. Results are compared to hazards evaluated using an existing Coastal Risk Assessment Framework (CRAF), at two locations of the Catalan coast already identified as hotspots for storm-induced erosion and/or flooding.


Introduction
As coastal regions continue to grow in population size and density, they are also increasingly exposed to coastal storm-induced hazards, such as inundation and erosion (IPCC, 2012;IPCC, 2013). The projections of rising sea levels, the existence of background coastal erosion and other concerns related to climate change, such as increases in the magnitude and/or frequency of storms (Lionello et al., 2008, Conte andLionello, 2013;IPCC, 2014) or changes in the directionality of incoming waves (Casas-Prat and Sierra, 2013), highlight the need for hazard assessments that consider these factors. Applications of these assessments include designing coastal management plans, which will require a specific chapter on coastal risks, as recognized in the protocol of Integrated Coastal Zone Management in the Mediterranean plan (UNEP/MAP/PAP, 2008). This may include the need to identify coastal hotspots at regional scale, to compare assessments for different scenarios, and the development of early warning systems (EWS) among other elements (e.g., Ciavola et al., 2011aCiavola et al., , 2011bVojinovic et al., 2014;Oumeraci et al., 2015;Van Dongeren et al., 2018). The problem to be solved with these approaches is complex due to its multidimensionality (both in terms of multiple and interdependent hazards and vulnerabilities) and their multiscale nature (both in space and time). Moreover, hazard estimation comes with a number of associated uncertainties that can become especially large when dealing with future projections (Vousdoukas et al., 2018).
A widely used approach in risk assessment is the Source-Pathway-Receptor (SPR) model (Sayers et al., 2002;Narayan et al., 2014). This is a conceptual model which describes how a given risk propagates across a given domain from the source to the receptors. In this particular case, it is applied to storm-induced hazards (erosion and inundation) assessment, and the problem is schematized in terms of a source (storms), the pathway (beach or coastal morphology) and receptors (elements of interest) at the coast. To this end and under the common scarcity of direct observations, storm-induced hazards are usually assessed by using predictive models which are fed with information on both the source and the pathway. However, the implementation of a suite of complex numerical models at a high resolution (e.g., 1-10 m spatially, 1 s-1 h temporally) at the regional scale (e.g., hundreds of km) may not be feasible or time/cost efficient. Furthermore, the probabilistic definition of the hazard component should preferably be performed based on the response (Garrity et al., 2006;Callaghan et al., 2008;Sanuy et al., 2020) due to the nonlinear and multidimensional dependencies involved in the driving processes (see e.g., Hawkes et al., 2002;Masina et al., 2015;Lin-Ye et al., 2016). This implies the simulation of the erosion and inundation induced by a large set of conditions characterizing the existing storm climate. Additionally, model error is an important source of uncertainty that should be included in the assessment. Simple parametric models which use only bulk information on the source (e.g. peak Hs, Tp, direction and duration) and the pathway (e.g., slope, grain size, or beach height and width) can be a suitable option to generate sufficient hazard estimations at a reduced computational expense (e.g., Stockdon et al., 2006;Mendoza and Jim� enez, 2006). The drawback is the simplification of the storm characteristics to a set of parameters and the simplification of the morphology to the bulk definition of a 1-D cross-shore profile. When a long coastal stretch (~1 km) is represented by the bulk characteristics of a single profile, the pathway is treated deterministically, as neither spatial (alongshore) nor temporal (seasonal) variability is included in the assessment.
All this makes that any robust methodology for hazard assessments needs both to account for the existing variability in source and pathway and to tackle model uncertainties. However, this has to be done in a reasonable cost-effective manner, and it is in this aspect where Bayesian Networks (BNs) have demonstrated their versatility and utility in efficiently combining multiple variables to predict system behaviour. They have been used e.g. in groundwater flow predictions supporting decision making (Fienen et al., 2013), habitat protection against natural hazards (Palmsten et al., 2013;Gieder et al., 2014), coastal vulnerability and shoreline evolution to sea level rise (Gutierrez et al., 2011;Plant et al., 2016), postevent hazard and damage assessment (Van Verseveld et al., 2015;Poelhekke et al., 2016) and disaster risk reduction assessments Sanuy et al., 2018). BNs can handle nonlinear systems by means of conditional probability tables, which represent the interdependencies of the variables involved, and are solved through data training. They are not computationally expensive, can work with data from different sources (e.g., modelled, observed, or even opinions) and have a simple and intuitive graphical structure that is easily understood by nontechnical users (Uusitalo, 2007). Notably, they can be used to represent the SPR scheme through the dependency relations that physically, or even psychologically, exist between the different steps (e.g., Straub, 2005;J€ ager et al., 2018) and thus, they can be adapted to assess many kinds of natural hazards and impacts on many kinds of receptors.
The general aim of this work is to develop a methodology for the assessment of storm-induced erosion and inundation at regional scale with the purpose of hotspot identification and characterization under present conditions and future scenarios. Within this context, the development and implementation of a BN-based hazard assessment methodology is presented. It is based on the Coastal Risk Assessment Framework (CRAF) phase 1 methodology , which has been validated and applied across different study sites within the RISC-KIT project Ferreira et al., 2018), thereby indicating its robustness. The CRAF phase 1 combines a probabilistic treatment of the source (storms) with a deterministic treatment of the pathway (coastal morphology), without including model error uncertainties. The here presented new methodology is able to deal with the small-scale variability of coastal morphology (20-30 profiles per kilometre) and with model uncertainties by using Monte-Carlo simulations based on known model errors. The method is applied under current conditions and given scenarios defined in terms of SLR and observed rates of background erosion, where Monte-Carlo simulations are used to incorporate their associated variability. The proposed methodology, while maintaining a relatively simple structure, represents an advantage with respect to other existing approaches where the treatment of the source is deterministic or do not perform the statistics focusing on the hazard variables (e.g., use of hurricane categories as levels, Stockdon et al., 2007;use of event approach, e.g., Villatoro et al., 2014;Armaroli and Duo, 2018; use of a homogenous representation of the source in Poelhekke et al., 2016). Moreover, by accounting for the spatial variability of the coastal morphology, it also represents an improvement with respect to existing methodologies adopting a deterministic treatment of the pathway (Callaghan et al., 2008;Bosom and Jim� enez, 2011;Ballesteros et al., 2018. Notably, none of the referenced works included model errors in the analysis. This methodology is applied at two sectors of the Catalan coast (NW Mediterranean), and obtained results are compared to those of the semideterministic CRAF phase 1 . One sector is a 2 km long fully rigidized coastline composed by a rip-rap revetment protecting a coastal railway. The other one is a 3 km long sandy coastal fringe protecting a low lying deltaic area with campsites in the hinterland. Both locations have been identified as coastal hotspots to storm impacts in Jim� enez et al. (2018).
The structure of the paper is as follows: (i) the second section describes the data used, (ii) the third section presents the BN-based proposed methodology, (iii) the fourth section shows and example of the application of the method, followed by (iv) the discussion and comparison with CRAF phase 1 and (v) the main conclusions.

Data
To perform this analysis, data on waves and water levels to characterize the forcing, and on coastal morphology to characterize the pathway, are required. The assessment of coastal inundation and erosion hazards requires a long-term series of wave conditions and water levels that have the appropriate spatial and temporal coverage. This work uses hindcast waves from the Downscaled Ocean Waves dataset (Camus et al., 2013) derived from the Global Ocean Waves (Reguero et al., 2012) and hindcast surges from the Global Ocean Surge dataset (Cid et al., 2014), which were obtained at a node located in front of the Tordera Delta at ~20 m depth and cover the period 1950-2014. Simultaneous mean water levels were available at the same resolution of the GOS dataset.
The coastal morphology has been represented by using LIDARderived topography from the Institut Cartogr� afic i Geol� ogic de Catalunya. The data were provided as high-resolution digital elevation models (DEMs) with 1m � 1m grid cells and a vertical precision of 2-3 cm (Ruiz et al., 2009). This was complemented with data obtained from a topo-bathymetric survey provided by the Ministry of Agriculture, Fish, Food and Environment. Data up to 1.5 m water depth were obtained by total station topographic surveying and, data down to 50 m water depth were recorded with a dual-frequency echosounder with a nominal accuracy of 7 cm.
The long-term shoreline evolution is obtained from a 25-year shoreline dataset (Jim� enez et al., 2019) and the global SLR projections correspond to the RCP8.5 IPCC (2014).

The BN-CRAF methodology
In the following, the hazard assessment general scheme is introduced while presenting the different steps involved in the BN-based methodology and comparing with the corresponding steps of the previously developed CRAF phase 1 methodology. Then, methods followed in each step are described in detail.

The methodological scheme
The general methodological scheme follows CRAF phase 1  which was applied to the Catalan coast in Jim� enez et al. (2018).
In CRAF phase 1 (Fig. 1A), the source is characterized by identifying all the storms in a long (e.g. 60 years) time series of waves and water levels. The pathway (morphology) is characterized by a single cross-shore profile for each ~1 km coastal sector. This sector-characteristic profile can be calculated by simple beach profile averaging, or by choosing the worst-case profile (e.g. lowest dune and narrowest beach). The second option was applied in Jim� enez et al. (2018). Then, parametric models are used to estimate hazard magnitudes (e.g. total water level at the coast for inundation, and shoreline retreat for erosion) for each storm. Obtained hazard values are fitted by an extreme probability distribution, and finally, erosion and inundation magnitudes associated with selected probabilities/return periods are transformed to hazard indicators (see e.g. Ferreira et al., 2017) within a scale ranging from 0 to 5, being 0 no hazard and 5 the most hazardous situation. These indicators combine information of the hazard magnitude with basic information on the relative exposure of the hinterland to each hazard (e.g. beach width for erosion and, dune/berm height for inundation). The final output is, thus, a single intensity value per sector for each hazard associated with a given return period.
In the here presented BN-based methodology (Fig. 1B), the source is characterized as in CRAF phase 1, but the pathway is characterized by using a set of profiles (~10-30, Fig. 2) to properly capture the existing morphological variability within each sector. Hazards magnitude is computed by means of parametric models for each storm and profile. To account for the uncertainty associated with model errors, Monte-Carlo simulations (Hastings, 1970) are used. Thus, for each profile and storm, multiple erosion and inundation hazard values are calculated, which are later transformed to hazard indicators using the same scale as in CRAF phase 1. These hazard estimations together their corresponding source (storm characteristics) and pathway (profile characteristics) data are used to train the BN. The final output are probabilistic distributions of the different erosion and inundation hazard indexes at each sector, which consider model uncertainties and preserve the information on conditional dependencies with source and pathway.
In addition to detecting hotspots along the coast, the BN methodology is used to assess hazard values under given scenarios associated with future global SLR and local background erosion projections, where the Monte-Carlo approach is applied to simulate multiple conditions per scenario.

Source and pathway characterization
For comparison purposes, the areas to be studied are divided into ~1 km sectors, with the same spatial division as in the CRAF-phase1 application . Each sector is represented by a number of profiles (Fig. 2). Profile selection is aimed to capture the existing morphological variability in each sector, in terms of slope, beach height (i.e. dune or upper berm or embankment height), and beach width (Fig. 3b). The shape of the first part of the hinterland is also taken into account for profile selection. Notably, rigidized coastal sectors such as Cabrera South (CS) and Cabrera North (CN) require a lower number of profiles (16 and 9 respectively) than natural beach sectors with larger variability such as Malgrat North (MN) and S'Abanell (SA) with 20 and 32 profiles respectively. Malgrat South (MS) is more homogeneous and was characterized with 15 profiles.
Coastal storms have been identified in the dataset (see location of the node in Fig. 2) by means of a double threshold P.O.T analysis as similarly conducted in Sanuy et al. (2020). The 0.98 quantile (Hs ¼ 2 m) was used as the first threshold to characterize storms by controlling the event duration. The 0.995 quantile (Hs ¼ 2.6 m) was then applied to the subset to identify the extreme events. The resulting dataset is composed of 179 storms (~3 storms per year), in which the complete hourly evolution of the different wave conditions is stored, i.e., significant wave height (Hs), peak period (Tp), storm surge and wave direction. The storm duration is calculated as the time during which Hs exceeds the selected threshold. Fig. 3a shows scatterplots of the main characteristics of the dataset.

Hazard indicators
In this work, the hazards have been estimated similarly to that of Bosom and Jim� enez (2011) and applied similarly to Jim� enez et al. (2018). The inundation potential is calculated using, as proxy, the total water level at the coast (TWL), characterized by the wave-induced run-up and the surge level. The run-up models used are the Stockdon et al. (2006) model [eq (1)] for sectors composed of sandy beaches and the EurOtop (Allsop et al., 2007) model [eq (2)] for rigidized sectors with protected revetments. These models use significant wave height in deepwaters (Hs) or breaking conditions (Hb), wave length in deepwaters (Lo) and beach slope (tanβ, calculated from the 1 to 1 m elevations) to estimate the vertical level potentially reached by waves at the coast. The EurOtop (Allsop et al., 2007) model also includes factors accounting for wave directionality (γ β ) and surface friction (γ f ).
Ru ¼ 1:1 0:35 tan βðHs LoÞ 1=2 þ ½Hs Loð0:563 tanβ 2 þ 0:004Þ� The erosion hazard has been estimated by using the Mendoza and Jim� enez (2006) model, which has been widely applied to estimate storm-induced erosion potential (see e.g., Armaroli and Duo, 2018;Silva, 2019). The model is given by the following: ( 3) where ΔV is the eroded volume at the beach face (Fig. 4), Tp is the wave period at the peak and w f is the sediment fall velocity.
where ΔX is the beach retreat calculated as a function of the eroded volume, the beach height (ZB) and the pivot point depth (d*), (Fig. 4).
Finally, hazard estimations are converted into a 0-5 hazard intensity scale, as in Jim� enez et al. (2018). For the inundation hazard, the TWL is compared against the profile to obtain the potential reach of the flooding by means of the bathtub approach. This magnitude is then compared with the beach width (BW) to derive the hazard index (Table 1). The erosion index, in a different manner, is derived by comparing the resulting post-event BW (i.e., the pre-event BW subtracting the storm-induced retreat) with a reference width. Following Jim� enez et al. (2018), the reference width of the different levels is based on the average beach retreat for the 10-year return period at the study site ΔX 10 , with a value of 7.5 m (Table 1).

Future scenarios
The assessment of the possible evolution of hazards requires having information on future morphology (pathway) and forcing (source) conditions. In the present application, this is done by estimating how such elements will vary under a given scenario by a given time horizon.
As an example, to assess how hazard will evolve at mid-term scale (from years to few decades), we have considered the current decadal-scale coastline evolution to build the future morphology (beach width) at selected time horizons while using Mendoza and Jim� enez (2006) to assess the corresponding-storm induced retreats. When the analysis is extended at the long-term scale (several decades to centuries), we have to consider potential changes in forcing conditions. If storminess is not expected to change significantly (see e.g., Somot et al., 2019 for the Mediterranean basin or Bender et al., 2010 for Atlantic hurricanes), the main variable affecting analysed hazards is the long-term change in sea level. In such a case, this is done by considering the future mean sea level under a given SLR scenario by the selected time horizons. In the case that long-term changes in storminess are expected, storm properties must be modified accordingly.
In this study, mid-term background erosion along the sandy beach sectors has been estimated by analysing shoreline position changes from aerial photographs available at different time frequencies that cover the last 25-30 years. The estimated average retreat for the 3 sandy beach To introduce SLR-induced effects, the IPCC (IPCC, 2014) RCP8.5 scenario was considered in this study, which accounts to 0.30 m and 0.63 m of global mean SLR by the years 2050 and 2100, respectively, with respect to the period 1986-2005.

Monte-Carlo approach
All hazard models and future projections are affected by errors and uncertainties. The available knowledge on these errors, or reasonable assumptions on the uncertainties, may be included in this hazard assessment. Known model root mean squared error (RMSE), and confidence ranges (Table 2) have been included in hazard estimations by means of Monte-Carlo simulations (Hastings, 1970). Other existing uncertainties may also be included following a similar procedure.
Uncertainty in the run-up has been estimated as follows. For the Stockdon et al. (2006) model, we use the RSME of the formula as an estimator of the standard deviation of the variable . It is assumed that the variable follows a normal distribution with the mean calculated by [eq (1)] and a standard deviation equal to the model RSME (0.32 m). For the EurOtop (Allsop et al., 2007) model, the knowledge of the standard deviation of the coefficient in [eq (2)] ( Table 2) is included in the assessment along with an estimation of the uncertainty of the actual value of the revetment friction coefficient.
Uncertainty of erosion potential (ΔV) has been similarly estimated. First, the RMSE of the eroded volume given by the Mendoza and Jim� enez (2006) model is used. Then, the uncertainty on the real height of the eroded profile (i.e., distance from submerged pivot point to beach top) is included by assessing the variability of the pivot point position assessed from simulations with XBeach 1D under different storm intensities (d*, Table 2).
Confidence ranges on projected conditions under future scenarios are included in the assessment. In this case, the 95% confidence interval of the SLR estimates (IPCC, 2014) is used to derive the standard deviation of the variable, assuming a normal distribution. The standard deviation of shoreline evolution rates within each sector (Jim� enez et al., 2019) was used to simulate multiple values of beach width per future scenario (Table 2).

Table 1
The hazard index levels as a function of the beach width (BW) against inundation reach or beach retreat. All units are in meters, and ΔX 10 stands for average profile retreat associated with TR ¼ 10 years for the site, which has a value of 7.5 m.

Table 2
Ranges of the uncertain variables included in the assessment. Note that total number of cases is always 100 for hazard variables and 20 for future scenario variables The total number of cases can be the result of the simulation of a single variable or of the combination of simulations from 2 variables. Therefore, following the Monte-Carlo example, normally distributed values of the uncertain variables for each combination of storm and profile characteristics are simulated. This extends uniformly to the dataset maintaining the natural variability of the storm climate (source) and the morphology (pathway). The total number of cases introduced in the BN training is 100 values of erosion and inundation hazards per profile and storm combined with 20 values of MSL or beach width per future scenario.
In this application, Oi represents the hazard-estimators, such as inundation reach or beach retreat, and Pi represents the parent conditions for such hazard's results. Basically, Pi is a set of variables containing the storm and scenario morphological characteristics, e.g., Hs, duration, direction, surge, period, slope, beach height and beach width. By this method, the hazard outputs at each profile are mapped through the conditional probability tables to their parent inputs in the BN, and the sector results are obtained by weighting each profile according to their contribution length over the sector. Source variables, which are known to present interdependencies can also be interconnected in the BN, e.g., establishing Hs, Tp, direction and duration as parents of storm surge. Fig. 5 shows the BN structure used in this work. Arrows denote dependency relations between the variables. The forcing and morphologic variables are parents of the preliminary hazard output (i.e., Ru, and eroded volume). Ru, eroded volume, and morphology are parents of the main hazard outputs (i.e., inundation reach and beach retreat). Finally, these hazards are combined with beach width to obtain the final hazard index.
As previously mentioned, the uncertainties on a number of variables are included in the BN by means of Monte-Carlo simulations. Fig. 6 shows a graphical example of how this information goes into the BN and propagates to the output. Simulated variables are assumed to be Gaussian. For each profile, storm and scenario, the BN hazard output is a discretized probability distribution accounting for the considered uncertainties (see the cases in Table 2). The total number of cases is a function of the number of profiles, number of storms, number of uncertain variables considered, and number of simulations per uncertain variable. As an example, a sector with 20 profiles and 150 storms would have 100 Monte-Carlo outputs of TWL for each specific profile, storm and MSL condition. Since each future SLR scenario is defined with 20 Monte-Carlo values of MSL, this leads to a total number of 6*10 6 cases per time projection to feed the BN.

Current state hazard profile
The first direct result produced by the BN is the unconstrained distribution of all the variables. Data training fills the conditional probability tables with all simulated cases, being each case a set of the source (waves and water level parameters), pathway (beach characteristics) and hazard variables (Fig. 5). When no conditioning is applied to the BN, it shows the prior probabilities, i.e. the marginal distributions of each variable without information on the other variables ([eq (5)]). The distributions of the calculated hazard indexes (Fig. 7-B) represent the probabilities of each hazard level at each sector, both considering the variability of the source (storms) and the variability of the sector morphology. This consideration allows a more representative and detailed sector intercomparison than from the semideterministic approaches ( Fig. 7-A), which aims for a single hazard value for the sector associated with a given return period (e.g., Jim� enez et al., 2018). Note that results are presented here with the same sector definition (~1 km) as in Jim� enez et al. (2018) for comparison purposes, but the framework allows integration at any scale.
The percentages in the hazard profile results can be interpreted in two ways. For example, looking at the erosion distribution at SA as follows: (i) from the perspective of an incoming event of unknown characteristics, there is a 6% probability of hazard index 5 occurring at some location within the sector, or (ii) from the perspective of an average incoming event, 6% of the total length of the sector is estimated to be affected by an intensity 5 erosion hazard. This dual interpretation is a consequence of the combination in the BN of the storm-climate and the morphological variabilities which allows comparing differences between sectors under the same storm-climate with similar representative profiles but with different morphological variability, which in the case of a deterministic treatment of the pathway would produce the same result.
Focusing on the inundation hazard in Tordera Delta, both MN and MS sectors present percentages >5% (a typical threshold for statistical significance) of index 5, although they show different distributions. In SA, where the topography is higher, the maximum inundations are less frequent but intermediate hazard indexes present higher probabilities, which is a consequence of its narrow beach. In MN the response is the opposite due to its morphology, characterized by wider beaches and a lower hinterland. Since the inundation reach is approximated by using the bathtub approach, when water exceeds the berm height a large extension is potentially inundated, and therefore, the hazard index jumps from low (0-1) to extreme values (5). This dual state behaviour is reinforced in MS as a consequence of its even wider and more homogeneous beaches. Identifying the different hazard responses is of key importance if a hazard profile is to be combined with vulnerability data to obtain a risk profile. In Cabrera de Mar, both sectors have similar inundation hazard distributions, while the semideterministic approach (Fig. 7-A) showed a different index between the sectors because of its representation by single-transect and no inclusion of uncertainties. Notably, using a single transect fixes a single beach height and hinterland shape, and the inundation reach upon the profile chosen for CS for the TR ¼ 100 yr was higher than 40 m and lower than 60 m. However, an even higher profile but with lower hinterland topography may exist, which combined with the inclusion of model errors in the analysis (higher run-ups) may result in the possibility of inundation reaches higher than 60 m leading to an intensity 5 inundation hazard. A similar effect can be observed comparing CRAF phase 1 results with the BNbased distributions for inundation in MS.
Moving to the erosion result at Tordera Delta ( Fig. 7-B), it can be observed how the hazard intensity decreases from north (SA) to south (MS). In addition, the use of multiple profiles and the inclusion of model errors shows that erosion intensity 5 is actually probable at SA (6%), whereas it was not shown in the semideterministic approach (Fig. 7-A). In addition, SA and MN show a quite different profile, whereas in Fig. 7-A they appear as equivalent sectors. SA shows a generalized vulnerability to erosion, having significant probabilities throughout all intensity levels. Meanwhile, MN has a large probability of no-hazard at 75%, and the remaining probabilities are mainly at medium intensity levels, with some residual probability of occurrence of a top-level erosion episode. This outcome is a consequence of SA having narrower and more uniform beaches along the costal stretch while MN has various areas (wide and narrow) that behave differently, which are also detectable with the BN (see results on backward propagation and hazard source below).

Identification of hazard sources
One of the main advantages of using BN is its capacity to assess correlations between variables by means of the information stored in the conditional probability tables. Thus, the BN can also be used to answer questions about the sector behaviour because conditional probability tables are built preserving the relationships between variables connected in Fig. 5. By means of backpropagation, i.e. imposing values in output variables (hazards) and assessing the distributions obtained in their parent variables (source and pathway), the BN can assess storm or morphological characteristics associated with induced hazard indexes of interest.
To illustrate the potential of BNs in this aspect, the case of the erosion hazard at MN is presented in Fig. 8, where the number at the left side of the boxes represent the ranges in which each variable is discretized, and numbers next to black bars represent the discretized probability density at each bin. This case shows a probability of hazard at � 4 of 6%, with a residual of 3% associated to an index 5. The prior distributions represent unconstrained probabilities of some variables of the sector. As an example, the beach width distribution shows a sector mainly with a 40-100 m wide beach, with some extension (~22%) of a relatively narrower coast.
The BN can be constrained to show how the source and the morphological variables change, e.g. when constraining the maximum erosion hazard level (i.e., 4 and 5), obtaining the conditional probabilities related to that specific case (posterior distributions) and identifying the combination of variables that result in such hazard values. As observed, sectors that have a beach width between 10 and 40 m (representing only ~22% of the sector's locations) concentrate the ~86% of the probability of getting a 4 to 5 erosion intensity.
Focusing on the hydrodynamic variables, it can be observed how the duration of the event is the main driver of high intensity erosion episodes at the sector (Fig. 8, posterior). Higher Hs and Tp also favor a large erosion but the changes between prior and posterior are less intense.
Since the erosion hazard is quite sensitive to the beach width in the analysed sector, the BN can be used to statistically assess the expected erosion hazard for existing width values along the area. Thus, Fig. 9 shows the mean erosion hazard index and its 95% confidence range, together the probability of occurrence of an intensity 5 by constraining beach widths to a range between 15 and 70 m. As it can be observed, a minimum width of 40 m is needed in MN to fully avoid an erosion index 4 to 5 (probability of index 5 < 0.05 and confidence interval out of intensities 4 to 5). With a value of BW of 70 m or more, all index values become zero, which is also important information for other uses of the coast, such as the recreational use of the beach. Note that the results presented in Fig. 9 are for current conditions and do not include information on the temporal evolution of the study site.

Hazard assessment at future scenarios
When the scenario-datasets are introduced in the BN, the hazard profile at different time horizons is obtained. Semideterministic indexes (Fig. 7-A) in future scenarios will only show the changes in the sector hazard's single-level, and once it reaches its maximum (index 5), it will stop showing changes unless new levels are created. Alternatively, the BN-approach allows the evaluation of changes in the probabilities of the different hazard intensities at different time horizons. To illustrate this, two main results are presented: erosion hazard mid-term horizons due to background erosion in the Tordera Delta assuming that past shoreline trends remain in the future (Fig. 10), and long-term horizons of the inundation hazard due to SLR at Cabrera de Mar (Fig. 11).
In the Tordera Delta, SA is the sector showing the current highest index of the erosion hazard, but MN shows a faster progression towards intense erosion events in future scenarios. Notably, in 5 years, SA still shows a general situation with larger frequencies of medium-high erosion episodes, but MN reaches a similar frequency for hazard index 5. At the 10-year horizon, MN is clearly in a worse situation than SA for index 5 frequencies but still better for probabilities of 0-4 erosion intensities. At the 20-year horizon, MN is the most vulnerable sector to erosion with almost all events causing at least an index 3 situation, with 86% frequencies of index 5. In contrast, MS remains unaltered across scenarios showing a no-impact profile. (Fig. 10).
In Cabrera de Mar, both CN and CS show similar responses to SLR at the 2050 and the 2100 RCP 8.5 scenarios. CN shows a higher increase in frequency of medium-intensity events whereas CS keeps on showing slightly larger probabilities of index 5 situations. In this case, this distinction is crucial, since index 3 intensities already reach the infrastructure behind the revetment also causing significant disruptions which makes the northern sector more vulnerable (in general terms) than the southern one (Fig. 11).

Discussion
In this work, a BN-based methodology for regional storm induced hazard assessment has been presented and applied at two sectors on the Catalan coast to illustrate its potential.
Obtained results highlight the benefits of estimating a hazard probability distribution for each sector. The explicit inclusion of uncertainties (i.e. including model and parameter errors when producing results) and the extensive coverage of the morphological characteristics showed the potential highest-intensity hazards not detected in CRAF phase 1 (e.g., erosion at SA or inundation at CS, Fig. 7). In other words, the methodology has been successful in hotspot identification, as in , but allowing detailed sector intercomparison and reducing hazard underprediction due to the omission of uncertainties. In order to illustrate this, Fig. 12 shows an additional comparison between results obtained here and those using CRAF phase 1. The BN-method outputs the complete curve of Ru vs the hazard index. This is obtained by constraining the BN at different Ru levels and assessing the distribution of the hazard indicator. On the other hand, the application of CRAF phase 1 results in a single hazard index value associated to a given return period (illustrated by a point in the graph). At MN, both results are equivalent (see also Fig. 7) and the BN output allocates the Ru associated to TR ¼ 100 yr to a hazard index 5 with a ~95% of probability. At CS, the CRAF phase 1 output is a hazard index 4 for TR ¼ 100 yr (which has a 0% probability of reaching index 5 for the associated Ru, which corresponds to 5.6m), while the BN-method outputs some probability of having index 5. This difference between both approaches is due to the configuration of the area, which is controlled by the embankment height and the topography of the hinterland. The consequence is the existence of an abrupt increase in the expected mean hazard index and the probability of occurrence of hazard intensity 5 (and associated 95% confidence ranges) for Ru values higher than 5.5 m. The inclusion of model uncertainty in the BN assessment produces values for Ru up to 6.5 m which will fall into the index 5 category. This explains the difference in hazard intensities obtained using both methods showed in Fig. 7. In both cases, it can be observed how the BN-method gives more detailed information on the relation between Ru values and expected hazard indexes (Fig. 12) than the CRAF phase 1 single value, showing explicitly the expected variability due to the morphology and model uncertainties.
Accounting for this variabilities and uncertainties is especially important when performing mid-term and long-term hazard assessments. When assessing hazards at future horizons, changes in the frequencies of low-intermediate hazard intensities can be as important as those of the extreme (large TR) (Figs. 10 and 11). Thus, obtained changes in the hazard distributions give complete information about the sector response. This allows a preliminary identification of tipping points for coastal adaptation, as information of intermediate hazard indexes at different time horizons is also needed for such purpose. This also allows identification of sectors with different responses over time but with nearly the same current hazard profile (Fig. 11). Mid-term and long-term scenarios are based on modelled projections with associated uncertainties that are often estimated with a given error or expected range (e.g., expected SLR for a given horizon and the associated 95% confidence range). This is included in the hazard assessment, and it should be considered for sector intercomparison since it prevents hazard underprediction that may arise from single-value future estimations.
The methodology can also produce results of sector morphology and forcing characteristics conditioned to a given hazard intensity (Figs. 8 and 9). Thus, the BN, through back-propagation, indicates the potential causes of given hazard index values, giving insights into which processes are important or even to what remediation measures could be effective. Moreover, the BN can produce hazard distributions conditioned to specific storm conditions (i.e., known values of Hs, Tp, duration, direction and surge) in real time. Thus, since these variables are produced by operational hydrodynamic forecasting systems, the tool can perform as a preliminary regional EWS without the need of time-consuming, detailed, morphodynamic, process-based modeling.
In this work, we have applied the developed methodology by integrating hazards at ~1 km spatial units. This was done mainly for comparison purposes with the results obtained in Jim� enez et al. (2018). However, the BN can integrate results at any spatial and scale depending on the physical-process or management questions to be addressed.
One limitation of the presented methodology is that a long record of storms is needed to produce robust results. As opposed to CRAF phase 1, no extreme analysis is done and therefore, no extrapolation of storm conditions is performed to consider events of higher magnitude than recorded ones. However, as highlighted in Fig. 12 and from the direct comparison of obtained results, for the analysed record length (60 years), this seems to have a smaller impact than not accounting for model errors and morphological variability in the sense that CRAF phase 1 results associated to a TR ¼ 100 yr underestimate hazard indexes when compared to the BN-CRAF using 60 years of storm events.
In the current application, some assumptions have been imposed for the sake of simplicity and test the output. However, they are independent of the method, and they could easily be substituted by other choices. Thus, future scenarios due to background shoreline retreat have been built assuming that measured past trends will keep unchanged in the future. However, they could be substituted by time-varying rates or results from other models. Moreover, process-based models (e.g., 1D-Xbeach, Roelvink et al., 2009;SBeach, Larson and Kraus, 1989;LIS-FLOOD, Bates and De Roo, 2000;Bates et al., 2005) could be used to quantify erosion and inundation hazards instead of simple parametric models. The uncertainty associated with numerical modeling can be included through an ensemble of simulation outputs from different parametrizations or different models, similar to model ensembles used to transfer the uncertainties in future projections of SLR to the inundation maps (Purvis et al., 2008). Following this concept, many different sources of uncertainty (e.g., seasonality of the coast morphology), which were omitted in the present application could be included, with the only consequence being an increase in the size of the case-dataset feeding the BN. The current application was executed on a regular laptop, meaning that there is still room, computationally speaking, to explore new elements to include in the analysis, or increase the number of simulations (e.g. >20 for future projections), increasing the confidence of the results. In this sense, the approach could also be extended with socioeconomic data in order to translate hazard profiles into impact profiles.

Conclusions
A Bayesian Network-based methodology for storm-induced hazards assessment at regional scale has been developed. It has successfully identified coastal erosion and inundation hotspots in two sectors of the Catalan coast (NW Mediterranean) similarly to what was previously done with the CRAF-phase 1 method while giving further information.
The method accounts for the contribution of different sources of variability (i.e., variability in storms and coastal morphology), as well as uncertainties associated with the hazard models (i.e., model errors, parameter uncertainties and confidence ranges on future scenarios). As a result, instead of resulting in a single value, the Bayesian Network-based method provides hazard distributions which will indicate the probability of occurrence of hazards at any intensity.
The method can be easily used to forecast changes in storm-induced hazard levels by applying it under given future scenarios, based on changes in morphology and/or storm properties. In such a case, the method will predict the expected change in the probability of occurrence for different ranges of hazard indexes.
The inherent capability of Bayesian Network's to assess the interdependence between variables is used to identify source (storms) and pathway (beach morphology) characteristics responsible for given hazard distributions. When this is applied under different scenarios permit to identify harmful combinations taking place at specific locations and at given time horizons. From the management standpoint this will provide essential information for making decisions on selecting coastal adaptation alternatives.
The application of the method at two known hotspots of the Catalan coast highlighted concerning scenarios of erosion (Tordera Delta) and inundation (Cabrera de Mar) in 10-20 years from baseline, while allowing a better comparison and characterization of the ~1 km sectors within hotspots.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.