A Spatial Structure Variable Approach to Characterize Storm Events for Coastal Flood Hazard Assessment

: Over the last decades, the evaluation of hazards and risks associated with coastal ﬂooding has become increasingly more important in order to protect population and assets. The general purpose of this research was to assess reliable coastal ﬂooding hazard maps due to overﬂow and wave overtopping. This paper addresses the problem of deﬁning credible joint statistics of signiﬁcant wave heights Hs and water levels ζ , focusing on the selection of the sample pair that characterizes each sea storm, to evaluate the occurrence probability of extreme events. The pair is selected maximizing a spatial structure variable, i.e., a linear combination of Hs and Relaix, F., Zammit, P.S. Satellite cells are essential for skeletal muscle, speciﬁc to each point of the area at risk. The structure variable is deﬁned by the sensitivity of the ﬂooding process to Hs and ζ , as found by analyzing a set of inundation maps produced through a Simpliﬁed Shallow-Water numerical model (SSW). The proposed methodology is applied to a coastal stretch in the Venetian littoral (Italy), by means of a 30 year-long time series recorded at the “Acqua Alta” oceanographic research tower, located in the Northern Adriatic Sea in front of the Venetian lagoon. The critical combination of Hs and ζ forming the structure variable is presented in a map, and it can be related to the topography and the presence of mitigation measures. The return period associated with the two recent large storms that occurred in this area in 2018 and 2019 is also investigated. The proposed procedure gives credible occurrence probabilities for these events, whereas other approaches would consider them extremely unlikely.


Introduction
Coastal flooding is frequently associated with two mechanisms occurring during a sea storm: overflow and wave overtopping. In the first case, the freeboard of the local levee/defense structure is steadily below the water level, whereas in the second case the flood volumes depend on the oscillating nature of the waves. They are often threatened as separate mechanisms although waves and water levels clearly affect each other. The local water level in front of the coast is a combination of offshore sea level and wave setup, influenced by the breaking process. Similarly, the wave characteristics are affected by the water level during propagation, and hence these phenomena are interrelated.
The main motivation of this paper is related to the challenges posed by the appropriate definition of sample pairs to be used for bivariate (significant wave heights H s − water levels ζ) statistical analysis, in the framework of coastal flooding hazard assessment.
Bivariate statistical methods are significantly affected by the approach used for the selection of just two lumped values (H s , ζ) characterizing each stormy event.
Several authors (Coles and Tawn [1], Zachary et al. [2], Masina et al. [3], Li et al. [4], Hawkes et al. [5]) analyzed engineering problems considering a sample formed by the load maxima within the storm and assuming that these maxima are simultaneous. This procedure overlooks the correlation between wave heights and water levels within the event and systematically considers the worst possible scenario in terms of coastal flooding. Zachary et al. [2] also examined the approach where one variable is considered of primary between the incident significant wave heights Hs and the offshore sea levels ζ is investigated (for example in [2,3,12,13]).
One of the key aspects to define the joint statistics that arises in the discussion of Coles and Tawn [1] is the importance of defining a sample formed by independent extreme events. Sampling and dependence are closely linked (Mazas and Hamm [14]), therefore, (i) a declustering scheme is necessary to ensure independence and (ii) a criterion for the selection of the parameters that describe the drivers' combination within the cluster (event) is crucial to correctly account for the correlation structure.

