<figure><img src="image_1.png" /></figure>

HAL
open science

# Simulation of extreme functionals in meteoceanic data: Application to surge evolution over tidal cycles

Nathan Gorse, Olivier Roustant, Jérémy Rohmer, Déborah Idier

<figure><img src="image_2.png" /></figure>

## To cite this version:

Nathan Gorse, Olivier Roustant, Jérémy Rohmer, Déborah Idier. Simulation of extreme functionals in meteoceanic data: Application to surge evolution over tidal cycles. 2025. hal-05212699v2

HAL Id: hal-05212699

https://hal.science/hal-05212699v2

Preprint submitted on 19 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Simulation of extreme functionals in meteoceanic data: Application to surge evolution over tidal cycles

Nathan Gorse¹, Olivier Roustant¹, Jérémy Rohmer², and Déborah Idier²

¹UMR CNRS 5219, Institut de Mathématiques de Toulouse, INSA, Université de Toulouse, France

²BRGM, F-45060 Orléans, France

August 19, 2025

## Abstract

We investigate the influence of time-varying meteoceanic conditions on coastal flooding under the prism of rare events. Focusing on conditions observed over half tidal cycles, we observe that such data fall within the framework of functional extreme value theory, but violate standard assumptions due to temporal dependence and short-tailed behavior. To address this, we propose a two-stage methodology. First, we introduce an autoregressive model to eliminate temporal dependence between cycles. Second, considering the model residuals, we adapt existing techniques based on Pareto processes. This allows us to build a simulator of extreme scenarios, by applying inverse transformations. These simulations depend on an initial time series, which can be randomly selected to tune the desired level of extremes. We validate the simulator performance by comparing simulated times series with observations, through several criteria, based on principal component analysis, extreme value analysis, and classification algorithms. The approach is applied to the surge data, on the Gâvres site, located in southern Brittany, France.

## 1 Introduction

Coastal flooding frequently makes headlines due to the severe damage it causes and its long-lasting impacts. Notable examples include Storm Johanna in 2008 and Storm Xynthia in 2010, which affected vast areas and significantly raised public awareness of this risk. Coastal flooding is particularly concerning because coastal zones are densely populated and host high-value activities [Dawson et al., 2009]. This concern is expected to grow in the coming decades, as projections estimate that between 2.1 and 2.9 billion people will live in near-coastal areas by 2100 [Reimann et al., 2023]. Consequently, improving the forecasting of such events has become a critical objective for public authorities [Toimil et al., 2017].
Numerical hydrodynamic simulators provide essential insights into the complex relationship between coastal flooding and extreme meteoceanic conditions. In this study, we focus on extreme surge-induced flooding [Chaumillon et al., 2017] but with the added complexity that the governing conditions are time-varying and evolve over tidal cycles (3 hours around the tidal peak). In our context, the challenge lies in the fact that the inputs to the numerical models are time-dependent functions, i.e., functional inputs, and our goal is to generate extreme time series. To this end, we adopt a functional extreme value analysis framework which has been applied in particular to extreme windstorms, heavy spatial rainfall and temperature [de Fondeville and Davison, 2022, Ferreira and de Haan, 2014].

To simulate extreme data, we aim to construct a probabilistic model for these time series. However, applying the theory of functional extremes typically requires two main assumptions:

1. The observations must be independent and identically distributed (i.i.d.).

2. The time series should exhibit regular variation [Hult and Lindskog, 2005], implying heavy-tailed marginal distributions.

Unfortunately, these assumptions are hardly met in many situations of interest for coastal and ocean engineering; see for instance the study by Mackay et al. [2021] for the problem of dependence, and see that of Sando et al. [2024] for an example where the extreme distributions of waves do not exhibit regular variation. This is also the case for the surge time series analysed in our case in France (full details are provided below in Sect. 2), which motivated this study.

To overcome this, we propose to proceed in two steps:

1. Temporal dependence is addressed using an autoregressive model, which captures inter-cycle dependence while preserving intra-cycle structure.

2. To recover heavy-tailed marginals, we apply the approach of Opitz et al. [2021], combining a semi-parametric model of the empirical distribution with a Fréchet marginal transformation.

Following this preprocessing stage, the transformed time series satisfy the requirements for functional extreme modeling. We then build a probabilistic model following the framework in Dombry and Ribatet [2015], which relies on a polar coordinates representation [Kokoszka et al., 2019] whose components asymptotically follow a Pareto process.

This probabilistic model allows us to simulate extreme time series by first sampling from that Pareto process, then applying inverse transformations to recover time series in the original space. These simulations depend on an initial time series, which is required to invert the autoregressive model. The desired extremeness level can be tuned by the choice of this reference series, enabling a range of applications – from data-like simulations, following the same law as extreme observations, as sequences of consecutive extreme events.

A special care is paid to assess the quality of our simulated extremes by proposing a series of tools including Principal Component Analysis (PCA), extreme value analysis and
two-sample classification tests [Lopez-Paz and Oquab, 2016]. The overall performance of the method is evaluated using the storm surge case study.

This article is organized as follows. In Section 2, we describe the context of our case study and analyze the characteristics of the observations. Section 3 details the two-stage methodology leading to the probabilistic model on functional data. Section 4 explains how this model is employed to simulate extreme data. Section 5 applies the full framework to the storm surge case study and evaluates the simulation results.

# 2 A motivating case study

## 2.1 Data description

We study the site of Gâvres, located on the French Atlantic coast, which is particularly sensitive to coastal flooding. The town has been impacted several times in its history like in 1904 or in 1924. Despite the evolution of its coastal defenses, about 120 houses were flooded during storm Johanna (08/03/2008 – 10/03/2008) Idier et al. [2020]. Thus, it is a relevant study case to analyze extreme meteoceanic conditions. To this end, we use the hindcast database presented in Idier et al. [2020].

The small town lies in a macrotidal area such that coastal floods are strongly controlled by tides and occur around the high tide. Hence we concentrate on meteoceanic conditions over half tidal cycles, occurring within ±3 hours around high tide. Since we use a fixed time step of 10 minutes, we have several time series of length 37 for each tidal cycle. We consider the univariate case with an analysis of the surge (S) over the period (1979-2016) with approximately 27000 time series that will be what we call “observation” in what follows. Traditionally, analyses focus on time series that follow a consistent pattern. For instance, in Tendijck et al. [2024], all the time series considered are centered around a peak of significant wave height. On the contrary, the time series in our dataset do not exhibit a specific shape (Fig. 1). We will use the notation $X_M^t$ to describe the value obtained at time $t$ for the $M^{th}$ tidal cycle. Thus, $X_M^1, \dots, X_M^{37}$ represents the time series associated with the $M^{th}$ tidal cycle.
## 2.2 Two issues with the standard assumptions

<figure><img src="image_4.png" /><figcaption>Figure 1: Sampling of 20 observations from 1979 to 2016 for surge S (left) and for tide (right).</figcaption></figure>

We study the time series that distinguish themselves by an extreme behavior. This trait refers to the framework of extreme functionals, which is based on some assumptions discussed in the next section.

## 2.2 Two issues with the standard assumptions

**Correlated observations.** We generally assume that the observations (here time series of length 37) are i.i.d [Dombry and Ribatet, 2015], which implies that, for all $t = 1, \dots, 37$, the series $(X_1^t, \dots, X_N^t)$ with $N$ the number of observations should constitute an i.i.d sample. However, as illustrated by the autocorrelogram ACF and the partial autocorrelogram PACF at time $t = 1$ (Fig. 2), we observe that the time series are correlated for each value of $t$. The correlation level is around 0.10 when the lag in terms of tidal cycle index equals 50.

<figure><img src="image_5.png" /><figcaption>Figure 2: Time series analysis of $S$ at the first time step of each event ($t = 1$) with a lag in terms of tidal cycle index: (a) ACF; (b) PACF</figcaption></figure>

In addition to these strong correlations, we observe a seasonality of the variable with
## 2.2 Two issues with the standard assumptions

higher values observed during winter months (September-March) (see Appendix 6). In conclusion, the independence assumption is rejected.

