*3.2. Marginal Distributions (Step 2)*

Figure 3 illustrates an example of the marginal probability distribution fitting *Fi*(*x*) of each variable *Xi* for (Hs; SWL; U) using the GPD defined as:

$$F\_{\vec{i}}(\mathbf{x}) = \begin{cases} \tilde{F}\_{\vec{i}}(\mathbf{x}) & \mathbf{x} \le u\_{\vec{i}} \\\ 1 - \left(1 - \tilde{F}\_{\vec{i}}(u\_{\vec{i}})\right) \left[1 + \xi\_{\vec{i}} \frac{\left(\mathbf{x} - u\_{\vec{i}}\right)}{\sigma\_{\vec{i}}}\right]\_{+}^{-\frac{1}{\xi\_{\vec{i}}}} & \mathbf{x} > u\_{\vec{i}} \end{cases} \tag{1}$$

where *ui* is a high threshold to choose, *σi* and *ξi* are the GPD parameters.

**Figure 3.** Example of the marginal distributions' fitting. (**a**) Location of the points considered. NWW3 560 provides U (m/s) and Hs (m), MARS 8 provides SWL (m). Generalized Pareto Distribution (GPD) and associated 70% (red dashed lines) and 95% (red thin lines) CI for (**b**) Hs, (**c**) SWL, and (**d**) U. Color bars give the directions Dp and Du (◦).

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]. Uncertainties 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.

### *3.3. Defining the Multivariate Scenarios (Steps 3–5)*

### 3.3.1. Dependency Modeling and Monte Carlo Simulations

First, the variables *Xi* for (Hs; SWL; U) are transformed onto the Gumbel scales using the standard probability integral transformation to give the variables *Yi*. If *Y*−*<sup>i</sup>* denotes a vector for all variables excluding *Yi*, a multivariate non-linear regression model is applied to *Y*−*i*|*Yi* > ν:

$$Y\_{-i} = aY\_i + Y\_i^b \mathcal{W} \quad \text{for } \mathcal{Y}\_i \succ \mathbf{v} \text{ and } \mathcal{Y}\_i \succ \mathcal{Y}\_{-i} \tag{2}$$

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 *Yi*). 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 variable 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.** 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.
