Deriving the 100-Year Total Water Level around the Coast of Corsica by Combining Trivariate Extreme Value Analysis and Coastal Hydrodynamic Models

: As low-lying coastal areas can be impacted by ﬂooding caused by dynamic components that are dependent on each other (wind, waves, water levels—tide, atmospheric surge, currents), the analysis of the return period of a single component is not representative of the return period of the total water level at the coast. It is important to assess a joint return period of all the components. Based on a semiparametric multivariate extreme value analysis, we determined the joint probabilities that signiﬁcant wave heights (Hs), wind intensity at 10 m above the ground (U), and still water level (SWL) exceeded jointly imposed thresholds all along the Corsica Island coasts (Mediterranean Sea). We also considered the covariate peak direction (Dp), the peak period (Tp), and the wind direction (Du). Here, we focus on providing extreme scenarios to populate coastal hydrodynamic models, SWAN and SWASH-2DH, in order to compute the 100-year total water level (100y-TWL) all along the coasts. We show how the proposed multivariate extreme value analysis can help to more accurately deﬁne low-lying zones potentially exposed to coastal ﬂooding, especially in Corsica where a unique value of 2 m was taken into account in previous studies. The computed 100y-TWL values are between 1 m along the eastern coasts and a maximum of 1.8 m on the western coast. The calculated values are also below the 2.4 m threshold recommended when considering the sea level rise (SLR). This highlights the added value of performing a full integration of extreme offshore conditions, together with their dependence on hydrodynamic simulations for screening out the coastal areas potentially exposed to ﬂooding