**Departure from the regular variation hypothesis.** The simulation method based on regular variations requires that $X_M^t$ is heavy-tailed distributed for each value of $t$. To properly define this concept, we recall some prerequisites in extreme value theory. Let $U$ be a univariate random variable and assume that $U$ lies in the attraction domain of a so called Generalized Extreme Value (GEV) law [Balkema and De Haan, 1974, Ferreira and de Haan, 2014]. This means that, for a $n$-sample $U_1, \dots, U_n$ following the law of $U$, there exists constants $a_n > 0, b_n$ such that

$$M_n = \max_{i \in [[1, n]]} \frac{U_i - b_n}{a_n}$$

converges in distribution to a GEV law. In the Peak-Over-Threshold (POT) approach, under the same assumption, the exceedances of $U$ over a high threshold $u$ follow a Generalized Pareto distribution (GPD). Its expression is

$$\mathbb{P}(U \le x | U > u) = \begin{cases} 1 - \left(1 + \gamma \frac{x-u}{\sigma}\right)^{-\frac{1}{\gamma}} & \text{if } \gamma \ne 0 \\ 1 - \exp\left(-\frac{x-u}{\sigma}\right) & \text{otherwise} \end{cases} \quad (1)$$

where $x > u$ verifies $1 + \gamma \frac{x-u}{\sigma} > 0$. We denote this law GPD($u, \sigma, \gamma$) and $F_{GPD}$ the corresponding cumulative distribution function (cdf), where $\gamma, \sigma$ are respectively shape and scale parameters. We say that $U$ is heavy-tailed distributed if $\gamma > 0$.

Following the peak-over-threshold approach [Coles et al., 2001], the extremeness of a time series is quantified by a homogeneous functional $\ell$, which maps the time series to a scalar value. This function is referred to as a cost functional [Dombry and Ribatet, 2015] or a risk functional [de Fondeville and Davison, 2022]. A time series is classified as extreme if the corresponding value of $\ell$ exceeds a specified threshold. Using the $L^2$ norm as the function $\ell$, we analyze the extremes of $\ell(X_M)$.

As we estimate the evolution of $\gamma$ according to the number of exceedances, we analyze the sign of the shape parameter $\gamma$ (Fig. 3). The Hill estimator explored in De Haan and Resnick [1998] (dotted line) is based on a property of extreme observations under heavy-tail hypothesis ($\gamma > 0$). If the assumption holds, the estimator converges. Since the plot shows no stable region for the Hill estimator, we conclude that the condition $\gamma > 0$ is not met. Hence, the regular variations assumption is not satisfied.
<figure><img src="image_6.png" /><figcaption>Figure 3: Estimated value of the GPD shape parameter $\gamma$ with various methods (solid line: moments estimator, dotted line: Hill estimator, dashed line: MLE estimator)</figcaption></figure>

# 3 Construction of a probabilistic model

As discussed in the introduction, we aim at building a simulator accounting for temporal dependence and short-tailed behaviour of the observations. Our methodology has two stages. The first stage is pre-processing and aims at solving the two issues raised in the previous section. It is divided in two steps. First, we obtain independent data by considering the residuals of an autoregressive model (“whitening” step). Second, starting from these residuals, we recover heavy-tail distributions by using marginal transformations. In a second stage, we construct a probabilistic model based on a polar coordinates representation. The flowchart of the whole methodology is summarized in Fig. 4. We further detail each step in the following.

<figure><img src="image_7.png" /><figcaption>Figure 4: Construction of the polar coordinate representation from raw data (corresponding sections indicated near the boxes)</figcaption></figure>

## 3.1 Preprocessing stage

### 3.1.1 Recovering the independence case: Whitening with an auto-regressive model.

Starting from observed time series, we aim at removing the temporal dependence between successive time series. We first detrend the data by simply writing

$$X_M^t = \alpha_t M + c_t + \eta_M^t = \alpha_t M + \tilde{X}_M^t \quad (2)$$

where $\alpha_t$ represents the slope, $c_t$ is the intercept and $\eta_M^t$ is the residual term. The detrended time series is denoted by $\tilde{X}_M^t$.
## 3.2 Modeling stage based on a polar coordinate representation

Then, we impose a minimal duration denoted $\Delta$ between each event, which implies from now on that we only consider the observations with $M = 1 + \Delta m$ with $m = 1, 2, \dots$. Finally, for each value of $t$, we model the temporal dependence, with respect to $M$, by an autoregressive model, leading to

$$ \tilde{X}_M^t = \beta_0^t + \sum_{i=1}^p \beta_i^t \tilde{X}_{M-\Delta i}^t + \varepsilon_M^t \quad (3) $$

where $\beta_0^t$ is a constant term, $(\beta_i^t)_{i=1,\dots,p}$ are the model parameters and $\varepsilon_1^t, \varepsilon_{1+\Delta}^t, \dots, \varepsilon_N^t$ are the model residuals. These residuals form a white noise, i.e. are i.i.d. random variables, that we will now consider.

### 3.1.2 Recovering heavy-tailed distributions: Using a Fréchet marginal transformation.

Starting from the residuals $\varepsilon_M^t$ obtained in the previous step, we aim at obtaining a heavy-tailed distribution. A simple way to proceed is to use the transformation $\mathcal{T}: X \mapsto -1/\log(F_X(X))$ where $X$ is a continuous random variable with cdf $F_X$. Indeed, $\mathcal{T}(X)$ follows a Fréchet distribution (which is heavy tailed). However, in our case, the cdf of $\varepsilon_M^t$ must be estimated. Following Opitz et al. [2021], we approximate it by a mixture of the empirical cdf in the non-extreme part and a GPD distribution in the extreme region. Denoting $F^{\text{emp}}$ the empirical cdf and $F_{\text{GPD}}$ the distribution function of a GPD variable (1), $\hat{F}_{\varepsilon^t}$ has the form

$$ \hat{F}_{\varepsilon^t}(\varepsilon) = \begin{cases} F_t^{\text{emp}}(\varepsilon) & \text{if } \varepsilon < u^t \\ (1 - p_u) + p_u F_{\text{GPD}}(u^t, \sigma^t, \gamma^t)(\varepsilon) & \text{otherwise} \end{cases} \quad (4) $$

where $p_u$ is a parameter and $u^t$ verifies $p_u = 1 - \hat{F}_{\varepsilon^t}(u^t)$. In the sequel, we will denote $Z_M^t = \mathcal{T}(\varepsilon_M^t)$.

## 3.2 Modeling stage based on a polar coordinate representation

Since the pre-processing stage enables us to get back to the regular variations assumptions, we can now model the law of extreme time series. We use the framework introduced by Dombry and Ribatet [2015], in which functional data are modeled using a polar decomposition of the form

$$ f = \ell(f) \times A(f). \quad (5) $$

where $f = Z_M$ in our case. The radius $\ell(f)$, chosen here as the $L^2$ norm of $f$, allows to define extreme functionals with a scalar threshold $u_\ell$. Under the assumption of heavy-tailed distributions, the components $(\ell(f), A(f))$ become asymptotically independent, conditional on $\ell(f) > u_\ell$ as $u_\ell \to \infty$, and the conditional process $f/u_\ell$ $| \ell(f) > u_\ell$ converges in distribution to a $\ell$-Pareto process. As a consequence, we have

$$ \ell(f) \perp A(f) | \ell(f) > u_\ell, \quad u_\ell \to +\infty \quad (6) $$
Thus, we must select a convenient threshold $u_\ell$ as defined in Eq. (6). One selection method is based on the convergence of the angular component $\Theta_M = f/\ell(f)$ of the $\ell$ Pareto process. The main idea is to analyze the behavior of $(\Theta_M | \ell(Z) > u_\ell)$ as $u_\ell$ increases gradually [Clémençon et al., 2024].

After choosing a sufficiently high threshold $u_\ell$, we obtain the polar coordinate representation $(\Theta_M, \ell(Z_M) | \ell(Z_M) > u_\ell)$, which is used to simulate new extreme time series.

# 4 Simulation of extreme time series