Declustering Scheme for Sea Storm Identification
A sea storm is defined as a meteorologically induced disturbance to the local maritime conditions (i.e., waves and water levels) that has the potential to significantly alter the underlying morphology and expose the backshore to waves, currents and/or inundation (Harley [15]).
In order to specify the beginning and the duration of a sea storm, a declustering scheme can be used. For statistical convenience, the declustering procedure should identify independent sea storms. The identification can in principle be assessed using different variables and different thresholds. Frequently, a meteorological dataset based on wind records is either not available or not easily manageable and the identification involves the analysis of significant wave height H s and sea-level ζ.
For the purposes of this study, it is considered that sea storms are necessarily associated with wave storms, and therefore the proposed declustering scheme is only based on the values of H s . A wave storm is defined as a sequence of wave states during which H s is above a given threshold, named h CRIT . Usually, the threshold h CRIT is related to the average significant wave height < H s > calculated from its time-series, so it depends on the characteristics of the recorded sea states (e.g., h CRIT = 1.5 < H s >, Boccotti, [16]). In order to preserve independency of the data, some technicalities are also considered. It is still the same storm if H s falls below the threshold but only for a short time interval ∆t CRIT (Boccotti, [16], Mendoza et al. [17]), as shown in the central panel of Figure 1. For the choice of ∆t CRIT , it is suggested to analyze the autocorrelation function of H s and find the microscale λ (Taylor [18]) (λ is defined as the intercept of the parabola that matches the curvature of the autocorrelation at the origin with the time axis). For simplicity, it was assumed ∆t CRIT = λ. Figure 1 shows an example of wind speed (upper panel), significant wave height (central panel) and sea level (lower panel) measurements. The grey boxes highlight two independent wave storms, defined by the exceedance of h CRIT , that clearly point out the occurrence of sea storms in the wind time series. The sea level ζ is broken down into two sub-processes: the astronomical tide (ζ A ) and the residual (ζ R ). The latter (ζ R ) represents the sum of different contributions (Pasquali et al., [19]), such as storm surge induced by wind and pressure, basin seiching (Bajo et al. [20]), etc. For instance, in the illustrative figure, seiches are evident in the oscillations that appear in the ζ R series after day 9 and their amplitude fades after 4 days.

Sampling Procedure Within Sea Storms
During a sea storm, different combinations of simultaneous ζ and Hs can cause ing.
The typical approach would be to select a "a priori" combination such as:  Maximum Hs and simultaneous value of ζ;  Maximum ζ and simultaneous value of Hs;  Maxima of Hs and ζ during the storm regardless of their concomitance.
These approaches overlook the correlation between the variables or consider o be of primary importance. A method to overcome this sampling criticism, preservin correlation between ζ and Hs, consists of defining a third variable function of these dr Similar to Mazas and Hamm [Error! Reference source not found.], instead of consid an output variable (i.e., the overtopping discharge), a function that weights the eff the variables on the flooding can be used for the sampling: r = (ζ − < ζ >) + a (Hs − < Hs >) The selection of the coefficient a is based on the following considerations.

Sampling Procedure within Sea Storms
During a sea storm, different combinations of simultaneous ζ and H s can cause flooding. The typical approach would be to select a "a priori" combination such as:

•
Maximum H s and simultaneous value of ζ; • Maximum ζ and simultaneous value of H s ; • Maxima of H s and ζ during the storm regardless of their concomitance.
These approaches overlook the correlation between the variables or consider one to be of primary importance. A method to overcome this sampling criticism, preserving the correlation between ζ and H s , consists of defining a third variable function of these drivers. Similar to Mazas and Hamm [2], instead of considering an output variable (i.e., the overtopping discharge), a function that weights the effect of the variables on the flooding can be used for the sampling: The selection of the coefficient a is based on the following considerations.
Since the flooding is site-specific and depends on the topography, presence of defense structure, etc., the coefficient should not be uniquely defined, but related to the critical combination of H s and ζ that cause flooding.
For the evaluation of coastal flooding vulnerability maps, failure is defined in each grid cell of the map when the water depth h FLOOD exceeds a certain value h LIM , which is chosen in this study as 0.2 or 0.5 m, identifying two limit states.
The function that divides the safe region (not flooded) from the unsafe region (flooded) is hence computed through a "limit state" function g (H s , ζ) = 0 defined as: Water 2021, 13, 2556