Introduction
Coastal flooding often occurs because of several dynamic components that are dependent on each other (wind, waves, water levels-tide, atmospheric surge, currents); hence, the analysis of the return period of a single component is not representative of the return period of the total water level at the coast [1].As such, it is important to estimate the joint return period, taking into account the dependency between all of the components.Marcos et al. [2] even show that return periods of extreme sea levels are underestimated (by a factor of 2 or higher) in 30% of the coasts (at global scale), if the dependency is neglected.
The dependency can also be based on a spatial linkage assumption.Galiatsatou and Prinos [3] took into account the degree of spatial dependence between the sites, in order to estimates extreme storm surges along the Dutch coast on the North Sea.Some authors considered extreme quantiles of metocean data (significant wave height, wave period, and wind speed) estimated jointly, using a subset exceeding appropriately defined thresholds [4].Moreover, several joint probability approaches applied to sea conditions have been used by [5][6][7][8][9], for example.The reader might also be interested in the review proposed by [10].One of the limitations encountered within the classical joint probability approaches are related to the restrictive assumptions regarding the dependence structure when extrapolating to extremes.Idier et al. [11] mentioned that one of the difficulties with such approaches could be the number of offshore parameters, which can complicates the detection of the critical conditions.Heffernan and Town [12] developed a semiparametric approach (hereinafter called H&T04 methods), which overcomes the limitations imposed by dimensionality.This approach was first applied to air pollution data.Since then, several authors applied the H&T04 methods to metocean offshore conditions: [13][14][15][16].
In the present study, we are interested in coastal flooding all around Corsica Island in the Mediterranean Sea.To our best knowledge, a multivariate extreme value analysis has never been implemented along these coasts, though some studies undertook statistical analysis over the Mediterranean Sea or around Corsica, e.g., [17].In particular, in 2004, the Western European Union funded a Wind and Wave Atlas of the Mediterranean Sea, providing both bivariate and univariate statistics, and spatial distribution of statistical quantities for wind and wave parameters [18].This atlas was built 20 years ago with a 10-year dataset consisting of less than 1000 data points (spatial resolution of the outputs ~50 km), distributed throughout the whole Mediterranean Sea using dynamical models which have been updated since the beginning of the 21st century.Other initiatives in this basin should be noted.For example, the WaveForUs (Wave climate and coastal circulation Forecasts for public Use) platform provides 3-day meteorological and sea state forecasts in the Mediterranean Sea since 2013 with a 0.15 • spatial resolution [19].To date, these archives are too short to be used as such in local coastal flooding assessments.
To be more accurate at local scale, a solution could be to apply dynamical downscaling using a complex nested model chain.However, this is particularly delicate in the Mediterranean Sea because of specific issues mentioned in [20]: especially the orography (for winds and pressure), and the complex bathymetry, and limited fetch extension (for waves).Several authors proposed an alternative in order to reduce the complexity and the computational cost of such a numerical chain [21,22].Indeed, based on the forcing event selection, which consist in: (1) determining the joint probabilities that significant wave heights Hs, wind intensity at 10 m above the ground U, and still water level SWL exceed jointly imposed thresholds; (2) defining the extreme combinations; (3) populating the coastal hydrodynamic models at a local scale with the extreme combinations as input.For instance, the authors in reference [23] used dozens of scenarios of offshore forcing conditions (significant wave height and total water level) associated with the same joint probability of exceedance (here with a 100-year return period) as inputs of their modeling chain.The appealing feature of this approach is its computational simplicity and efficiency, because it is only based on a limited number of scenarios, i.e., on a limited number of long-running numerical simulations.
From an operational point of view, our primary motivation is based on the fact that, since 2013, a unique value of 2 m relative to the French national topographic reference (NGF) was taken into account to map entire zones potentially exposed to coastal flooding in Corsica (e.g., http://carto.geo-ide.application.developpement-durable.gouv.fr/429/risques_naturels_02a.map (accessed on 2 September 2020)).The reader may note that all data and results presented in this study are relative to the NGF.For the local stakeholders and coastal managers, it is important to assess if this 2 m value is relevant in their region [24].The main purpose of the present study is to show how multivariate extreme value analysis can help to better define low-lying zones potentially exposed to coastal flooding.
In this paper, we have first described the site (location, and context), and data used in Section 2. We then exposed the methods applied to constitute extreme scenarios based on triplets (Hs; SWL; U) for significant wave height, still water level, and wind speed.The strategy used to propagate offshore conditions to the coastline has also been presented in this section.In Section 3, we have applied the proposed strategy all around Corsica, considering extreme scenarios under both, current and climate change conditions.The results are presented in Section 4. Finally, a discussion of our results, and the limitations of the method are presented in Section 5.

Site Description and Data Overview
The study site is Corsica Island, located in the Mediterranean Sea (Figure 1).Since 2013, the approach implemented for the whole French Mediterranean coastline consists in mapping coastal sectors whose altitude is below 2 m.These sectors are considered potentially exposed to coastal flooding.Indeed, in past studies, the 100-year return period event was computed using the Sète tide gauge located at least 400 km from the coasts studied (red symbol in Figure 1).Moreover, the analysis did not take into account the combination of wind, waves, and water level as a joint return period.For coastal management, it is necessary to assess if this 2 m value is valid in our study area.
value analysis can help to better define low-lying zones potentially exposed to coastal flooding.
In this paper, we have first described the site (location, and context), and data used in Section 2. We then exposed the methods applied to constitute extreme scenarios based on triplets (Hs; SWL; U) for significant wave height, still water level, and wind speed.The strategy used to propagate offshore conditions to the coastline has also been presented in this section.In Section 3, we have applied the proposed strategy all around Corsica, considering extreme scenarios under both, current and climate change conditions.The results are presented in Section 4. Finally, a discussion of our results, and the limitations of the method are presented in Section 5.

Site Description and Data Overview
The study site is Corsica Island, located in the Mediterranean Sea (Figure 1).Since 2013, the approach implemented for the whole French Mediterranean coastline consists in mapping coastal sectors whose altitude is below 2 m.These sectors are considered potentially exposed to coastal flooding.Indeed, in past studies, the 100-year return period event was computed using the Sète tide gauge located at least 400 km from the coasts studied (red symbol in Figure 1).Moreover, the analysis did not take into account the combination of wind, waves, and water level as a joint return period.For coastal management, it is necessary to assess if this 2 m value is valid in our study area.Figure 1b highlights that the island is under the influence of different wind regimes coming from far as the Libeccio (can cross more than 1000 km) or coming from closer as the Mezzogiorno (~10 km).The waves induced by the wind can be swell when they cross long distance or wind waves if they stay close to their location of generation.Because of the distance from other coasts, we can assume that the west coast of Corsica can be impacted by swells, whereas the south and east coast mostly by wind waves.
In this area, the tidal range is from 0.2 to 0.4 m as it is a microtidal regime [25].For this reason, we consider a direct approach using the SWL instead of storm surges.Indeed, Haigh et al. [26] showed that for long return periods, when the ratio of the tidal to nontidal component is lower than 2, direct and non-direct approaches present good agreement.Kergadallan [27] showed that this condition (<2) is observed for the French Mediterranean coasts.Figure 1b highlights that the island is under the influence of different wind regimes coming from far as the Libeccio (can cross more than 1000 km) or coming from closer as the Mezzogiorno (~10 km).The waves induced by the wind can be swell when they cross long distance or wind waves if they stay close to their location of generation.Because of the distance from other coasts, we can assume that the west coast of Corsica can be impacted by swells, whereas the south and east coast mostly by wind waves.
In this area, the tidal range is from 0.2 to 0.4 m as it is a microtidal regime [25].For this reason, we consider a direct approach using the SWL instead of storm surges.Indeed, Haigh et al. [26] showed that for long return periods, when the ratio of the tidal to non-tidal component is lower than 2, direct and non-direct approaches present good agreement.Kergadallan [27] showed that this condition (<2) is observed for the French Mediterranean coasts.
In order to manage our statistical analysis, validate data, and run hydrodynamic simulations, we used different products listed at Table 1.The National Oceanic and Atmospheric Administration (NOAA) provides the NWW3 MED hindcast [28,29] consisting of sea states, simulated using the third generation wave model Wavewatch 3 (WW3) with the parametrization proposed by [30].Offshore wind conditions are also provided.
The SWL hindcast was generated by the BRGM using the MARS-2DH model [31].This dataset is available from 1979 to 2010 and gives an assessment of the regional hydrodynamics based on tidal components and meteorological reanalysis CFSR [32].
The bathymetry and topography used after in the nearshore numerical modeling chain were provided by the Shom [33] and IGN.
We used in-situ data from Candhis and ISPRA (Table 1) to correct the Hs in the NWW3_MED hindcast using a linear correction (see Appendix A).This correction reduces the bias and is relevant for both mean and Hs > 4 m.For extreme Hs, the bias reduction is ~1 m.In addition, we used the in-situ datasets provided by the Shom and CNES to validate the MARS_MED_ BRGM hindcast (Appendix A), and in-situ data coming from tide gauges provided by ISPRA and Météo-France in order to validate wind speed and direction.

Statistical and Numerical Strategy
The purpose was to determine triplets (Hs; SWL; U) of a 100-year joint exceedance return period and their covariates (Tp, Dp, Du) to populate coastal hydrodynamic models.In order to perform the statistical analysis, we used MATLAB ® software.The statistical method was implemented using the WAFO toolbox [34] and developments by Guanche Garcia [35].The different steps are illustrated in Figure 2a and exposed here: Selection of events and data set constitution: build a sample with a large number of independent triplets (Hs; SWL; U); 2.
Fitting of the marginal probability laws for each variable Hs, SWL, and U using the Generalized Pareto Distribution (GPD) with bootstrap simulation for Confidence Interval (CI) [36]; 3.
Dependency model adjustment: (a) Between the extreme variables Hs, SWL, U, and based on the semiparametric approach described by [12].The readers may also consult [37] for detailed information and application of this approach; (b) Between Hs et Tp using the empirical conditional distribution of wave steepness (St) knowing Hs described by [38]; 4.
Running Monte-Carlo simulations using the marginal laws and dependency models of a very large number of combinations (Hs; SWL; U) and their covariates Tp, Dp, Du with the same statistical characteristics as the observed data; 5.
Determine the triplets (Hs; SWL; U) of the 100-year joint exceedance return period; 6.
Dynamical downscaling by running SWAN model for the 10 m resolution nested boxes (in orange in Figure 2b).SWAN is here forced by the outputs of step 6; 8.
Assess the total water level (TWL) at the shoreline adding contributions coming from the models: wave setup provided by SWAN and wind setup by SWASH-2DH.We applied the statistical analysis for points located in each of the red boxes shown in Figure 2b.The triplets obtained and their covariates finally forced the SWAN and SWASH-2DH models.For more information on the SWAN and SWASH-2DH open source models, the reader is invited to refer to [39][40][41][42][43][44][45].
Moreover, as mentioned at step 6, some safety margins are taken into account in the numerical simulations.Following the French regulation [46], the water level must account the Sea Level Rise (SLR) potentially due to climate change.To do so, 0.20 m is added to the water level of the 100-year event defined offshore.In order to consider a "2100" future scenario, 0.60 m has to be added to the 100-year return level defined offshore by following the French regulation guideline.The values recommended in the guideline are based on We applied the statistical analysis for points located in each of the red boxes shown in Figure 2b.The triplets obtained and their covariates finally forced the SWAN and SWASH-2DH models.For more information on the SWAN and SWASH-2DH open source models, the reader is invited to refer to [39][40][41][42][43][44][45].
Moreover, as mentioned at step 6, some safety margins are taken into account in the numerical simulations.Following the French regulation [46], the water level must account the Sea Level Rise (SLR) potentially due to climate change.To do so, 0.20 m is added to the water level of the 100-year event defined offshore.In order to consider a "2100" future scenario, 0.60 m has to be added to the 100-year return level defined offshore by following the French regulation guideline.The values recommended in the guideline are based on the Fifth Assessment Report (AR5), chapter 13 [47] published by the Intergovernmental Panel on Climate Change (IPCC).Note that the value of 0.60 m is in agreement with the mean SLR value estimated for the Mediterranean basin in [48,49], when considering the pessimistic scenario of Representative Concentration Pathway (RCP) 8.5.Moreover, [46] recommends taking into account uncertainties adding a margin of 0.25 m to the 100-year return level if all the uncertainties cannot be estimated.There are many sources of uncertainties, and even if many of them showed small-to-moderate degree of uncertainty (confidence interval given with marginals, correction on Hs), and that parametric tests regarding the statistical approach showed moderate impact as well (diagnosis on thresholds in H&T04), residual uncertainties remain.These residual uncertainties correspond to altimetric variations of Earth's crust, local vertical ground motions (e.g., due to groundwater pumping), impacts of remote contributions in a semi-enclosed basin such as the Mediterranean Sea or complex oceanic processes in the Strait of Gibraltar.To account for them, a margin of 0.25 m is applied following the French regulators' recommendations of [46].We have considered two cases: (i) Present-day conditions: 100-year triplets/covariates + 0.2 m SLR + 0.25 m uncertainties; (ii) "2100" future conditions: 100-year triplets/covariates + 0.6 m SLR + 0.25 m uncertainties.

Multivariate Extreme Value Analysis
In this section, we provide further implementation details on the steps 1-5 of the procedure described in Section 2.2, which aims at estimating the extreme offshore conditions to be used as forcing inputs of the SWAN-SWASH 2DH modeling chain (see Figure 2a).These steps are illustrated with the box around the Gulf of Porto (Figure 2b)

Selection of Events and Dataset Set Up (Step 1)
The first step of the procedure consists in selecting independent combinations (Hs; SWL; U).To do so, we first performed a block maxima selection to define each wave event (statistically independent) [50].Indeed, based on a 40-year historical research in archives, press and in-situ data; and in consultation with the deconcentrated services of the state, we assumed that Hs is the primary driver of coastal flooding around the coast of Corsica.Following [23,51], we defined a block as a 3-day period with 1.5 days between each Hs peak in order to ensure the independence between each peak.Then, for each Hs peak, the SWL and U maxima were sought in a 12-h window centered on the Hs peak.This 12 h window allows the data of SWL to be scanned over approximately one tidal cycle.Each Hs peak is associated to the corresponding covariates, peak period Tp and peak direction Dp.Likewise, each value of U is associated with its direction Du.Thus, several combinations (Hs; Tp; Dp; SWL; U; Du) corresponding to 30 years of common data with ~100 events per year were selected in 7 different zones around Corsica.

Marginal Distributions (Step 2)
Figure 3 illustrates an example of the marginal probability distribution fitting F i (x) of each variable X i for (Hs; SWL; U) using the GPD defined as: where u i is a high threshold to choose, σ i and ξ i are the GPD parameters.
we have reduced them by using time series from numerical models with a common period of 30 years.

Dependency Modeling and Monte Carlo Simulations
First, the variables  for (Hs; SWL; U) are transformed onto the Gumbel scales using the standard probability integral transformation to give the variables  .If  denotes a vector for all variables excluding  , a multivariate non-linear regression model is applied to  | > ν: where  and  are vectors of parameters, ν a threshold, and W is a vector of residuals.Each fitted model describes the dependence between the variables (except  ).The models were adjusted using the maximum likelihood method and assuming that the residuals W are Gaussian.Using the diagnostic tools described by H&T04 [12], we here have selected ν = 0.95 (expressed as a probability of non-exceeding this threshold).Examples of the parameters a and b, and of the selected thresholds ν are given in Appendix D for the western and eastern region of Corsica.
Once the parameters of the dependency relationship were estimated, a 10,000-year period was simulated applying a Monte Carlo method.We obtained a set of fictive sea states, wind, and SWL conditions as illustrated in blue in Figure 4.The simulated variables Yi are then transformed back to the original scales.Finally, the output is a large sample of sea and wind conditions (119,300 for the dataset NWW3 560/MARS 8) where at least one The data points considered are located close to the Gulf of Porto (see Figure 3a).Here the point NWW3 560 provides U (m/s) and Hs (m), and the point MARS 8 gives the SWL (m).For SWL, the threshold u = 0.3 m, the 100-year return level is 0.43 m.The threshold is set to 5.32 m, and the 100-year Hs is 8.016 m for Hs, and set to 16.92 m/s for U leading to a 100-year U of 22.55 m/s.Different thresholds have been chosen for each zone around Corsica (see Appendix B).To do so, we applied several analysis of mean residual life plots and analysis of quantitative values resulting from statistical tests Chi-squared, Kolmogorov-Smirnov, or quantile-quantile graphs, based on [50] (not shown).
In addition, three methods for estimating the GPD parameters have been tested: the method of moments (MOM), the weighted method of moments (PWM), and the maximum likelihood method (MLE) [52].Here (Figure 3b,d), the MOM was used to determine the parameters of the laws for Hs and U and the PWM method for SWL (Figure 3c), because these methods were the most conservative.Directions Dp and Du associated to Hs and U, respectively, are also reported in color.In this zone, the extremes are globally due to winds and waves with directions between 250 and 300 • .Note that we chose to keep all directions Dp and Du, not only extremes coming from a given sector, because: (i) critical Dp for present day conditions are expected to change in the future due to future geomorphological changes on sandy coasts; (ii) only keeping data coming from a "given" sector can considerably reduce the sample size and increase the uncertainty in the statistical methods (marginals and dependence model).For instance, only 17% of data are kept if the dataset is split following directions between 200 • and 250 • for the NWW3 465/MARS 26 dataset (see Appendix C).Finally, the uncertainty associated with the distribution fitting is translated into the form of confidence intervals given using the parametric bootstrap [53].Uncertain-ties related to the duration of the dataset may exist, but in our case, we have reduced them by using time series from numerical models with a common period of 30 years.First, the variables X i for (Hs; SWL; U) are transformed onto the Gumbel scales using the standard probability integral transformation to give the variables Y i.If Y −i denotes a vector for all variables excluding Y i , a multivariate non-linear regression model is applied to where a and b are vectors of parameters, ν a threshold, and W is a vector of residuals.Each fitted model describes the dependence between the variables (except Y i ).The models were adjusted using the maximum likelihood method and assuming that the residuals W are Gaussian.Using the diagnostic tools described by H&T04 [12], we here have selected ν = 0.95 (expressed as a probability of non-exceeding this threshold).Examples of the parameters a and b, and of the selected thresholds ν are given in Appendix D for the western and eastern region of Corsica.
Once the parameters of the dependency relationship were estimated, a 10,000-year period was simulated applying a Monte Carlo method.We obtained a set of fictive sea states, wind, and SWL conditions as illustrated in blue in Figure 4.The simulated variables Y i are then transformed back to the original scales.Finally, the output is a large sample of sea and wind conditions (119,300 for the dataset NWW3 560/MARS 8) where at least one variable is extreme (i.e., exceeding a defined threshold) with respect to the individual marginal distributions, as well as the dependency relationship between variables.

Joint Probability Contours
Since we have our large sample obtained by Monte Carlo simulation, we can access the joint exceedance isocontours [22].These 3D isocontours (x; y; z) are in the space of the three variables (Hs; SWL; U) where each point of the contour has the same probability of being exceeded: where λ is the number of events per year (in this case λ = 122) and T the return period considered in years.Figure 5a gives a 2D representation of the isocontours.The values on the isocontours are the SWL values from 0.2 to 0.4 m.From Figure 5a, it is then possible to identify combinations (Hs; SWL; U) of 100-year joint exceedance return periods.Figure 5b is a 3D representation of 100-year joint exceedance isocontours for the dataset NWW3 560/MARS 8.The black dots represent the combinations taken to force the hydrodynamic models.We have selected 34 scenarios for the dataset NWW3 560/MARS 8.These were selected along the isocontours by focusing preferably in the regions where the joint effect of the different drivers is the largest, i.e., where the convexity of the isocontours is the

Joint Probability Contours
Since we have our large sample obtained by Monte Carlo simulation, we can access the joint exceedance isocontours [22].These 3D isocontours (x; y; z) are in the space of the three variables (Hs; SWL; U) where each point of the contour has the same probability of being exceeded: where λ is the number of events per year (in this case λ = 122) and T the return period considered in years.Figure 5a gives a 2D representation of the isocontours.The values on the isocontours are the SWL values from 0.2 to 0.4 m.From Figure 5a, it is then possible to identify combinations (Hs; SWL; U) of 100-year joint exceedance return periods.Figure 5b is a 3D representation of 100-year joint exceedance isocontours for the dataset NWW3 560/MARS 8.The black dots represent the combinations taken to force the hydrodynamic models.We have selected 34 scenarios for the dataset NWW3 560/MARS 8.These were selected along the isocontours by focusing preferably in the regions where the joint effect of the different drivers is the largest, i.e., where the convexity of the isocontours is the largest.
where λ is the number of events per year (in this case λ = 122) and T the return period considered in years.Figure 5a gives a 2D representation of the isocontours.The values on the isocontours are the SWL values from 0.2 to 0.4 m.From Figure 5a, it is then possible to identify combinations (Hs; SWL; U) of 100-year joint exceedance return periods.Figure 5b is a 3D representation of 100-year joint exceedance isocontours for the dataset NWW3 560/MARS 8.The black dots represent the combinations taken to force the hydrodynamic models.We have selected 34 scenarios for the dataset NWW3 560/MARS 8.These were selected along the isocontours by focusing preferably in the regions where the joint effect of the different drivers is the largest, i.e., where the convexity of the isocontours is the largest.

Covariates
Once the (Hs, SWL, U) combinations are identified each Hs must be associated to their peak period (Tp).The peak period (Tp) was simulated following the approach described by [37] and adapted from [22].In this approach, a regression model of steepness (S t ), conditional on Hs is used: This model was used to simulate the steepness (and Tp) from Monte-Carlo simulations of Hs (blue dots in Figure 6).In order to associate Tp to Hs, we take values around the median of the simulated periods for each of the Hs considered.Table 2 gives the combinations Hs, Tp, U, and SWL selected from the NWW3 560/MARS 8 dataset, and representing the conditions offshore the Gulf of Porto.
For the remaining covariates, namely peak directions (Dp) associated to waves, and wind directions (Du), we have made the choice to impose the dominant direction(s) leading to extreme Hs and U. Figure 3b,d showed that Dp and Du associated to extreme waves and winds can be carried by one dominant sector (directions between 250 and 300 • ). Figure 7 illustrates wave (Figure 7b) and wind roses (Figure 7c) at points NWW3 512 and NWW3 544 (see locations in Figure 7a) giving the dominant peak directions (Dp), and wind directions (Du).These points give two dominant directions of extremes for Hs and U.The most penalizing values in terms of effect on the coast come from the west-south-west for the waves and south-west for the winds for the point NWW3 512.We forced the hydrodynamic models (SWASH-2DH and SWAN) with these directions.If there were two penalizing directions, regarding the exposure of the coastline, the directions were taken into account by running the hydrodynamic models twice depending on the sector of extreme Dp and Du.We did this double analysis for the point NWW3 544 as well and at every point presenting two penalizing directions.

S = (4)
This model was used to simulate the steepness (and Tp) from Monte-Carlo simulations of Hs (blue dots in Figure 6).In order to associate Tp to Hs, we take values around the median of the simulated periods for each of the Hs considered.Table 2 gives the combinations Hs, Tp, U, and SWL selected from the NWW3 560/MARS 8 dataset, and representing the conditions offshore the Gulf of Porto.

Results
Figure 8 shows the 100-year joint exceedance isocontours for all the points analyzed around Corsica.Different extreme regimes are here outlined: the dependence between the extreme triplets on the western part appear to be less marked than for the eastern part (as indicated by the convexity and range covered by the isocontours).This means that the 100-year triplets on the western coast deviate less from the 100-year return levels estimated using the marginal (i.e., without accounting for the dependence).This is clearly not the case on the eastern part as illustrated for instance by the Bastia analysis (black box) which shows large convexity of the isocontours.This 100-year joint exceedance isocontours' spatialization also highlights higher values of Hs (>6 m) and U (>18 m/s) offshore on the western part of the island compared to the eastern part of Corsica.This is in agreement with the different wind regimes (and induced waves) presented in Figure 1b.Indeed, on the west side of the island, wind and wave features (Hs > 6 m, U > 18 m/s, Tp > 12 s, 500 km < fetch < 1300 km) suggest that extreme scenarios may be driven by swells, whereas offshore conditions on the east part of Corsica mostly by wind waves and high

Results
Figure 8 shows the 100-year joint exceedance isocontours for all the points analyzed around Corsica.Different extreme regimes are here outlined: the dependence between the extreme triplets on the western part appear to be less marked than for the eastern part (as indicated by the convexity and range covered by the isocontours).This means that the 100-year triplets on the western coast deviate less from the 100-year return levels estimated using the marginal (i.e., without accounting for the dependence).This is clearly not the case on the eastern part as illustrated for instance by the Bastia analysis (black box) which shows large convexity of the isocontours.This 100-year joint exceedance isocontours' spatialization also highlights higher values of Hs (>6 m) and U (>18 m/s) offshore on the western part of the island compared to the eastern part of Corsica.This is in agreement with the different wind regimes (and induced waves) presented in Figure 1b.Indeed, on the west side of the island, wind and wave features (Hs > 6 m, U > 18 m/s, Tp > 12 s, 500 km < fetch < 1300 km) suggest that extreme scenarios may be driven by swells, whereas offshore conditions on the east part of Corsica mostly by wind waves and high still water levels.For example, on the eastern part of the island the combination (Hs = 2.2 m; Tp = 6.8 s; U = 12.1 s; SWL = 0.4 m), represented by a blue star in the red box in Figure 8, is identified as an extreme scenario to force the hydrodynamic models.The SWAN and SWASH-2DH simulations (steps 6-7) provide the contributions to compute total water levels (TWL) (step 8) at the shoreline for each offshore condition combination tested (see Section 2.2 for details on the methods).Figure 9 illustrates the distribution of TWL along the shoreline (a) for current conditions, and (b) for future conditions (2100).For current conditions, the static TWL along the shoreline in the study area is between 0.80 and 1.8 m.For "2100" future conditions, the values are between 1.2 and 2.2 m.Clearly, outside the area on the west of the Cap Corse in the north of the island, this value is hardly reached even under future conditions.Moreover, we note that values of TWL are higher on the west coast rather than the east coast, especially in the north west and south west, for both current and future conditions.The SWAN and SWASH-2DH simulations (steps 6-7) provide the contributions to compute total water levels (TWL) (step 8) at the shoreline for each offshore condition combination tested (see Section 2.2 for details on the methods).Figure 9 illustrates the distribution of TWL along the shoreline (a) for current conditions, and (b) for future conditions (2100).For current conditions, the static TWL along the shoreline in the study area is between 0.80 and 1.8 m.For "2100" future conditions, the values are between 1.2 and 2.2 m.Clearly, outside the area on the west of the Cap Corse in the north of the island, this value is hardly reached even under future conditions.Moreover, we note that values of TWL are higher on the west coast rather than the east coast, especially in the north west and south west, for both current and future conditions.
Finally, we tested the sensitivity of the results to the choices of the offshore conditions' combinations: this is generally below 0.05 m and can be considered not significant.Only a few combinations show differences greater than 0.05 m.In particular, Figure 10a shows that combinations 19 and 27 lead to the highest total water levels at the shoreline in the Gulf of Porto considering the specific directions Dp = 270 • and Du = 240 • .Moreover, the directions Dp and Du of waves and winds influence the scenarios leading to maximum water level at the shoreline.Indeed, when applying other directions (even with small changes) Dp = 240 • and Du = 240 • to force the SWAN and SWASH-2DH models, the scenario leading to the highest total water levels is different: combination 17 (Figure 10b).The differences are particularly marked in the southern Gulf of Porto where total levels at the shoreline are higher when applying Dp = 270 • and Du = 240 • (Figure 10c) instead of Dp = 240 • and Du = 240 • (Figure 10d).Finally, we tested the sensitivity of the results to the choices of the offshore conditions' combinations: this is generally below 0.05 m and can be considered not significant.Only a few combinations show differences greater than 0.05 m.In particular, Figure 10a shows that combinations 19 and 27 lead to the highest total water levels at the shoreline in the Gulf of Porto considering the specific directions Dp = 270° and Du = 240°.Moreover, the directions Dp and Du of waves and winds influence the scenarios leading to maximum water level at the shoreline.Indeed, when applying other directions (even with small changes) Dp = 240° and Du = 240° to force the SWAN and SWASH-2DH models, the scenario leading to the highest total water levels is different: combination 17 (Figure 10b).The differences are particularly marked in the southern Gulf of Porto where total levels at the shoreline are higher when applying Dp = 270° and Du = 240° (Figure 10c) instead of Dp = 240° and Du = 240° (Figure 10d).

Discussion and Conclusions
The key motivation of this study was to show how a multivariate extreme value analysis improve the estimate of low-lying zones potentially exposed to coastal flooding in Corsica.Indeed, in this zone, coastal sectors with altitudes below 2 m (and 2.4 m under "2100" future conditions), are considered potentially exposed to flooding.In order to asses

Discussion and Conclusions
The key motivation of this study was to show how a multivariate extreme value analysis improve the estimate of low-lying zones potentially exposed to coastal flooding in Corsica.Indeed, in this zone, coastal sectors with altitudes below 2 m (and 2.4 m under "2100" future conditions), are considered potentially exposed to flooding.In order to asses if this value was relevant or not, we have first determined the joint probabilities that significant wave heights (Hs), wind intensity at 10 m above the ground (U), and water level (WL) exceed jointly imposed thresholds by applying a semiparametric multivariate extreme value analysis.Covariate peak direction (Dp), the peak period (Tp), and the wind direction (Du) were also accounted for.Once the (Hs, Tp, SWL, U) combinations were identified, we ran the coastal hydrodynamic models, SWAN and SWASH-2DH forced by these extreme scenarios.For current conditions, the results show that the TWL all along the shoreline in the study area are between 0.80 and 1.8 m, compared to 2 m currently applied on every low-lying zone.For "2100" future conditions, the values are between 1.2 and 2.2 m, compared to 2.4 m currently applied on every low-lying zone.In conclusion, the value 2 m (2.4 m when considering future conditions) seems to be overestimated regarding, the methods applied and the results of our study.This highlights the benefit of performing a full integration of extreme offshore conditions together with their dependence in hydrodynamic simulations for screening out the coastal areas potentially exposed to flooding around Corsica.More generally, the method can be adapted to areas where only topobathymetric studies have been carried out to map potential flooding areas.
There are however several aspects that might nuance this finding.First, uncertainty in the multivariate extreme value analysis was only partly integrated.Indeed, although we have taken into account some of the uncertainties (CI given with marginals, correction on Hs, and diagnosis on thresholds in H&T04), a full propagation of all the sources of uncertainty has not been performed.Thus, a simple-but-efficient approach consisting in affecting a 0.25 m margin on the TWL results showed that even in this situation, the 2 m threshold was hardly reached in several places around the island.Second, we did not apply the method in some areas like some cliffs or areas where topobathymetric data where insufficient at the time of the study.This constitutes a line for future research.Furthermore, the uncertainties also lie in the variables studied and the choice of dominant variable.Indeed, other variables could be studied, such as precipitation, river discharges, or even storm duration the same way as Hs [54].A sensitivity analysis on these variables could be conducted on this area as in [55].We can also note that we defined an independent wave event as a 3-day window (see Section 3.1), but [56] note that 25% of extreme wave events are consecutive events that hit the same location in less than a day.This could be relevant, to be taken into consideration, especially for "2100" future conditions (the 25% rate could be higher).
There are other points that we still need to work on.Indeed, we used the hydrodynamic models in a static way; it means that the propagation from offshore to nearshore does not take into account the dynamics of the phenomenon (dynamics of the event, flow velocity, etc.).With our static approach, the dynamic of an event is not taken into account and therefore the volume of water is only limited by the capacity of the low lying zone to fill up to an altitude corresponding to the total water level at the shoreline.Even if we attached importance to the dependency between several variables leading to flooding, the main component to determine low-lying zones potentially exposed to coastal flooding remains the topography.It could be relevant to consider a dynamic approach, which is more computationally demanding, but more precise in terms of overflow, breach flooding, volume of water assessment, or flood phenomena coming from rivers.Regarding methodological developments, future work could also concentrate on the comparison with alternatives approaches, such as a response approach [21], potentially aided by metamodeling techniques [57].
Furthermore, we looked at "2100" future conditions and, following the regulator's recommendations [46], we applied a 0.6 m value by 2100 to account the SLR in the hy-drodynamic models.Even if this value was derived from the RCP 8.5 scenario "likely range", it is interesting to note that, regarding national SLR planning, France has a low amount of SLR used in planning compared to neighboring countries.Indeed, the amount of SLR used in planning in France is less than 1 m by 2100, whereas in Greece, it is about 1 m by 2100, and in Belgium, it is almost 2 m when considering 2100 high-end SLR [58].Thiéblemont et al. showed that under the RCP 8.5, the 0.6 m value could be possibly doubled or tripled along the France coastline when considering high-end estimates of different components of sea-level projections [59].Besides, the recent Special Report on the Ocean and Cryosphere in a Changing Climate (SROCC), then the Sixth Assessment Report AR6 emphasized that there is a substantial likelihood that SLR will be outside the likely range [60,61].We note that these reports also identify deep uncertainties related to rare or indirectly deduced phenomenon, as the marine ice cliff instabilities taken into account in some global climate models.
Finally, the purpose here was the assessment of the total water level at the shoreline in order to inform on the low-lying zones in Corsica potentially exposed to flooding.The next steps of the method were to determine areas that were geomorphologically homogenous based on orthophotos and fieldwork analyses, in order to complete the mapping [62].Specific information on the swash were also computed using the SWASH model in its 1D configuration on some beaches identified as potentially impacted by swash and mechanical shocks due to waves, in current and/or future conditions.

Figure 1 .
Figure 1.Location of Corsica Island and Sète (red symbol).(a) The map is adapted from Google Earth.(b) Wind provenance and crossed distance.The map is adapted from ROL Corse (Coastal Observatory Network of Corsica).

Figure 1 .
Figure 1.Location of Corsica Island and Sète (red symbol).(a) The map is adapted from Google Earth.(b) Wind provenance and crossed distance.The map is adapted from ROL Corse (Coastal Observatory Network of Corsica).

Figure 2 .
Figure 2. (a) Recap of our multivariate extreme value analysis.(b) Nested boxes used all around Corsica Island for the numerical simulations; 50 m resolution grids are in red and 10 m resolution boxes are in orange.

Figure 2 .
Figure 2. (a) Recap of our multivariate extreme value analysis.(b) Nested boxes used all around Corsica Island for the numerical simulations; 50 m resolution grids are in red and 10 m resolution boxes are in orange.

3. 3 .
Defining the Multivariate Scenarios (Steps 3-5) 3.3.1.Dependency Modeling and Monte Carlo Simulations J. Mar.Sci.Eng.2021, 9, x FOR PEER REVIEW 9 of 22variable is extreme (i.e., exceeding a defined threshold) with respect to the individual marginal distributions, as well as the dependency relationship between variables.

Figure 4 .
Figure 4. Illustration of Monte Carlo simulation for variables Hs, U, and SWL based on 30 common years of declustered data (black dots).Simulated data (10,000 years simulated) are in blue for the dataset NWW3 560/MARS 8.

Figure 4 .
Figure 4. Illustration of Monte Carlo simulation for variables Hs, U, and SWL based on 30 common years of declustered data (black dots).Simulated data (10,000 years simulated) are in blue for the dataset NWW3 560/MARS 8.

Figure 5 .
Figure 5. (a) Representation of the 100-year joint exceedance isocontours for the dataset NWW3 560/MARS 8.The color of the isocontours represent the SWL values from 0.2 to 0.4 m.(b) Simplified 3D representation of the 100-year joint exceedance isocontours for the dataset NWW3 560/MARS 8.The black dots represent the selected combinations given as input of the hydrodynamic models.

Figure 6 .
Figure 6.Relationship between Tp and Hs: declustered data (black dots), Monte Carlo simulation results (blue dots) for the dataset NWW3 560/MARS 8.The median is plotted in red.

Figure 6 .
Figure 6.Relationship between Tp and Hs: declustered data (black dots), Monte Carlo simulation results (blue dots) for the dataset NWW3 560/MARS 8.The median is plotted in red.

J
. Mar. Sci.Eng.2021, 9, x FOR PEER REVIEW 13 of 22 still water levels.For example, on the eastern part of the island the combination (Hs = 2.2 m; Tp = 6.8 s; U = 12.1 s; SWL = 0.4 m), represented by a blue star in the red box in Figure8, is identified as an extreme scenario to force the hydrodynamic models.

Figure 8 .
Figure 8.The 100-year joint exceedance isocontours for all of the boxes (in color) analyzed around Corsica, the colorbar corresponds to SWL (m).

Figure 8 .
Figure 8.The 100-year joint exceedance isocontours for all of the boxes (in color) analyzed around Corsica, the colorbar corresponds to SWL (m).

J 22 Figure 9 .
Figure 9.Total water level at the shoreline obtained: (a) for current conditions; (b) for future conditions (2100).

Figure 9 .
Figure 9.Total water level at the shoreline obtained: (a) for current conditions; (b) for future conditions (2100).

Figure 10 .
Figure 10.Scenarios leading to maximal water level at the shoreline in Gulf of Porto and Gulf of Girolata: (a) for directions Dp = 270° and Du = 240°; (b) for directions Dp = 240° and Du = 240° in current conditions.Total water level at the shoreline obtained in Gulf of Porto and Gulf of Girolata: (c) for directions Dp = 270° and Du = 240°; (d) for directions Dp = 240° and Du = 240° in current conditions.

Figure 10 .
Figure 10.Scenarios leading to maximal water level at the shoreline in Gulf of Porto and Gulf of Girolata: (a) for directions Dp = 270 • and Du = 240 • ; (b) for directions Dp = 240 • and Du = 240 • in current conditions.Total water level at the shoreline obtained in Gulf of Porto and Gulf of Girolata: (c) for directions Dp = 270 • and Du = 240 • ; (d) for directions Dp = 240 • and Du = 240 • in current conditions.

Table 1 .
Data used in the study, spatial resolution, temporal resolution, and sources.

Table 2 .
Combinations Hs, Tp, U, and SWL selected from the Gulf of Porto case.

Table 2 .
Combinations Hs, Tp, U, and SWL selected from the Gulf of Porto case.