We propose a two-step simulator of new extreme observations based on the probabilistic model. In the first step (Sect. 4.1), we model the law of the polar coordinates and use properties of Pareto processes to generate simulations of $Z_M$, denoted $Z_{sim}$. In the second step (Sect. 4.2), we apply reverse transformations to $Z_{sim}$ to reconstruct time series of interest. We describe the complete method in Fig. 5, where the forward-pointing arrows represent the transformation step introduced in Sect. 3 whereas the backward-pointing arrows indicate the simulation phase. As previous sections describe the construction of the model, we now focus on the simulator derived from this representation.

<figure><img src="image_8.png" /><figcaption>Figure 5: Two-stage methodology to simulate extreme time series (corresponding sections indicated near the boxes). The simulation starts from the left with the polar coordinate representation resulting from the forward sense (described in Sect. 3)</figcaption></figure>

## 4.1 Modeling of polar coordinates

**Functional analysis of the angular component.** The simulation of extreme time series is based on the polar decomposition given in Eq. (5). Simulating the radial component $\ell(f)$ is straightforward as its asymptotic distribution is explicit (a Pareto distribution). In contrast, simulating the angular component $A(f)$ is much more challenging.

Several simulation methods have been proposed in the literature [de Fondeville and Davison, 2018], which rely on a spectral representation of the coordinate-wise maximum of $Z_M^\ell$ using Gaussian processes [Dombry and Ribatet, 2015]. They all depend on the choice of a parametric model for the variogram. In contrast, we use a dimensionality reduction method, here PCA, before fitting parametric models. We follow Clémençon et al. [2024] by decomposing
## 4.2 Applying the reverse transformations

$\Theta_M$ in a orthonormal basis with

$$ \Theta_M^t = \mu^t + \sum_{j=1}^{37} \langle \Theta_M - \mu, \nu_j \rangle \nu_j^t = \mu^t + \sum_{j=1}^{37} C_j \nu_j^t \quad (7) $$

where $(\nu_j)_{j \ge 1}$ are eigenvectors sorted by decreasing order of the eigenvalues $(\lambda_j)_{j \ge 1}$, $(C_j)_{j \ge 1}$ the PCA coordinates and $\mu^t = \mathbb{E}(\Theta_M^t)$.

To reduce the dimensionality of the problem, we retain only the first $J$ eigenvectors and denote by $R_J$ the ratio of total inertia explained by the truncated basis. We truncate the PCA basis when adding another eigenvector does not significantly decrease $(1 - R_J)$. This dimensionality reduction step allows us to estimate only the joint distribution of the scores, enabling the simulation of new extreme observations.

**Modeling the law of the angular component.** We account for the dependence between PCA coordinates by using the copula [Rüschendorf, 2009] C defined as

$$ F(C_1, \dots, C_J) = C(F_1(C_1), \dots, F_J(C_J)) = C(v_1, \dots, v_J) \quad (8) $$

where $F$ is the joint distribution function of the vector $(C_1, \dots, C_J)$ and $F_i$ is the distribution function of the $i^{th}$ marginal.

We approach the joint law of $(v_1, \dots, v_J)$ by fitting a parametric copula model. As copula models are defined for bivariate settings, we use vine copulas [Czado and Nagler, 2022] when $J > 2$, which use a tree representation of the dependence structure where each edge corresponds to a bivariate copula. Then, the law of the whole vector is modeled by combining the pairwise models in a hierarchical manner. To simulate new scores $(C_{sim,1}, \dots, C_{sim,J})$, we apply $F_i^{-1}$ for $i = 1, \dots, J$ to new values $(\tilde{v}_1, \dots, \tilde{v}_J)$ sampled from the model.

With this method, we generate $n_{sim}$ new angles $\Theta_{sim}$ by rewriting Eq. (7), which gives:

$$ \Theta_{sim}^t = \mu^t + \sum_{j=1}^{J} C_{sim,j} \nu_j^t \quad (9) $$

As $\ell(\Theta_M)$ must equal 1, we return $\Theta_{sim} = \Theta_{sim}/\ell(\Theta_{sim})$. This simulation step makes it possible to simulate new extreme time series, which is discussed in the next section.

## 4.2 Applying the reverse transformations

**Over a threshold decomposition and marginal transformation.** Thanks to the independence established in Eq. (6), we generate new extreme time series $Z_{sim}$ by simulating new values of the radius $\ell(f)$ and by combining them with the simulations $\Theta_{sim}$. As $Z_M^t \sim \text{Frechet}(1)$, we impose that $Z_{sim}^t > 0$ for every value of $t$. To this aim, we apply a rejection sampling: if a curve is strictly above 0, we accept it; otherwise, we discard it and we simulate both a new $\Theta_{sim}$ and a new radius.

The final step consists in applying the reversed transformation $T^{-1}$ to the accepted time series, yielding extreme time series $\varepsilon_{sim}$, which can be used to simulate new realizations of the time series $\tilde{X}_M$.
## 4.3 Implementation.

**Inverting the autoregressive model.** In the case where $p = 1$, we obtain simulated time series for $(\tilde{X}_{\text{sim}}^t)_{t=1,...,37}$, depending on initial time series $(\tilde{X}_{M-\Delta}^t)_{t=1,...,37}$ with

$$ \tilde{X}_{\text{sim}}^t = \beta_0^t + \beta_1^t \tilde{X}_{M-\Delta}^t + \varepsilon_{\text{sim}}^t \quad (10) $$

The values and the shape of $\tilde{X}_{\text{sim}}$ depend on the levels reached by $\tilde{X}_{M-\Delta}$. In this context, we can use an extreme $\tilde{X}_{M-\Delta}$ to produce consecutive extremes. As our objective is to simulate data-like extreme time series, we must choose an appropriate time series $\tilde{X}_{M-\Delta}$ given $\varepsilon_{\text{sim}}$. Thus, we can adopt a conditional or an unconditional approach to sample $X_{M-\Delta}$. If we choose a conditional sampling, we account for the joint distribution of the couple $(X_{M-\Delta}^t, \varepsilon_M^t)$ as these variables may still exhibit dependence, which should be preserved during the sampling process. A simple approach to incorporate this dependence is to sample from the empirical distribution of $(\ell(\tilde{X}_{M-\Delta}) | \ell(\varepsilon_M) = x)$. For each value $x$ of $\ell(\varepsilon_M)$, we define a symmetric window centered at $x$ that include 20 neighboring data points. Then, we draw randomly with replacement one point from this window.

Finally, following Eq. (2), we incorporate the trend of $X_M^t$ over the period by using the equation:

$$ X_{\text{sim},M}^t = \alpha_t M + \tilde{X}_{\text{sim}}^t \quad (11) $$

where $M$ is known thanks to the sampling of the $\tilde{X}_{M-\Delta}$. Another approach would be to predict today's forcing conditions by using the current index. Yet, we focus on the detrended time series to facilitate direct comparisons between our simulations and the observed extreme time series. We apply our method in the next section to the case study and discuss the validation of the simulator.

## 4.3 Implementation.

The method is coded in [R Core Team, 2023] and we used several packages. POT [Ribatet and Dutang, 2022] and extRemes [Gilleland and Katz, 2016] enabled us to study exceedances at each time. We based the marginal transformation code on [Opitz et al., 2021]. Besides, we used FactoMineR [Lê et al., 2008] for PCA and VineCopula [Nagler et al., 2023] to approach the law of copulas. We also use the package [Genz and Bretz, 2009]. Finally, we used e1071 [Meyer et al., 2023] and randomForest [Liaw and Wiener, 2002] to construct support-vector machines and random forest classifiers.

# 5 Application to the case study

In this section, we apply the methodology detailed in Sect. 3 and Sect. 4 to the available database from the site of Gâvres. Using the reading directions of Fig. 5, we present some implementation details and results for the key steps of the procedure.
## 5.1 Forward sense: Construction of the polar coordinate representation

### 5.1 Forward sense: Construction of the polar coordinate representation

Following Fig. 4, we apply a two-step method to obtain the representation used for generating new extreme time series. Since our data do not respect the framework assumptions (Sect. 2.2), the first step consists in retrieving the framework setting whereas the second step enables to produce the polar representation.