of 14
The h FLOOD reached in each grid cell of the investigated area can be calculated with a numerical model for a set of pairs of H s and ζ, representative of the typical range for the site. For the present study, the h FLOOD values are computed with the Simplified Shallow Water (SSW) model developed by Favaretto et al. [21]. This raster-based inundation model solves a simplified form of the Shallow-Water Equations (suited for GPU acceleration) for each cell (pixel) of the domain. Several validations show that this model is capable of simulate wet/dry transitions (Favaretto et al. [22]) and real flood events (Favaretto et al. [21]). Figure 2 shows an example of the values assumed by h FLOOD (H s , ζ) for an illustrative pixel, located on the beach. In the right panel of the figure, dots are the h FLOOD (H s , ζ) reached for 90 SSW simulations carried out with different H s and ζ, the black line is the "limit state" considering a critical value of 0.2 m and the grey area is the unsafe region. grid cell of the map when the water depth hFLOOD exceeds a certain value hLIM, which is chosen in this study as 0.2 or 0.5 m, identifying two limit states.
The function that divides the safe region (not flooded) from the unsafe region (flooded) is hence computed through a "limit state" function g (Hs, ζ) = 0 defined as: g(X) = g(Hs, ζ) = hLIM − hFLOOD(Hs, ζ) = 0 (2) The hFLOOD reached in each grid cell of the investigated area can be calculated with a numerical model for a set of pairs of Hs and ζ, representative of the typical range for the site. For the present study, the hFLOOD values are computed with the Simplified Shallow Water (SSW) model developed by Favaretto Figure 2 shows an example of the values assumed by hFLOOD (Hs, ζ) for an illustrative pixel, located on the beach. In the right panel of the figure, dots are the hFLOOD (Hs, ζ) reached for 90 SSW simulations carried out with different Hs and ζ, the black line is the "limit state" considering a critical value of 0.2 m and the grey area is the unsafe region. The coefficient a is reasonably the sensitivity of the limit state to the variable Hs relatively to ζ, and it is found based on the slope of the "limit state" function for each point of the vulnerability map. Actually, the slope should be evaluated at the point that is most likely to induce flooding. This point can be defined exactly once the joint statistics of Hs and ζ is completely known, and it can be assessed with subsequent iterations. In the first iteration, all data points are used without declustering, and the pairs [Hs, ζ] are considered independent. The proposed procedure follows the First Order Reliability Method (FORM The coefficient a is reasonably the sensitivity of the limit state to the variable H s relatively to ζ, and it is found based on the slope of the "limit state" function for each point of the vulnerability map. Actually, the slope should be evaluated at the point that is most likely to induce flooding. This point can be defined exactly once the joint statistics of H s and ζ is completely known, and it can be assessed with subsequent iterations. In the first iteration, all data points are used without declustering, and the pairs [H s , ζ] are considered independent. The proposed procedure follows the First Order Reliability Method (FORM [23]). The variables [H s , ζ] are transformed to equivalent independent standard normal random variables U = (U H , U ζ ), and the "limit state" function g(H s , ζ) is linearized (in similitude with the approach used in Martinelli et al. [24]). The left of Figure 3 shows the standard space of the "limit state" computed with the SSW model and linearized. The point that has the highest failure probability is the one with the shortest distance from the limit state to the origin in the standard space. At this point in the physical space, the slope of the "limit state" (a) is computed (Figure 3 Right). The value r LIM indicates the limit of the structure variable r over which failure occurs. Reference source not found.]). The left of Figure 3 shows the standard space of the "limit state" computed with the SSW model and linearized. The point that has the highest failure probability is the one with the shortest distance from the limit state to the origin in the standard space. At this point in the physical space, the slope of the "limit state" (a) is computed (Figure 3 Right). The value rLIM indicates the limit of the structure variable r over which failure occurs. The red segment highlights the point with the shortest distance from the origin. Physical space (right): limit state function derived from back-transformation in the standard space. The dashed line highlights the slope of the "limit space" and the red line is the structure variable limit over which failure occurs. In both spaces also the assumed joint probability density function of UH and U and of Hs and ζ is drawn.
The spatial structure variable r is therefore the most likely combination of Hs and ζ that determines the maximum flooding in that particular pixel. During each sea storm, the declustering is carried out identifying the maximum r, as shown in Figures 4 and 5.  The dashed line highlights the slope of the "limit space" and the red line is the structure variable limit over which failure occurs. In both spaces also the assumed joint probability density function of U H and U ζ and of H s and ζ is drawn.
The spatial structure variable r is therefore the most likely combination of H s and ζ that determines the maximum flooding in that particular pixel. During each sea storm, the declustering is carried out identifying the maximum r, as shown in Figures 4 and 5.
The maximum r for each storm provides the critical pairs of H s and ζ for each location, i.e., for each value of the coefficient a. In theory, these pairs could be used iteratively to define more approximate joint statistics and assuming the actual correlation structure. However, we consider this step unnecessary.
The samples, r, detected in each sea storm, were used to assess a univariate statistical analysis. For simplicity, regions with similar values of a were grouped. Once the univariate variable is defined (Equation (1)), several statistical methods may be applied to carry out extreme value analysis for the marine environment, for instance, Boccotti [16,25], Goda [26] and Laface et al. [27]. The Generalized Pareto Distribution (GPD) was chosen, being suitable to describe a sample of independent, identically distributed random variables conditioned by a threshold. The Maximum Likelihood (ML) estimation procedure was used to select the probability density function parameters (the shape and scale parameters). The choice of the threshold can be based on two empirical techniques (Coles et al. [28]): one is based on the interpretation of the mean excess (ME) plot, which should be approximately linear in the proximity of the appropriate threshold, and the other is based on the stability of shape and scale parameter of GPD in the vicinity of the threshold. The latter method was used, trying to maintain a low threshold where possible.
from back-transformation in the standard space. The dashed line highlights the slope of the "limit space" and the red line is the structure variable limit over which failure occurs. In both spaces also the assumed joint probability density function of UH and U and of Hs and ζ is drawn.
The spatial structure variable r is therefore the most likely combination of Hs and ζ that determines the maximum flooding in that particular pixel. During each sea storm, the declustering is carried out identifying the maximum r, as shown in Figures 4 and 5.  The maximum r for each storm provides the critical pairs of Hs and ζ for each location, i.e., for each value of the coefficient a. In theory, these pairs could be used iteratively to define more approximate joint statistics and assuming the actual correlation structure. However, we consider this step unnecessary.
The samples, r, detected in each sea storm, were used to assess a univariate statistical analysis. For simplicity, regions with similar values of a were grouped. Once the univariate variable is defined (Equation (1)), several statistical methods may be applied to carry out extreme value analysis for the marine environment, for instance, Boccotti  Pareto Distribution (GPD) was chosen, being suitable to describe a sample of independent, identically distributed random variables conditioned by a threshold. The Maximum Likelihood (ML) estimation procedure was used to select the probability density function parameters (the shape and scale parameters). The choice of the threshold can be based on two empirical techniques (Coles et al. [Error! Reference source not found.]): one is based on the interpretation of the mean excess (ME) plot, which should be approximately linear in the proximity of the appropriate threshold, and the other is based on the stability of shape and scale parameter of GPD in the vicinity of the threshold. The latter method was used, trying to maintain a low threshold where possible.