**Pre-processing stage: Whitening and marginal transformation.** First, by applying the whitening stage described in Sect. 3.1.1, we detrend the time series $X_M^t$ with Eq. (2) and consider the detrended time series $\tilde{X}_M^t$. Since significant correlations remain, we retain only one observation out of 3, corresponding to $\Delta = 3$. The following graphic (Fig. 6) illustrates the database construction, where the grey rectangles represent the tidal cycles that are excluded from the analysis.

<figure><img src="image_9.png" /><figcaption>Figure 6: Summary of the database construction: selection of winter time series and one cycle out of 3</figcaption></figure>

We now analyze the pair of correlated variables $(\tilde{X}_M^t, \tilde{X}_{M+\Delta h}^t)$. The slow decay of correlation levels on Fig. 2a and the significant peak in Fig. 2b at $h = 1$ suggest the use of the models Eq. (3) with order $p = 1$. Thus, we apply for each value of $t$ an autoregressive model of order 1 denoted AR(1), which reduces effectively the correlation value. We describe with more details the effect of these models in Appendix 7.1.

Besides, since we study $\varepsilon_M^t$'s exceedances for a given $t$, we assess extremal dependence as defined in Coles et al. [1999]. To this end, we apply the statistical test proposed by Thomas [2001] whose null hypothesis is the asymptotic dependence of the residuals $\varepsilon_M^t$. The null hypothesis is rejected at confidence level 1% ($p$ value < $e-10$), indicating that the residuals $\varepsilon_M^t$ are asymptotically independent.

Then, as detailed in Sect. 3.1.2, we aim at obtaining heavy-tailed marginals for each value of $t$. To achieve this, following Opitz et al. [2021], we estimate the cdf of $\varepsilon_M^t$ and apply the transformation $\mathcal{T}$ to the residuals, obtaining Frechet marginals $Z_M^t$ for each value of $t$.

An intermediate step involves selecting the threshold parameter $u^t$ used in Eq. (4). To guide this choice, we use the property that for any ($w^t > u^t$) the exceedances ($\varepsilon_M^t > w^t$) should
## 5.1 Forward sense: Construction of the polar coordinate representation

follow a GPD law if $u^t$ is a suitable threshold. Accordingly, we define updated parameters $(\sigma', \gamma')$ that remain stable as $w^t$ increases. With the same objective, we analyze, for each value of $t$, the behaviour of the mean residual life (MRL) $\mathbb{E}(\varepsilon_M^t - u^t | \varepsilon_M^t > u^t)$ [Coles et al., 2001] as a function of $u^t$. Guided by these tools, we define $u^t$ as the 90% percentile of $\varepsilon_M^t$ since the MRL plot becomes linear beyond this point and the couple $(\sigma', \gamma')$ remains stable for all thresholds ($w^t > u^t$) (see Appendix 7.2).

To assess the quality of the marginal transformation, we compare empirical quantiles $Z_M^t$ with those of the Frechet distribution across various values of $t$. The levels obtained are close to the theoretical ones, demonstrating the consistency of the transformation.

The transformation induces a change in the tail behavior of the distribution of $\ell(\cdot)$: while $\ell(\varepsilon_M)$ is light-tailed distributed (see Fig. 7a), whereas $\ell(Z_M)$ displays heavy-tailed behavior (Fig. 7b), as confirmed by the convergence of the Hill estimator.

<figure><img src="image_10.png" /><figcaption>Figure 7: Estimation of $\gamma$ for: (a) $\ell(\varepsilon_M)$; (b) $\ell(\mathcal{T}(\varepsilon_M))$ (solid line: moments estimator, dotted line: Hill estimator, dashed line: MLE estimator)</figcaption></figure>

**Modeling stage based on a polar coordinate representation.** As we get back to the framework assumptions with the pre-processing step, we now choose the threshold $u_\ell$ that defines the set of extreme observations via the condition $\ell(Z_M) > u_\ell$ of Eq. (6). Following Sect. 3.2, the choice of $u_\ell$ depends on the convergence behavior of $\Theta_M$ when $u_\ell$ goes to infinity. In our case, $\Theta_M$ appears to converge in distribution when the number of extreme observations is between 200 and 300, corresponding to approximately 4% and 6% of the database (see Appendix 8.1).

As a consequence, we select the top 5% of the dataset, corresponding to approximately 7 extreme observations per year, and represent a sample of the resulting $n = 259$ time series on their original scale and on the Frechet scale (Fig. 8). The figure shows that applying the marginal transformation amplifies the variations between successive measures.
## 5.2 Backward sense: Simulation of extreme time series

<figure><img src="image_11.png" /><figcaption>Figure 8: Sample of 20 time series selected for the surge (left panel: original scale, right panel: Frechet scale with 95% confidence intervals in dotted lines)</figcaption></figure>

## 5.2 Backward sense: Simulation of extreme time series

The transformations of the previous section enable us to use a convenient probabilistic model for extreme time series. Following Fig. 5, we now describe the simulator based on this model, which uses the reverse transformations.

**Modeling the law of polar coordinates.** We use the polar decomposition defined in Eq. (5) as a basis for our simulation method. After selecting the $n = 259$ extreme observations, we first compute their angular component $\Theta_M$ and derive the corresponding PCA basis Eq. (7). Following Sect. 4.1, we select the first $J = 3$ eigenvectors of the basis as the ratio $(1 - R_J)$ slowly decreases beyond the third dimension (see Appendix C: Fig. 21). With this choice, the selected dimensions explains $R_J = 83\%$ of the variance.

We use vine copulas and model the joint distribution of the coordinate vector by fitting a range of copula families, namely elliptical and Archimedean copulas such as Tawn and Clayton models, to the bivariate vectors defined by the vine structure. Finally, for each pair of coordinates, we use the Akaike Information Criterion (AIC) [Akaike, 1974] to select the best copula and estimate the tail coefficients $(\lambda_+, \lambda_-)$ along with the Kendall's $\tau$. In our case, a $t$ copula is selected for each pair of coordinates. Please refer to Appendix 8.3 for additional information about the modeling.

The iso-density contours of the fitted law are consistent with the distribution of the data points. To go further, as we assume that the observations $(v_1, \dots, v_J)$ follow the fitted copula model, we perform some goodness-of-fit tests presented in Schepsmeier [2019]. We do not reject the null hypothesis when we use the bootstrap test based on the White statistic [White, 1982] ($p$ value $\approx 0.5$).

Using the composite copula model, we simulate new vectors $(v_{sim,1}, v_{sim,2}, v_{sim,3})$ and transform each coordinate back to its initial scale by applying $F_i^{-1}$ for $i = 1, \dots, 3$. We simulate $n_{sim} = 2,000$ triplets of coordinates from the model and we compare the simulated values with the coordinates of the extreme observations (Fig. 9). The marginals obtained in the
## 5.2 Backward sense: Simulation of extreme time series

simulations are consistent with the data distributions, indicating that our model is consistent with the dependence structure of the coordinate vector.

<figure><img src="image_12.png" /><figcaption>Figure 9: Coordinates obtained of 2,000 simulated and recorded extreme time series: (a) $(C_1, C_2)$; (b) $(C_2, C_3)$; (c) $(C_1, C_3)$</figcaption></figure>

We use the coordinates $C_{sim}$ to obtain the angular components $\Theta_{sim}$ by applying Eq. (11). The resulting patterns resemble that of the extreme observations although the simulated curves appear to be smoother.

**Simulating new extreme times series: Applying reverse transformations**

Following Sect. 4.2, we generate extreme time series $Z_{sim}$ by combining the simulated angles $\Theta_{sim}$ with simulations of $\ell(Z)$. Then, we apply the inverse transformation $T^{-1}$ to $Z_{sim}$ to return to the original scale, which produces simulated extreme residuals $\varepsilon_{sim}$.

Using Sect. 4.2, we now invert the estimated AR(1) models to return simulated time series $\tilde{X}_{sim}$. We consider two choices of $\tilde{X}_{M-\Delta}$: one observed during the extreme storm Johanna (10/03/2008) and another with a non-extreme behavior (27/10/1990), verifying $\ell(\mathcal{T}(\varepsilon_M)) < u_\ell$. The results, displayed in Fig. 10, show that the choice of $\tilde{X}_{M-\Delta}$ significantly affects both the values and the shape of $\tilde{X}_{sim}$. The left panel reveals notably high values, illustrating the effect of having consecutive extremes.
## 5.3 Consistency of simulations with extreme observed time series

<figure><img src="image_13.png" /><figcaption>Figure 10: Simulated time series obtained with two different types of $\tilde{X}_{M-\Delta}$ (dotted line) selected among extreme observations (left panel) or non-extreme ones (right panel).</figcaption></figure>

Since the variables $(\tilde{X}_{M-\Delta}^t, \varepsilon_M^t)$ are dependent, we apply the conditional sampling illustrated in Sect. 4.2 (see Appendix 8.4). Then, we generate new time series $\tilde{X}_{sim}$ by using $\varepsilon_{sim}$ and the selected $\tilde{X}_{M-\Delta}$. We see that the shape and the values of simulated and observed extreme time series are quite similar, suggesting that our simulations are consistent with the observations. However, we can use other methods to validate our simulator, which is discussed in details in the following.

<figure><img src="image_14.png" /><figcaption>Figure 11: Comparing simulated time series and extreme observations. Left panel: simulated sample of 100 extreme time series, right panel: sample of 100 extreme observations. The dotted lines represent the 95% confidence bands</figcaption></figure>

## 5.3 Consistency of simulations with extreme observed time series

We validate our simulator by comparing simulated time series with extreme observations. To this end, we use a set of diagnostic tools, using extreme value analysis, PCA and classification algorithms.
## 5.3 Consistency of simulations with extreme observed time series

### 5.3.1 Comparison of percentile levels.

First, we compare the empirical percentiles at each time point $t$ with those derived from the simulated time series (Fig. 12). To quantify uncertainty in the observations, we construct a bootstrap 95% confidence band by resampling the extreme time series 500 times. The simulated levels obtained are coherent with the observations as the simulations' percentiles lie within the data confidence band.

<figure><img src="image_15.png" /><figcaption>Figure 12: Percentiles obtained in the data and in the simulations (blue lines: data, orange dotted lines: simulations)</figcaption></figure>

### 5.3.2 Comparison of coordinates in a PCA basis.

The simulated extreme time series should exhibit shape similarities with the observed extremes. One simple way to compare the shape of the time series is to use dimensionality reduction methods and to compare the resulting coordinates. To this end, we analyse the angles of the time series as the simulations extrapolate the $L^2$ norm of the observations. We apply PCA to the normalized series $\tilde{X}/l(\tilde{X})$ and, to preserve potential differences, both simulated and observed series are further normalized by using the mean and the standard deviation of the observations.

Then, we compare the PCA coordinates of recorded time series with those of simulated time series (Fig. 13). The similarity in the score distributions suggests that simulated time series and extreme observations exhibit comparable behavior. Yet, this similarity is more pronounced in the first dimension, which accounts for the largest proportion of variability (here 60%). A KS test confirms that simulations' and observations' coordinates follow the same law in the first dimension ($p$ value $\approx 15\%$) but do not in the second dimension ($p$ value $\approx 0.8\%$).
5.3 Consistency of simulations with extreme observed time series

<figure><img src="image_16.png" /><figcaption>Figure 13: PCA coordinates for S in simulations (orange triangles) and extreme observations (blue dots)</figcaption></figure>

5.3.3 Comparison of the distribution upper tails.

Then, we compare the temporal dependence within each tidal cycle between observations and
simulations, with a focus on extremal dependence. Following [Opitz et al., 2021, Dell'Oro and
Gaetan, 2025], we assume that $\tilde{X}_M^t$ and $\tilde{X}_M^s$ are asymptotically dependent for all pairs $(t, s)$.
To quantify the dependence, we estimate the $\ell$-extremogram [de Fondeville and Davison,
2022] for both the recorded and simulated time series. We detail the approach in Appendix
9.1.

We focus on the values exceeding the 90% percentile and use 500 resamples to construct a
bootstrap 95% confidence band. The extremogram of the simulated time series lies within
the confidence band (Fig. 14) based on the observed extreme time series.

<figure><img src="image_17.png" /><figcaption>Figure 14: Extremogram for S for simulations (orange) and observations (blue) with confidence bands in dotted lines.</figcaption></figure>

Moreover, the density obtained at each time *t* should be an extrapolation of the observed
data. Directly using the quantiles of the simulations may produce very large values. This
## 5.3 Consistency of simulations with extreme observed time series

occurs because we simulate $\tilde{X}_{sim}$ verifying $\ell(T(\varepsilon_M)) > u_\ell$, meaning their extremes can be larger than those observed.

In this setting, we compare our values with the empirical quantiles of the observations and the fitted law of exceedances. Following Solari et al. [2017], the results are expressed in terms of return period $P$ and corresponding return levels, with approximately 7 extreme observations per year (Fig. 15). All implementation details are provided in Appendix 9.2.

The highest simulated levels fall within the confidence band, indicating that our simulations look similar to the extreme observations. Thus, they can be considered a reliable extrapolation of the observed data toward higher values.

<figure><img src="image_18.png" /><figcaption>Figure 15: Return level for simulations (orange triangles) and observations (blue dots), for two different values of $t$. Left panel: $t = 19$ (tidal peak), right panel: $t = 13$ (one hour before). The dotted lines represent 95% confidence bands.</figcaption></figure>

### 5.3.4 Comparison of predictions with classification algorithms.

Finally, following the same idea as two-sample classification tests [Lopez-Paz and Oquab, 2016], we rely on a classification-based approach to determine whether the generated extreme times series differ from the original dataset. We define a binary classification problem, where the target variable $y$ equals 1 if the input $x$ is a simulated time series, 0 otherwise. The goal is to learn a function $g$ that, given the input $x$, predicts the correct class $y$. If the simulations behave like the extreme observations, such classifier should hardly identify the simulated time series, resulting in an accuracy close to 50%.

To address this question, we begin by randomly drawing $n = 259$ simulated time series without replacement (Step 1) since the number of simulated time series $n_{sim} = 2,000$ differs from $n$. We construct a dataset composed of these sampled time series and extreme observations, ensuring that simulated time series make up 50% of the data. Then, after splitting this dataset into training and testing subsets, we train our model on the training part (Step 2). Finally, we evaluate the model's overall accuracy on the test set (Step 3).

To account for the influence of the sampling used in Step 1, we repeat 100 times these three steps and provide a confidence interval of order 90% for the overall accuracy. We use three classifiers with this approach: a support-vector machine (SVM) with a radial kernel, a generalized linear model (GLM), and a random forest [Hastie et al., 2009].

We apply this procedure with different types of input features: the time series, its norm and
## 5.3 Consistency of simulations with extreme observed time series

its angle. If we provide the time series, all classifiers achieve an overall accuracy close to 50% (see Table 1), meaning that the ML models struggle to identify simulated time series. This suggests that our simulations are quite consistent with the observed data.

For example, the distributions of $\ell(\tilde{X}_{sim})$ and $\ell(\tilde{X}_{M})$ are quite similar (see Appendix 9.4), with differences appearing for the largest values as the simulations better explore the extremes of $\ell(S)$. In contrast, when using the angle as input feature, the models are better able to distinguish between simulations and observations but the accuracy still remains low on the order of 60% in average.

<table><caption>Table 1: Confidence intervals at order 90% for the global accuracy</caption><thead><tr><th>Input</th><th>SVM</th><th>GLM</th><th>RForest</th></tr></thead><tbody><tr><td>Hyper-params</td><td>kernel: radial</td><td>link: logit</td><td>Nb-trees: 500</td></tr><tr><td>X</td><td>44-56</td><td>40-53</td><td>48-61</td></tr><tr><td>ℓ(X)</td><td>43-55</td><td>43-54</td><td>43-56</td></tr><tr><td>A(X)</td><td>49-64</td><td>40-53</td><td>59-72</td></tr></tbody></table>

This approach also reveals the influence of $\tilde{X}_{M-\Delta}$ as the attained accuracy depends on the sampling method, which is discussed in Appendix 10. The classifiers identify more easily the simulated time series when an unconditional sampling is used. For instance, the overall accuracy exceeds 50% when the time series itself is used as input.