Case Study
The Venetian littoral faces the North Adriatic Sea and is characterized by two main wind (and correspondingly wave) regimes, Bora and Scirocco, which blow from the Northeast and Southeast, respectively (

Case Study
The Venetian littoral faces the North Adriatic Sea and is characterized by two main wind (and correspondingly wave) regimes, Bora and Scirocco, which blow from the Northeast and Southeast, respectively (Umgiesser et al. [29]). Waves up to 6 m and high storm surges are responsible for the flooding of this coastal area as well as Venice.
Recently, two extreme events occurred along the Venetian littoral in 2018 (Vaia storm described in [30,31]) and in 2019 ( [32,33]). The first event was characterized by an extreme significant wave height of 5.92 m and a sea level of 1.56 m relative to the local tidal reference (named ZMPS, Zero Mareografico Punta Salute), conversely, the 2019 event was characterized by an extreme sea level of 1.82 m ZMPS and a wave height of 3.5 m. Both events caused localized erosion and marine ingression along the Venetian coast.
Caorle is a coastal town located on the Northeast side of the Venetian littoral and is one of the 20 coastal cells identified by Ruol et al. [34]. The Caorle coastline is 5 km long and is confined between the lagoon's inlet of Falconera to the North and the mouth of the Livenza river to the South, both protected with jetties ( Figure 6). The cell can be subdivided into three main parts: (i) "Spiaggia di Levante" to the Northeast; (ii) "Murazzi" in the central part; and (iii) "Spiaggia di Ponente" to the Southwest. The historic town is in the central part, protected by a sea wall that stabilizes the shoreline and mitigates the risk of flooding. North of this sea wall, there is a cusp, where the historic church named "Chiesa della Madonna dell'Angelo" is situated. The entire cell is urbanized and no system of dunes is present. The economy is mainly based on tourism (∼4,300,000 visitors in 2018, i.e., the 9th city in Italy for the number of presences) and fishing activities. Natural subsidence affects the area, with a rate of −2 to −4 mm/year [35].
The Caorle coastline is 5 km long and is confined between the lagoon's inlet of Falconera to the North and the mouth of the Livenza river to the South, both protected with jetties ( Figure 6). The cell can be subdivided into three main parts: (i) "Spiaggia di Levante" to the Northeast; (ii) "Murazzi" in the central part; and (iii) "Spiaggia di Ponente" to the Southwest. The historic town is in the central part, protected by a sea wall that stabilizes the shoreline and mitigates the risk of flooding. North of this sea wall, there is a cusp, where the historic church named "Chiesa della Madonna dell'Angelo" is situated. The entire cell is urbanized and no system of dunes is present. The economy is mainly based on tourism (∼4,300,000 visitors in 2018, i.e., the 9th city in Italy for the number of presences) and fishing activities. Natural subsidence affects the area, with a rate of −2 to −4 mm/year [Error! Reference source not found.].

Available Data
For the area under investigation, data and information were collected, harmonized and stored in a Coastal GIS [34]. The topographic data, useful for the present study, The annual mean sea levels ζ MSL for the Venetian littoral [36] are available from 1870 to 2017 measured at "Punta della Salute" station, situated approximately in front of the St. Mark square in the city of Venice. The ζ MSL is the long-term trend, described as a function of time, whereas ζ* (defined as ζ* = ζ − ζ MSL ) is a stationary stochastic variable.

Results
Following the proposed methodology, a dataset of 90 flooding maps is assessed, considering an offshore significant wave height H s ranging from 0 m to 7 m and ζ* ranging from 0 m to 2.5 m.
The SSW model is used to evaluate each map. The grid is formed by regular cells (dx = dy = 1 m) and covers an area of~1.5 km × 4.9 km (7.5 × 10 6 cells); the Gauckler-Strickler coefficient Ks used is 50 m 1/3 /s. Two different boundary conditions were considered: one at the shoreline, i.e., sea levels and significant wave heights, and the other at the two mouths (Livenza and Falconera), i.e., only the sea levels. Each simulation covers only 6 h and the flooding maps include the maximum flooded depth reached during the simulations.
These resulting maps were used to evaluate the coefficient a and the r LIM value of the structure variable, considering h LIM = 0.2 m, i.e., a "nuisance flooding" (Moftakhari et al. [37]). This type of flooding refers to low levels of inundation that do not cause notable threats to people or extensive damages, but it can disturb daily activities, add strain on infrastructures (roads, sewers, drainage systems, etc.) and cause minor damages to private properties. The annual mean sea levels ζMSL for the Venetian littoral [Error! Reference source not found.] are available from 1870 to 2017 measured at "Punta della Salute" station, situated approximately in front of the St. Mark square in the city of Venice. The ζMSL is the longterm trend, described as a function of time, whereas ζ* (defined as ζ* = ζ − ζMSL) is a stationary stochastic variable.

Results
Following the proposed methodology, a dataset of 90 flooding maps is assessed, considering an offshore significant wave height Hs ranging from 0 m to 7 m and ζ* ranging from 0 m to 2.5 m.
The SSW model is used to evaluate each map. The grid is formed by regular cells (dx = dy = 1 m) and covers an area of ~1.5 km × 4.9 km (7.5 × 10 6 cells); the Gauckler-Strickler coefficient Ks used is 50 m 1/3 /s. Two different boundary conditions were considered: one at the shoreline, i.e., sea levels and significant wave heights, and the other at the two mouths (Livenza and Falconera), i.e., only the sea levels. Each simulation covers only 6 h and the flooding maps include the maximum flooded depth reached during the simulations.
These resulting maps were used to evaluate the coefficient a and the rLIM value of the structure variable, considering hLIM = 0.2 m, i.e., a "nuisance flooding" (Moftakhari et al. [Error! Reference source not found.]). This type of flooding refers to low levels of inundation that do not cause notable threats to people or extensive damages, but it can disturb daily activities, add strain on infrastructures (roads, sewers, drainage systems, etc.) and cause minor damages to private properties.     The time series of Hs and ζ* (introduced above) covers 31.5 years. From the time series of Hs, following the methodology described above, the storm identification is performed considering hCRIT equal to 1 m and selecting tCRIT equal to 12 h from the analysis of the autocorrelation function. The total number of independent storms identified is 1057, i.e., ~34 events per year. Four samples r can be derived, corresponding to the different zones highlighted in Figure 7, in order to establish coastal hazard maps.
The samples, r, are analyzed with a univariate approach, fitting the data to a Generalized Pareto Distribution (GPD). The MLE procedure is used to select the probability density function parameters. The threshold is selected in order to have the same number of extremes for each coefficient a.
The evaluation of the four statistical analyses allows to assess the hazard map through the evaluation of the probability P{r > rLIM}, i.e., in this case the probability that a The time series of Hs and ζ* (introduced above) covers 31.5 years. From the time series of Hs, following the methodology described above, the storm identification is performed considering h CRIT equal to 1 m and selecting ∆t CRIT equal to 12 h from the analysis of the autocorrelation function. The total number of independent storms identified is 1057, i.e., 34 events per year. Four samples r can be derived, corresponding to the different zones highlighted in Figure 7, in order to establish coastal hazard maps.
The samples, r, are analyzed with a univariate approach, fitting the data to a Generalized Pareto Distribution (GPD). The MLE procedure is used to select the probability density function parameters. The threshold is selected in order to have the same number of extremes for each coefficient a.
The evaluation of the four statistical analyses allows to assess the hazard map through the evaluation of the probability P{r > r LIM }, i.e., in this case the probability that a flooding level equal to 0.2 m is exceeded. Figure 9 shows the hazard map for the Caorle area, expressed as a function of the return period T R . ries of Hs, following the methodology described above, the storm identification is performed considering hCRIT equal to 1 m and selecting tCRIT equal to 12 h from the analysis of the autocorrelation function. The total number of independent storms identified is 1057, i.e., ~34 events per year. Four samples r can be derived, corresponding to the different zones highlighted in Figure 7, in order to establish coastal hazard maps.
The samples, r, are analyzed with a univariate approach, fitting the data to a Generalized Pareto Distribution (GPD). The MLE procedure is used to select the probability density function parameters. The threshold is selected in order to have the same number of extremes for each coefficient a.
The evaluation of the four statistical analyses allows to assess the hazard map through the evaluation of the probability P{r > rLIM}, i.e., in this case the probability that a flooding level equal to 0.2 m is exceeded. Figure 9 shows the hazard map for the Caorle area, expressed as a function of the return period TR.

Discussion on the October 2018 and November 2019 Events
The stormy event that occurred on 28 October 2018 developed as a consequence of an explosive cyclogenesis in the Western Mediterranean Sea (Cavaleri et al. [Error! Reference source not found.]). This led to intense winds blowing from the Southeast (Scirocco) that caused severe damages in the Dolomites area. Extreme waves and high-water levels were originated in the Northern Adriatic Sea, and the recorded significant wave height was up to 6 m.
On 12 November 2019, a severe stormy event occurred in the Venetian littoral, characterized by high waves (Hs up to 3.5 m) and an extreme sea level. The water level reached 1.82 m ZMPS at 22.50 UTC (measured at the oceanographic tower), the second-highest value ever measured in this area. This level was caused by a combination of an astronomical tidal peak with a severe storm surge generated by a strong wind (up to 30 m/s) and a sudden pressure drop down to 987 hPa. The wind direction quickly turned from 100° N to 230° N. During this event, several failures occurred along the Venetian littoral. Venice

Discussion on the October 2018 and November 2019 Events
The stormy event that occurred on 28 October 2018 developed as a consequence of an explosive cyclogenesis in the Western Mediterranean Sea (Cavaleri et al. [31]). This led to intense winds blowing from the Southeast (Scirocco) that caused severe damages in the Dolomites area. Extreme waves and high-water levels were originated in the Northern Adriatic Sea, and the recorded significant wave height was up to 6 m.
On 12 November 2019, a severe stormy event occurred in the Venetian littoral, characterized by high waves (Hs up to 3.5 m) and an extreme sea level. The water level reached 1.82 m ZMPS at 22.50 UTC (measured at the oceanographic tower), the second-highest value ever measured in this area. This level was caused by a combination of an astronomical tidal peak with a severe storm surge generated by a strong wind (up to 30 m/s) and a sudden pressure drop down to 987 hPa. The wind direction quickly turned from 100 • N to 230 • N. During this event, several failures occurred along the Venetian littoral. Venice city was completely flooded, causing damages to historical heritage, houses, and boats. The Venetian coast was affected by localized erosions and even by coastal flooding phenomena that caused damages to many touristic and productive facilities.
The left of Figure 10 shows the sea levels and the significant wave heights measured at the CNR oceanographic tower during these two extreme events. Two instants are highlighted with stars that indicate the chosen sample for an area with a equal to 0.17. The right of Figure 10 shows the development of the two events in the Hsζ plane and the limit state for a pixel characterized by a = 0.17 and r LIM = 1.1 m. The two events caused the failure since both exceeded the limit state.
The statistical analysis carried out on the structure variables, built with different a values and presented in the previous paragraph, allows for estimating the return period of the two events. The results are summarized in Table 1 in terms of chosen couples, values of the structure variable for the events (r 0 ), exceedance probability P (r > r 0 ), and return period T R (where λ is the number of events per year). Table 1 shows also, in the last three rows, the results carried out with three more traditional sampling procedures. In detail, the following samples were taken into account: (S1) Maximum H s and simultaneous value of ζ*; (S2) Maximum ζ* and simultaneous value of H s ; (S3) Maxima of H s and ζ* during the storm regardless of their concomitance. For these samples, copulas are used to describe the dependence structure associated with the joint distribution of the two variables. A thorough introduction to Copula modelling and a large selection of the most common families are provided in [38,39]. Similar to the analysis of De Michele et al. [40], the Gumbel-Hougaard distribution, i.e., an Archimedean extreme value copula, is selected as the best candidate to model the data. To estimate the model parameters θ, the maximum likelihood estimation method (MLE) guarantees that the dataset is the most probable under the assumed statistical model. The 2018 event was characterized by extreme waves and high sea levels; conversely, the 2019 event was characterized by extreme sea levels and high waves. For the 2018 event, the return period (T R ) based on sample S1 is very high. Analyzing the same event considering sample S2, the associated T R is low and the exceptionality is not confirmed. For the 2019 event, the T R based on sample S2 is very high and the T R based on sample S1 low. Both events appear to be very rare with the analysis on sample S3, assuming that the maxima occur at the same time in the event.
Water 2021, 13, x FOR PEER REVIEW 12 of 17 city was completely flooded, causing damages to historical heritage, houses, and boats. The Venetian coast was affected by localized erosions and even by coastal flooding phenomena that caused damages to many touristic and productive facilities.
The left of Figure 10 shows the sea levels and the significant wave heights measured at the CNR oceanographic tower during these two extreme events. Two instants are highlighted with stars that indicate the chosen sample for an area with a equal to 0.17. The right of Figure 10 shows the development of the two events in the Hs-ζ plane and the limit state for a pixel characterized by a = 0.17 and rLIM = 1.1 m. The two events caused the failure since both exceeded the limit state. The statistical analysis carried out on the structure variables, built with different a values and presented in the previous paragraph, allows for estimating the return period of the two events. The results are summarized in Table 1 in terms of chosen couples, values of the structure variable for the events (r0), exceedance probability P (r > r0), and return period TR (where λ is the number of events per year). Table 1 shows also, in the last three rows, the results carried out with three more traditional sampling procedures. In detail, the following samples were taken into account: (S1) Maximum Hs and simultaneous value of ζ*; (S2) Maximum ζ* and simultaneous value of Hs; (S3) Maxima of Hs and ζ* during the storm regardless of their concomitance. For these samples, copulas are used to describe the dependence structure associated with the joint distribution of the two variables. A thorough introduction to Copula modelling and a large selection of the most common families are provided in [Error! Reference source not found.,Error! Reference source not found.]. Similar to the analysis of De Michele et al. [Error! Reference source not found.], the Gumbel-Hougaard distribution, i.e., an Archimedean extreme value copula, is selected as the best candidate to model the data. To estimate the model parameters , the maximum likelihood estimation method (MLE) guarantees that the dataset is the most probable under the assumed statistical model. The 2018 event was characterized by extreme waves and high sea levels; conversely, the 2019 event was characterized by extreme sea levels and high waves. For the 2018 event, the return period (TR) based on sample S1 is very high. Analyzing the same event considering sample S2, the associated TR is low and the exceptionality is not confirmed. For the 2019 event, the TR based on sample S2 is very high and the TR based on sample S1 low. Both events appear to be very rare with the analysis on sample S3, assuming that the maxima occur at the same time in the event.  From the hazard map presented in Figure 9, the point x = 2.5 km, y = 0.3 km, i.e., the one corresponding to picture no. 2 of Figure 11, has a = 0.15 and r LIM = 1.4 m. From Table 1, the r 0 of November 2019 is larger and therefore flooding is expected. The associated return period in this area is T R = 47 years. In fact, marine ingression occurred along the Caorle coastline during November 2019: both the beaches were completely flooded and the waves reached the church at the cusp (pictures in Figure 11). From the hazard map presented in Figure 9, the point x = 2.5 km, y = 0.3 km, i.e., the one corresponding to picture no. 2 of Figure 11, has a = 0.15 and rLIM = 1.4 m. From Table  1, the r0 of November 2019 is larger and therefore flooding is expected. The associated return period in this area is TR = 47 years. In fact, marine ingression occurred along the Caorle coastline during November 2019: both the beaches were completely flooded and the waves reached the church at the cusp (pictures in Figure 11). The SSW model is used to simulate this event and it is initialized with the same grid and the same conditions described before. The simulation covers only 4 h around the peak of sea level. Figure 11 shows the maximum flooded level reached during the simulation. Since the flooding occurred during nighttime, only some qualitative information and some photographs, taken the morning after, are available ( Figure 11). However, these photographs show the damages and the residual flooding on the beach and at the church. The simulation results are consistent with this information. The SSW model is used to simulate this event and it is initialized with the same grid and the same conditions described before. The simulation covers only 4 h around the peak of sea level. Figure 11 shows the maximum flooded level reached during the simulation. Since the flooding occurred during nighttime, only some qualitative information and some photographs, taken the morning after, are available ( Figure 11). However, these photographs show the damages and the residual flooding on the beach and at the church. The simulation results are consistent with this information.

Conclusions
In many low-lying coastal regions, flooding hazard assessment is critically affected by the correct evaluation of the simultaneous occurrence of the significant wave height H s and the water level ζ. In turn, the bivariate statistics is significantly affected by the methodology used for the selection of the independent sample pair H s , ζ during a sea storm.
Following Mazas and Hamm [2], the suggested sampling procedure is to define a structure variable that considers the combination of H s and ζ that is the most critical for flooding. This study proposes a procedure to find a spatial structure variable, that is a linear combination of water level and wave height, specific for each location, based on the slope of the failure function at the design point.
As a consequence, the statistical analysis is specific to each point of the area at risk. For instance, for low-lying points located far from a coastal dune system, ζ is more important than H s , and hence the relative flooding occurrence is mainly affected by the maximum water levels statistics. Conversely, for an area located just behind a seawall where overtopping is relevant, the role played by H s is dominant, and the relative flooding occurrence is strongly affected by the maximum significant wave height statistics.
The proposed procedure was applied to a stretch of the Venetian littoral that was affected by marine ingression during two recent stormy events (occurred in 2018 and 2019). A map showing the chosen combination of the variables is presented together with the hazard map for the investigated coastline.
The paper draws attention to the (apparently trivial) need for a consistent procedure for sampling multivariate events from sea measurements. Based on real examples, it is shown that the method used to sample the bivariate event formed by H s and ζ significantly affects the return period of the critical storm and the sampling method should therefore be considered as an important aspect for the evaluation of the flooding hazard. The comparison between traditional sampling methods and the proposed one highlights that the latter gives more consistent results in terms of return periods.
The method can be implemented to coastal stretches subject to rather uniform environmental conditions. For instance, for the Venetian littorals, it was assessed that the length of each stretch is of order 3-12 km.
As further work, it is envisaged that application of the procedure to many real and schematic cases could point out some a priori guidelines to select the proper combination of water level and significant wave height only based on the morphological characteristics of the site, thus defining a priori the proper structure variable for the location of interest.