## Conclusions

Motivated by the modeling of surge-induced coastal flooding, our goal is to develop a simulator for extreme time series. To this end, we align our time series with the framework of regular variations. As the observations do not meet standard assumptions, we get back to the framework setting by applying an autoregressive model and marginal transformation to the time series. It allows us to obtain a probabilistic model and to develop a simulation method capable of simulating data-like extreme time series, following the same law as extreme observations, but also consecutive extremes.

We apply this method to a hindcast database Idier et al. [2020] from the site of Gâvres in Brittany, which is considered as a set of observations in our developments. To evaluate the quality of our simulations, we provide several approaches, including a classification task where the goal is to distinguish between observed and simulated time series. Our results support the validity of our simulation approach as machine learning algorithms struggle to reliably distinguish the two groups.

This paper focuses on simulating univariate extreme time series. Yet, real-world forcing conditions evolve in a multivariate context [Rohmer et al., 2022, Gouldby et al., 2014] as other variables such as wind speed also intervene at each time step. Therefore, extending the approach to the multivariate case [Kim and Kokoszka, 2024] represents a promising direction for effectively modeling extreme forcing conditions.
Acknowledgement

Our work has benefitted from the AI Interdisciplinary Institute ANITI. ANITI is funded by
the France 2030 program under the Grant agreement n°ANR-23-IACL-0002.
# 6 Complements on data analysis

We present here the evolution of $S$ at tidal peak, $t = 19$, according to the month considered (Fig. 16). The extremes of winter values are much higher than the extremes of summer values, which is explained by a higher frequency of storms in winter.

<figure><img src="image_19.png" /><figcaption>Figure 16: Seasonal evolution of the distribution of $S$ at tidal peak</figcaption></figure>

# 7 Whitening and marginal transformation

## 7.1 Effect of the autoregressive models

We present here the effect of applying the autoregressive models on $\tilde{X}_M^t$. We can first use the ACF and the PACF on the residuals $\varepsilon_M$ (Fig. 17).
## 7.1 Effect of the autoregressive models

<figure><img src="image_20.png" /><figcaption>Figure 17: Pearson's correlations observed for S at $t = 19$. Left panel: ACF, right panel: PACF</figcaption></figure>

**Correlations of extremes.** We recall the definition of the $\chi$ measure [Coles et al., 2001]. Let $U$ and $V$ be two variables with uniform distributions and we note $\chi(u) = P(V > u | U > u)$ with $(0 < u < 1)$. They are asymptotically independent if $\chi(u) \to 0$ when $u \to 1$. In this case, $\bar{\chi}$ measures the dependence between the two variables.

In our case, the pair $(U, V)$ corresponds for a given $t$ and $h$ to the pair $(\varepsilon_M^t, \varepsilon_{M+\Delta h}^t)$ after a rank transformation. At the tidal peak ($t = 19$), 0 belong to the confidence interval of $\chi(u)$ when $u > 0.9$ for $h = 1$ (Fig. 18). As a consequence, we cannot reject the hypothesis of asymptotic independence.

<figure><img src="image_21.png" /><figcaption>Figure 18: Measures of asymptotic dependence for the surge for ($t = 19, h = 1$) with confidence band (grey lines). The blue lines represent the theoretical bounds of the correlation coefficient.</figcaption></figure>
## 7.2 Marginal transformation

We describe the evolution of diagnostic statistics when $u^t$ is between the median of $\epsilon_M^t$ and its 98% quantile. We observe the evolution of the MRL and the updated parameter (Fig. 19). We represent with the red line the threshold chosen for the marginal transformation. We do have a stable couple $(\sigma', \gamma')$ (top left and top right) and the MRL (lower left corner) seems to be a linear function.

<figure><img src="image_22.png" /><figcaption>Figure 19: Graphic diagnostics for the residual of $S$ at tidal peak ($t = 19$)</figcaption></figure>

We show the evolution of the dispersion index in the lower right corner. This tool is based on the link between extreme values and Poisson point process. $u^t$ is a consistent threshold if the index is around 1, which is the case for the threshold we chose.

# 8 Complements on the simulation of extreme time series

## 8.1 Choice of the extreme individuals

The convergence of $\Theta_M$ is close in spirit to the stability of GPD parameters (Sect.5.1) as we look for a convergence period of $E(|\langle\Theta_M, h_j\rangle|)$. We use here as $h_j(t) = \sin(2\pi j t)$ with $j \in \{1, \dots, 8\}$. We consider the evolution of the mean of $|\langle f/\ell(f), h_j\rangle| \mid \ell(f) > u_\ell$ when we increase $u_\ell$. Figure 20 shows that there is a period of stability for the mean projection for several functions when we have between 200 and 300 extreme time series.
## 8.2 Choice of the number $J$

<figure><img src="image_23.png" /><figcaption>Figure 20: Evolution of the mean projection for $T(S)$</figcaption></figure>

## 8.2 Choice of the number $J$

As we increase the number $J$ used in the truncated expression of $\Theta_M$, we analyze the evolution $1 - R_J$ (Fig. 21).

<figure><img src="image_24.png" /><figcaption>Figure 21: Evolution of the unexplained inertia of $\Theta_M$ for $S$</figcaption></figure>

After choosing $J = 3$, we simulate $n_{sim} = 2,000$ angles and compare our simulations with the angles of the extreme observations (Fig. 22).
<figure><img src="image_25.png" /><figcaption>Figure 22: Comparison of the angles $\Theta_{\text{sim}}$ and $\Theta_M$ for S</figcaption></figure>

### 8.3 Details on the copula family

We describe here the parameters estimated for the copula (Fig. 2).

<table><caption>Table 2: Modeling of $A(f)$'s coordinates with copulas: the first row explains the modeling of the law of the first two coordinates while the last one presents the modeling of the conditional law of $(C_1, C_3)$ knowing $C_2$</caption><thead><tr><th>Tree</th><th>Couples-cond</th><th>Name</th><th>par</th><th>par2</th><th>τ</th><th>λ<sub>+</sub></th><th>λ<sub>-</sub></th></tr></thead><tbody><tr><td>1</td><td>1,2</td><td>t</td><td>-0.2</td><td>2</td><td>-0.2</td><td>0.1</td><td>0.1</td></tr><tr><td>1</td><td>2,3</td><td>t</td><td>0.3</td><td>2</td><td>0.2</td><td>0.3</td><td>0.3</td></tr><tr><td>2</td><td>1,3;2</td><td>t</td><td>-0.2</td><td>3.4</td><td>-0.1</td><td>0.1</td><td>0.1</td></tr></tbody></table>

We compare the fitted copula models with observed coordinates by looking at the relationship between the two dimensions. In Fig. 23 and Fig. 24, we compare the data observations with the iso-density curves of the model respectively in the uniform scale and in the data scale.

Figure 23 shows that we are able to capture the most of the dependence structure as the repartition of data points is coherent with the levels found. However, if we return to the original scale of the coordinates, we see on Fig. 24 the limits of our modeling as some data points do not coincide with the iso-density curves of our model.
## 8.4 Choice of $X_{M-\Delta}$

<figure><img src="image_26.png" /><figcaption>Figure 23: Iso-density curves (dark lines) of the copula with the data coordinates (blue points): (a) ($v_1, v_2$), (b) ($v_2, v_3$)</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>Figure 24: Iso-density curves (dark lines) of the joint law with the data coordinates (blue points): (a) ($C_1, C_2$), (b) ($C_2, C_3$)</figcaption></figure>

We also use statistical test to compare our model with the law of the observations. The null hypothesis $H_0$ is that the observations follow the model. We reject $H_0$ if we apply the White test under the asymptotic distribution of the test statistic ($p$ value $\approx 0.2\%$). However, when using the bootstrapped distribution, we do not reject the null hypothesis, which is consistent with the result of the Kolmogorov-Smirnov (KS) test on the empirical copula process.

## 8.4 Choice of $X_{M-\Delta}$

We apply a conditional sampling of $X_{M-\Delta}$ knowing $\epsilon_{sim}$. The resulting levels are quite coherent with those observed in the dataset (Fig. 25) and simulations provide as expected fill gaps where data observations are not available.
<figure><img src="image_28.png" /><figcaption>Figure 25: Relation between $\ell(\epsilon_M)$ and $\ell(\tilde{X}_{M-\Delta})$ for $S$ (blue dots: data points, orange triangles: coordinates obtained in the simulations)</figcaption></figure>

# 9 Complements on the consistency of simulations with the observations

## 9.1 Correlation of extremes

The extremogram is defined for extreme time series with

$$ \pi(s', s) = \lim_{q \to 1} \mathbb{P}(\tilde{X}_M^{s'} > u_q(s') | \tilde{X}_M^s > u_q(s)) \quad (12) $$

where $u_q(s) = F_s^{-1}(q)$ is the quantile function of order $q$ of $\tilde{X}_M^s$. We use $q = 0.9$ to compare our simulations with the extreme time series. This tool is based on the hypothesis of asymptotic dependence between the values of the same tidal cycle. We can see for example that the first and the last value of a time series are asymptotically dependent (Fig. 26) as $\chi$ does not converge to 0 when $u$ goes to 1.
## 9.2 Return level estimation

<figure><img src="image_29.png" /><figcaption>Figure 26: Measures of asymptotic dependence between $\tilde{X}_M^1$ and $\tilde{X}_M^{37}$ for the surge with confidence band (grey lines). The blue lines represent the theoretical bounds of the correlation coefficient.</figcaption></figure>

## 9.2 Return level estimation

To obtain the levels obtained in the simulations, we use the formula

$$ \mathbb{P}(A | B) = \mathbb{P}(A \cap C | B) + \mathbb{P}(A \cap \bar{C} | B) \quad (13) $$

with $A = (\tilde{X}_M^t \le x)$, $B = (\tilde{X}_M^t > u^t)$ and $C = (\ell(\mathcal{T}(\varepsilon_M)) > u_\ell)$. The first member on the right-hand side is directly estimated in the simulated time series whereas we use observations to estimate the second member. Then, we use the model quantiles as $x$ levels. The probability $p = \mathbb{P}(A | B)$ is associated to a return period. The link between the quantile order $p$ and the return period $P$ is given by the formula

$$ P = \frac{1}{\text{npy}(1 - p)} \quad (14) $$

where npy is the average number of extreme observations per year, chosen at 7 events per year. Thus, looking at Fig. 15 (left panel), if we observe one day a surge of 0.6m at tidal peak ($t = 19$ for every cycle), we will wait in average 20 years before seeing an event reaching this extreme level.

## 9.3 Extreme aspects

We can analyze the extreme values obtained for other values of $t$ in the simulations and in the data (Fig. 27). The levels obtained are consistent with the theoretical values.
9.4 Distribution of $\ell(\cdot)$ for simulated and observed extreme time series

<figure><img src="image_30.png" /><figcaption>Figure 27: Return level for simulations (orange triangles) and observations (blue dots), for two different values of $t$. Left panel: $t = 25$ (one hour after the tidal peak), right panel: $t = 31$ (two hours after the peak). The dotted lines represent 95% confidence bands.</figcaption></figure>

9.4 Distribution of $\ell(\cdot)$ for simulated and observed extreme time series

We describe here what we obtain for the $\ell(\cdot)$ function in the simulations and in the extreme time series (Fig. 28). The shape of the two distributions are quite similar but their extreme values are quite different. While the maximum value in the data is slightly above 0.6m, the maximum value obtained in the simulation is larger than 0.8m.

<figure><img src="image_31.png" /><figcaption>Figure 28: Distribution of $\ell(S)$ (orange: simulated time series, blue: recorded extreme time series)</figcaption></figure>

10 Choice of $\tilde{X}_{M-\Delta}$

Here, we present the effect of using an unconditional sampling of $\tilde{X}_{M-\Delta}$.
## 10.1 Comparison of percentile levels

### 10.1 Comparison of percentile levels

Firstly, we analyze the quantiles obtained for each value of $t$ in the simulated time series and in the extreme observations (Fig. 29). We see that for several orders the quantile obtained in the simulations do not belong to the confidence interval of the extreme observations.

<figure><img src="image_32.png" /><figcaption>Figure 29: Percentiles obtained in the data and in the simulations without the conditional sampling (blue lines: data, orange dotted lines: simulations)</figcaption></figure>

### 10.2 Comparison of coordinates in a PCA basis

We use the PCA decomposition of the angle $\tilde{X}/\ell(\tilde{X})$ of the extreme observations with the unconditional sampling of $\tilde{X}_{M-\Delta}$. Figure 30 shows that the quantiles of simulated coordinates do not behave like the observed coordinates as empirical quantiles are deviating from the bisector. As a consequence, the observations' and simulations' coordinates do not follow the same law for the first dimension with a $p$ value of 0.01% for the KS test.
10.3 Comparison of the distribution upper tails

<figure><img src="image_33.png" /><figcaption>Figure 30: Q-Q plot for $C_{sim,1}$ (orange triangles) on data scale (blue line) with the unconditional sampling of $\tilde{X}_{M-\Delta}$</figcaption></figure>

10.3 Comparison of the distribution upper tails

We compare the behaviour of extreme values for the simulated and recorded extreme time series when unconditional sampling is used. We see for instance that the $\ell$-extremogram estimated in the simulations is further from the data estimation with conditional sampling than what we obtain with conditional sampling (Fig. 31).

Consequently, the selection method does have an impact on the result.

<figure><img src="image_34.png" /><figcaption>Figure 31: Extremogram for S for simulations (orange) and observations (blue) with confidence bands in dotted lines. Values obtained with the unconditional sampling</figcaption></figure>

10.4 Comparison of predictions with classification algorithms

We focus on the performances of the classifiers (Table 3). The models are able to identify
the simulated time series as 50% does not belong to the confidence interval.
<table><caption>Table 3: Confidence intervals at order 90% for the global accuracy</caption><thead><tr><th>Input</th><th>SVM</th><th>GLM</th><th>RForest</th></tr></thead><tbody><tr><td>Hyper-params</td><td>kernel : radial</td><td>link : logit</td><td>Nb-trees : 500</td></tr><tr><td>X</td><td>53-65</td><td>49-62</td><td>51-65</td></tr><tr><td>1(X)</td><td>51-62</td><td>50-62</td><td>46-56</td></tr><tr><td>A(X)</td><td>54-65</td><td>47-58</td><td>62-72</td></tr></tbody></table>

References

Hirotugu Akaike. A new look at the statistical model identification. *IEEE transactions on automatic control*, 19(6):716-723, 1974. doi: 10.1109/TAC.1974.1100705.

August A Balkema and Laurens De Haan. Residual life time at great age. *The annals of probability*, 2(5):792-804, 1974. doi: 10.1214/aop/1176996548.

Eric Chaumillon, Xavier Bertin, André B Fortunato, Marco Bajo, Jean-Luc Schneider, Laurent Dezileau, John Patrick Walsh, Agnès Michelot, Etienne Chauveau, Axel Créach, et al. Storm-induced marine flooding: Lessons from a multidisciplinary approach. *Earth-Science Reviews*, 165:151-184, 2017. doi: 10.1016/j.earscirev.2016.12.005.

Stephan Clémençon, Nathan Huet, and Anne Sabourin. Regular variation in hilbert spaces and principal component analysis for functional extremes. *Stochastic Processes and their Applications*, 174:104375, 2024. doi: 10.1016/j.spa.2024.104375.

Stuart Coles, Janet Heffernan, and Jonathan Tawn. Dependence measures for extreme value analyses. *Extremes*, 2:339-365, 1999. doi: 10.1023/A:1009963131610.

Stuart Coles, Joanna Bawa, Lesley Trenner, and Pat Dorazio. *An introduction to statistical modeling of extreme values*, volume 208. Springer, 2001. doi: 10.1007/978-1-4471-3675-0.

Claudia Czado and Thomas Nagler. Vine copula based modeling. *Annual Review of Statistics and Its Application*, 9(1):453-477, 2022. doi: 10.1146/annurev-statistics-040220-101153.

Richard J Dawson, Mark E Dickson, Robert J Nicholls, Jim W Hall, Mike JA Walkden, Peter K Stansby, Mustafa Mokrech, Julie Richards, Jianguo Zhou, Jessica Milligan, et al. Integrated analysis of risks of coastal flooding and cliff erosion under scenarios of long term change. *Climatic Change*, 95(1):249-288, 2009. doi: 10.1007/s10584-008-9532-8.

Raphaël de Fondeville and Anthony C Davison. High-dimensional peaks-over-threshold inference. *Biometrika*, 105(3):575-592, 2018. doi: 10.1093/biomet/asy026.

Raphaël de Fondeville and Anthony C Davison. Functional peaks-over-threshold analysis. *Journal of the Royal Statistical Society Series B: Statistical Methodology*, 84(4):1392-1422, 2022. doi: 10.1111/rssb.12498.

Laurens De Haan and Sidney Resnick. On asymptotic normality of the hill estimator. *Stochastic Models*, 14(4):849-866, 1998. doi: 10.1080/15326349808807504.
REFERENCES

Lorenzo Dell'Oro and Carlo Gaetan. Flexible space-time models for extreme data. *Spatial Statistics*, 68:100916, 2025. doi: 10.1016/j.spasta.2025.100916.

Clément Dombry and Mathieu Ribatet. Functional regular variations, pareto processes and peaks over threshold. *Statistics and its Interface*, 8(1):9–17, 2015. doi: 10.4310/SII.2015. v8.n1.a2.

Ana Ferreira and Laurens de Haan. The generalized pareto process; with a view towards application and simulation. *Bernoulli*, pages 1717–1737, 2014. doi: 10.3150/13-BEJ538.

Alan Genz and Frank Bretz. *Computation of Multivariate Normal and t Probabilities*. Lecture Notes in Statistics. Springer-Verlag, Heidelberg, 2009. doi: 10.1007/978-3-642-01689-9.

Eric Gilleland and Richard W. Katz. extRemes 2.0: An extreme value analysis package in R. *Journal of Statistical Software*, 72(8):1–39, 2016. doi: 10.18637/jss.v072.i08.

B. Gouldby, F.J. Méndez, Y. Guanche, A. Rueda, and R. Mínguez. A methodology for deriving extreme nearshore sea conditions for structural design and flood risk analysis. *Coastal Engineering*, 88:15–26, 2014. ISSN 0378-3839. doi: 10.1016/j.coastaleng.2014.01. 012.

Trevor Hastie, Robert Tibshirani, and Jerome Friedman. The elements of statistical learning, 2009.

Henrik Hult and Filip Lindskog. Extremal behavior of regularly varying stochastic processes. *Stochastic Processes and their Applications*, 115(2):249–274, 2005. doi: 10.3150/ 19-BEJ1113.

Déborah Idier, Jérémy Rohmer, Rodrigo Pedreros, Sylvestre Le Roy, Jérome Lambert, Jessie Louisor, Goneri Le Cozannet, and Erwan Cornec. Coastal flood: a composite method for past events characterisation providing insights in past, present and future hazards—joining historical, statistical and modelling approaches. *Natural Hazards*, 101(2):465–501, 2020. doi: 10.1007/s11069-020-03882-4.

Mihyun Kim and Piotr Kokoszka. Extremal correlation coefficient for functional data. *arXiv preprint arXiv:2405.17318*, 2024. doi: 10.48550/arXiv.2405.17318.

Piotr Kokoszka, Stilian Stoev, and Qian Xiong. Principal components analysis of regularly varying functions. *Bernoulli*, 25(4B):3864–3882, 2019. doi: 10.3150/19-BEJ1113.

Sébastien Lê, Julie Josse, and François Husson. FactoMineR: A package for multivariate analysis. *Journal of Statistical Software*, 25(1):1–18, 2008. doi: 10.18637/jss.v025.i01.

Andy Liaw and Matthew Wiener. Classification and regression by randomforest. *R News*, 2 (3):18–22, 2002. doi: 10.32614/cran/package.randomforest.

David Lopez-Paz and Maxime Oquab. Revisiting classifier two-sample tests. *arXiv preprint arXiv:1610.06545*, 2016. doi: 10.48550/arXiv.1610.06545.
Ed Mackay, Guillaume de Hautescorcque, Erik Vanem, and Philip Jonathan. The effect of serial correlation in environmental conditions on estimates of extreme events. *Ocean Engineering*, 242:110092, 2021. doi: 10.1016/j.oceaneng.2021.110092. URL https://www.sciencedirect.com/science/article/pii/S0029801821014189.

David Meyer, Evgenia Dimitriadou, Kurt Hornik, Andreas Weingessel, and Friedrich Leisch. *e1071: Misc Functions of the Department of Statistics, Probability Theory Group (Formerly: E1071)*, TU Wien, 2023.

Thomas Nagler, Ulf Schepsmeier, Jakob Stoeber, Eike Christian Brechmann, Benedikt Graeler, and Tobias Erhardt. *VineCopula: Statistical Inference of Vine Copulas*, 2023.

Thomas Opitz, Denis Allard, and Gregoire Mariethoz. Semi-parametric resampling with extremes. *Spatial Statistics*, 42:100445, 2021. doi: 10.1016/j.spasta.2020.100445.

R Core Team. *R: A Language and Environment for Statistical Computing*. R Foundation for Statistical Computing, Vienna, Austria, 2023. URL https://www.R-project.org/.

Lena Reimann, Athanasios T Vafeidis, and Lars E Honsel. Population development as a driver of coastal risk: Current trends and future pathways. *Cambridge Prisms: Coastal Futures*, 1:e14, 2023. doi: 10.1017/cft.2023.3.

Mathieu Ribatet and Christophe Dutang. *POT: Generalized Pareto Distribution and Peaks Over Threshold*, 2022.

Jeremy Rohmer, Deborah Idier, Remi Thieblemont, Goneri Le Cozannet, and François Bacchoc. Partitioning the contributions of dependent offshore forcing conditions in the probabilistic assessment of future coastal flooding. *Natural Hazards and Earth System Sciences*, 22(10):3167-3182, 2022. doi: 10.5194/nhess-22-3167-2022.

Ludger Rüschendorf. On the distributional transform, sklar's theorem, and the empirical copula process. *Journal of statistical planning and inference*, 139(11):3921-3927, 2009.

Kosuke Sando, Ryota Wada, Jérémy Rohmer, and Philip Jonathan. Multivariate spatial and spatio-temporal models for extreme tropical cyclone seas. *Ocean Engineering*, 309:118365, 2024. doi: 10.1016/j.oceaneng.2024.118365.

Ulf Schepsmeier. A goodness-of-fit test for regular vine copula models. *Econometric Reviews*, 38(1):25-46, 2019. doi: 10.1080/07474938.2016.1222231.

Sebastián Solari, Marta Egüen, María José Polo, and Miguel A Losada. Peaks over threshold (pot): A methodology for automatic threshold estimation using goodness of fit p-value. *Water Resources Research*, 53(4):2833-2849, 2017. doi: 10.1002/2016WR019426.

Stan Tendijck, Philip Jonathan, David Randell, and Jonathan Tawn. Temporal evolution of the extreme excursions of multivariate k-th order markov processes with application to oceanographic data. *Environmetrics*, 35(3):e2834, 2024. doi: 10.1002/env.2834.
Michael Thomas. Statistical analysis of extreme values: with applications to insurance, finance, hydrology, and other fields. Birkhäuser Verlag, 2001. doi: 10.1007/978-3-7643-7399-3_6.

Alexandra Toimil, Iñigo J Losada, Pedro Díaz-Simal, Cristina Izaguirre, and Paula Camus. Multi-sectoral, high-resolution assessment of climate change consequences of coastal flooding. *Climatic Change*, 145(3):431-444, 2017. doi: 10.1007/s10584-017-2104-z.

Halbert White. Maximum likelihood estimation of misspecified models. *Econometrica: Journal of the econometric society*, pages 1-25, 1982. doi: 10.2307/1912526.