# Multivariate spatial and spatio-temporal models for extreme tropical cyclone seas

Kosuke Sando, Ryota Wada, Jérémy Rohmer, Philip Jonathan

► To cite this version:

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

HAL Id: hal-04668283

https://brgm.hal.science/hal-04668283v1

Submitted on 6 Aug 2024

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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Research paper

# Multivariate spatial and spatio-temporal models for extreme tropical cyclone seas

Kosuke Sando<sup>a</sup>, Ryota Wada<sup>a,*</sup>, Jérémy Rohmer<sup>b</sup>, Philip Jonathan<sup>c,d</sup>

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

## ARTICLE INFO

Keywords:

Multivariate extremes

Time-series

Tropical cyclone

Spatial extremes

## ABSTRACT

Estimates of extreme environments and responses of offshore structures for tropical cyclone conditions are typically made using time-series of ocean environmental data, hence helping to ensure safe structural design. However, estimates are often subject to large uncertainties because of the short length of available time-series. We propose a methodology to characterise extreme multivariate time-series for tropical cyclones, by extending the STM-E spatial extreme value model of Wada et al. (2018) to incorporate (a) storm peaks of multiple metocean variables, using the conditional extremes model of Heffernan and Tawn (2004) (leading to MSTM-E methodology), and additionally (b) time-series evolution around the storm peak, using a history-matching approach (leading to MSTM-TE). We use both MSTM-E and MSTM-TE to estimate the return values of multivariate extremes from synthetic cyclone data for a spatial neighbourhood of locations offshore Guadeloupe (in the Lesser Antilles). The comparison of storm peak analysis using MSTM-E against single location conditional model shows the benefit of MSTM-E in reducing return value variance without sacrificing bias, in both marginal and joint extremes. Moreover, characteristics of multivariate time-series realisations generated under fitted MSTM-TE models (with 200 years of data) are shown to be in good agreement with those of the original time-series data used to fit the model (with 1000 years of data).

## 1. Introduction

A quantitative description of extreme environments and corresponding structural responses is essential for design and reliability analysis; for example, we might adopt the *N*-year structural response as ultimate limit state for design, typically with *N* = 10<sup>m</sup>, *m* = 2, 3, 4. In ocean basins exposed to tropical cyclones (TCs), extreme responses are produced by one or more of severe winds and waves. To establish a good statistical model for extreme conditions, a reasonable amount of representative data are needed. Since extreme tail behaviour of winds and waves under cyclones is typically not strongly informed by the behaviour of typical non-extreme winds and waves, only data for extreme conditions are useful to characterise joint extremes (Wada et al., 2016). Moreover, the frequency of TC events is low compared to that of extra-tropical storms. In addition, the spatial extent of TCs is smaller than that of extra-tropical storms (Jonathan and Ewans, 2007). It is apparent therefore that the size of relevant sample data for empirical inference of extreme TC conditions is often limited, making precise estimation of quantities such as return values challenging. For

ocean basins with offshore oil and gas production, hindcast simulation data of the ocean environment are sometimes available. The GOMOS (Gulf of Mexico Oceanographic Study) hindcast for the Gulf of Mexico is an exceptionally long example, covering the period of 1900–2008, including 379 hurricane events (Cardone and Cox, 2013). However, most ocean basins lack such long-duration hindcasts. Thus, design for TC conditions typically requires characterisation of extreme events with return periods longer than the period of available hindcasts or measurements. This further motivates the use of extreme value models for the tails of distributions. Estimation of univariate metocean extremes using extreme value methods has been studied widely, for example by Caires and Sterl (2005), Teena et al. (2012) and Wada et al. (2016). The challenges of modelling TCs, and in particular the parameter uncertainties from inference using small samples is noted by Jonathan and Ewans (2007) and Wada and Waseda (2020). There are several approaches to tackle the extrapolation problem for TCs. In the current work, we simulate under an extreme value model for the tail of the joint distribution of hindcast data. An alternative approach is so-called

<footnote><sup>a</sup> Graduate School of Frontier Sciences, The University of Tokyo, Tokyo, Japan</footnote>

<footnote><sup>b</sup> BRGM, Orleans, France</footnote>

<footnote><sup>c</sup> Shell Research Limited, London, United Kingdom</footnote>

<footnote><sup>d</sup> Department of Mathematics and Statistics, Lancaster University, Lancaster, United Kingdom</footnote>

<footnote>* Corresponding author.
E-mail address: r_wada@k.u-tokyo.ac.jp (R. Wada).</footnote>
**stochastic simulation.** For example, Bloemendaal et al. (2020) generate 10,000 years of synthetic TCs from estimates of probability densities for key parameters of TC evolution, such as direction, translation speed, and wind speed, derived empirically from historical storm track data. The accuracy of the stochastic simulation depends on the adequacy of the description of cyclone parameters, and whether fitted densities are informative for extreme cyclones.

Offshore structures in general are exposed to a range of environmental phenomena, including waves, winds, currents, sea ice, etc. Structural response is often driven by multiple environmental variables. For example, the pitch of a vessel is known to depend on significant wave height $H_S$ and period $T$. Therefore, reliable quantification of extreme structural responses requires a statistical model for the joint variation of metocean drivers. Multivariate extreme value theory is concerned with the characterisation, estimation and extrapolation of tails of multidimensional distributions. Estimation of joint extremes is more challenging than that of univariate extremes, requiring estimation of all marginal extreme value distributions as well as the dependence between variates (Jonathan and Ewans, 2013). We note the large statistical literature on multivariate extremes (reviewed by Dutfoy, 2021; Walshaw, 2023) and spatial extremes (reviewed by Davison et al., 2012; Tawn et al., 2018 and Huser and Wadsworth, 2022). Statistical modelling of multivariate and spatial extremes is complicated by the fact that, in contrast to the univariate situation, no general limiting distribution for extremes exists (e.g. Jonathan and Ewans, 2013). The key assumption of max-stability underpinning univariate extreme value theory implies that the general form of the limiting distribution for extreme values does not depend on the underlying distribution from which the sample is drawn. A less general form of max-stability motivates much of the existing methodology on multivariate and spatial extremes, but the models generally require large samples for fitting. Furthermore, it is not generally clear whether the restrictive assumptions of the methodology are appropriate for the sample (Huser et al., 2024), which potentially can lead to large bias in the estimation of return values and related quantities. There are specific examples of bivariate model fits to oceanographic data used for design using $H_S$ and peak period $T_p$, including Haver (1987), Ferreira and Soares (2002) and Dong et al. (2013). These studies are typically conducted in non-cyclone dominated ocean basins; application to cyclone-dominated regions would be more challenging. Tendijck et al. (2023b) attempts to characterise such models asymptotically, and provide some justification for their use. In this work we adopt the conditional extremes model of Heffernan and Tawn (2004) to characterise the conditional distribution of extremes of a multivariate quantity given that at least one of its components takes an extreme value. The conditional extremes model provides a pragmatic regression-type methodology for applied multivariate extreme value analysis, and admits both asymptotic dependence and asymptotic independence between variables. In particular, its modelling assumptions are relatively easily confirmed using model fit diagnostics.

Various approaches to utilise data from different spatial locations to increase the effective sample size for statistical analysis have been discussed in the ocean engineering literature. A straightforward approach is pooling of data from a neighbourhood of locations (Heideman and Mitchell, 2009), assuming they all correspond to observations from the location of interest (over an extended period of time), whilst ignoring the spatial dependence present. Estimates of extreme quantiles and their uncertainties from an extreme value model fit need to be carefully interpreted, since the model wrongly assumes independent observations. There is a large statistical literature on spatial extremes, which can as noted above be thought of as a particular extension of multivariate extremes to accommodate spatially-dependent observations. These models tend to be unwieldy: problematic to estimate reliably and hence of limited practical use; Ross et al. (2017) presents one of the few applications of spatial extremes in ocean engineering. Another approach to exploit spatial data, particularly for TCs, is track-shifting (Vickery et al.,

2000) or stochastic generation models (Bloemendaal et al., 2020). Here, a large number of sample data sets are synthesised with shifted cyclone tracks or random selections of storm features from past hurricanes. Of course, the validity of the track-shifted or stochastic cyclone data depends on the quality of the track modelling algorithms or stochastic process employed. In this work, we extend the simple spatial model (STM-E) for extreme waves in TC seas introduced by Wada et al. (2018) to the multivariate context. The underpinning STM-E model assumes that the spatial distribution of significant wave height data for a given cyclone can be described in terms of two independent quantities: (a) the peak severity or Spatio-Temporal Maximum (STM) of the cyclone, and (b) the spatial Exposure (E) of the cyclone, which quantifies the maximum influence of the cyclone at a given location, expressed as the ratio of maximum significant wave height at that location (over the period of the cyclone) to the STM. The assumption of independence implies that the characteristics of Exposure do not change with STM. If shown to be a reasonable assumption for a given sample of data, we can then exchange STM and Exposure between different storms, leading to a mechanism to simulate extreme storms over extended periods. Using the conditional extreme model, we extend STM-E to characterise extremes of multiple variables such as significant wave height and wind speed. The resulting methodology is referred to as MSTM-E (Multivariate Spatio-Temporal Maximum and Exposure).

Typically, extreme value models are fitted to observations assumed conditionally independent in time. Therefore, it is usual to estimate extreme value models for observations of maxima or peaks of storm events, which can be considered approximately independent. Simulations of time-series under such a model can then be achieved by first simulating storm peak events, then coupling those events to time-series (possibly adjusted, see Hansen et al., 2020) of historical storm evolution for storms with similar peak characteristics to the simulated peak events. In this way, we achieve simulations of extreme multivariate time-series without explicitly adopting a model for time-series evolution. We again note the statistical literature of extremes of time-series (reviewed by Chavez-Demoulin and Davison, 2012) and specific models for multivariate time-series of serially-correlated extreme values, such as those discussed by Tendijck et al. (2023a), but these are research tools which have not yet been stress-tested for practical application. In this work, we further extend the MSTM-E methodology outlined above to the time-series context incorporating a history-matching approach similar to that of Hansen et al., 2020, providing a method for simulating coupled time-series of multiple oceanographic variables corresponding to severe TC conditions. The extended approach is referred to a MSTM-TE (Multivariate Spatio-Temporal Maximum and Temporal Exposure).

The importance of capturing temporal evolution of storms for reliable design has been noted by Brown et al. (1984) and Anderson et al. (2019) in the context of wind turbines and coastal defences respectively. MSTM-TE thus extends MSTM-E to incorporate the temporal evolution of the multivariate environment throughout the extreme event. Whereas MSTM-E considers the maximum value for each metocean variable during the extreme event only, MSTM-TE considers the full multivariate time-series for the extreme TC event. As a result, MSTM-TE offers a number of advantages in estimating the extreme response. Specifically, (a) the component maxima for an MSTM event are typically not contemporaneous, whereas an MSTM-E analysis effectively assumes that they are. Hence, return values and similar estimates from an MSTM-E model may tend to be conservative relative to those from MSTM-TE; and (b) TC duration (captured by MSTM-TE but not MSTM-E) influences the maximum response observed (Mackay et al., 2021).

## Aims and outline

This paper seeks to improve the estimation of multivariate extreme environmental conditions for TC conditions. By exploiting spatial information wisely, it has previously been demonstrated (Wada et al.,
2018, 2022) that STM-E improves the precision of univariate return value estimation, relative to site-specific estimation, without introducing bias. We anticipate that, by exploiting spatial information similarly, the new multivariate models (MSTM-E, MSTM-TE, referred to together as MSTM-(T)E) will provide better understanding of extremal dependence between metocean variables, and hence more precise estimates of return values, conditional return values and environmental design contours, relative to site-specific estimation. When the design of a marine structure requires knowledge of the multivariate temporal structure of TCs, it is reasonable to expect that MSTM-TE will provide a better basis for design, given that it can be well estimated from data. We will quantify the relative performance of different models by comparing estimates for extreme quantiles with known estimates from long-term TC simulations over the Caribbean Sea.

The article is structured as follows. Section 2 introduces the motivating application in the region of Guadeloupe in the Caribbean Sea. Section 3 then outlines the multivariate spatial MSTM-E and MSTM-TE methodology, including recommendations on diagnostics to support good model fit. Section 4 then applies MSTM-E to the Caribbean data together with the illustration of MSTM-TE application. Section 5 draws conclusions and considers some of the wider implications of the work. On-line Supplementary Material (referenced below with prefix “S”) provides supporting arguments and evidence; see the Acknowledgements section for the appropriate hyperlink.

## 2. Motivating application

The study area (60.5–62.5W°, 15.0N–17.0N°) for the current work is the French overseas region of Guadeloupe, an archipelago located in the southern part of the Leeward Islands (see Fig. 1). This region has been impacted by many devastating TCs in the past; Hugo in 1989 and Maria in 2017 being the more damaging recent ones. Analysis of the HURDAT database (Landsea and Franklin, 2013) reveals that approximately 0.6 cyclones per year passed within 300 km of the study area on average for the period 1970–2019.

For clarity of discussion in this work, we refer to four spatial domains. The first is the *Leeward region* corresponding to the full domain on which TCs are simulated, shown in the left hand panel of Fig. 1. The second is the *Guadeloupe region* corresponding to the subdomain of the Leeward region on which MSTM-(T)E models are developed, shown in the right hand panel of Fig. 1. The third and fourth regions are the *west subregion* and *east subregion*, lying to the west and east of longitude 61.5W respectively, corresponding to a binary partition of the Guadeloupe region. As will be explained later in Section 4.1, it will also be convenient to identify cyclones whose maximum $H_s$ occurs at a location in the east (west) subregion. These will be referred to as *east cluster* (*west cluster*) events.

The data used in this study correspond to the output of a TC simulator for the Guadeloupe region, available at https://carib-coast.brgm.fr/, provided by the French Geological Survey (BRGM). The data provide time-series of significant wave height and wind speed in particular on a mesh of locations around Guadeloupe, corresponding to 685 TCs passing within 300 km of the islands over a period of 1000 years.

The wave simulator is based on data for key TC characteristics (track positions over time, and time evolution of radius, maximum wind speed and central atmospheric pressure) from stochastically generated TCs (Bloemendaal et al., 2020). TC characteristics are translated into two-dimensional surface wind and pressure fields over the entire cyclone course (Holland, 1980). These fields provide input to a Wavewatch 3 model (Tolman, 2014, henceforth WW3) generating offshore waves propagating to the Guadeloupe coastline. WW3 version 4.18 with source term package described by Ardhuin et al. (2010) is adopted, with discretisation into 32 wave frequencies and 36 directions. The computational grid covers a large part of the Antilles (50.0–84.7W°, 8.4–22.1N°) with minimum resolution of about 300 m at the coast and a maximum resolution of about 20 km offshore.

Bathymetry is computed from GEBCO (https://www.gebco.net), and from SHOM (https://data.shom.fr/donnees) for the French West Indies.

A brief exploratory analysis of the data, including plots of locations of STM events, histograms of STM magnitudes, and TC tracks is provided in section S.2 of the Supplementary Material.

## 3. Methodology

This section describes the methodology underpinning MSTM-E and MSTM-TE, both which can be thought of as extensions of the STM-E model of Wada et al. (2018), achieved by incorporation of the conditional extremes model of Heffernan and Tawn (2004) (for multivariate extremes in MSTM-E) and also a time-series history matching procedure (for multivariate time-series extremes in MSTM-TE). For completeness, an outline of the original STM-E model is provided in the Supplementary Material.

### 3.1. Overview

The modelling procedure is described sequentially in the subsections below. The first step is the isolation of so-called multivariate spatio-temporal maxima (MSTM) and exposure (denoted E; TE for temporal exposure) sets for tropical cyclones; these entities are defined in Section 3.2. In essence, the MSTM is a vector of maximum values of metocean variables associated with a TC. The associated exposure sets describe the spatial characteristics of the TC (e.g. in time and space) relative to its MSTM. The key assumption for both MSTM-E and MSTM-TE models, which must be carefully justified in application, is that it is reasonable to characterise the distribution of (a) extreme MSTM and (b) exposure sets independently of each other. Estimating the distribution of extreme MSTM requires a multivariate extreme value model, since extremes of the metocean variables of which it is comprised will be correlated in general. In the current work, we adopt a combination of marginal extreme value models and the conditional extremes model of Heffernan and Tawn (2004) to achieve this. These are outlined in Section 3.3. The full modelling procedure is summarised diagrammatically in Fig. 2. The modelling procedure is empirical in nature. It is therefore critical to assess, for the application at hand, whether modelling assumptions appear to be valid. Specifically, reasonable inference from extreme value models requires careful specification of threshold values; we describe sensitivity analyses of inferences to threshold choice in Section 3.5. Another critical model assumption is the independence of (M)STM and (T)E; in Section 3.6, we discuss the methods used to assess this. Once estimated, we use the fitted MSTM model to generate new realisations of MSTM vectors, which are then combined with historical occurrences of exposure sets to provide spatial and spatio-temporal realisations of TCs. Using these, we can estimate the distributional characteristics of all metocean variables of interest, including marginal return values for individual metocean variables, joint and conditional return values for multiple metocean variables etc., whilst capturing the spatial and spatio-temporal dependence of TCs in a pragmatic but principled manner. The procedure for merging realisations of MSTM with historical exposure sets is described in Section 3.7.

The approach to assessing the performance of MSTM-E and MSTM-TE models against competitor models and long TC simulations, by comparing estimates of return values and environmental contours, is outlined in Section 3.8.

### 3.2. Isolating MSTM and exposure sets

In general, consider a data sample for *D* metocean variables describing the evolution of TCs in space and time. We refer to the observation of variable *d* (*d* = 1, 2, ..., *D*) for cyclone *n* (*n* = 1, 2, ..., *N*) at location *m* (*m* = 1, 2, ..., *M*) and relative time point *t* (*t* = 1, 2, ..., *L*<sub>n</sub>) as $\tilde{x}_{n,m}^d(t)$ (with “tilde notation” always reserved for functions of time). Note that
<figure><img src="image_4.png" /><figcaption>Fig. 1. Leeward Islands and Guadeloupe. Left: the full spatial domain of the Leeward Islands (referred to for clarity as the Leeward region), with red box highlighting the Guadeloupe region of interest. Two tracks from cyclone events used for Fig. 3 are also plotted; note that the Leeward region on which cyclones are simulated, is larger than the Guadeloupe region of interest in this work. Right: magnified view of the Guadeloupe region including the islands of Guadeloupe, with simulation grid locations (black dots) and five locations (red crosses) to be considered in detail. In order, the latitude–longitude coordinates of Locations 1 to 5 are (61.87W°,15.96N°), (61.07W°,16.51N°), (61.06W°,15.98N°), (61.86W°,16.51N°), (61.47W°,16.07N°). The vertical blue line at longitude 61.5W is used to partition the set of cyclones into west and east clusters.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 2. A diagrammatic overview of MSTM-E and MSTM-TE estimation procedures. Red boxes are actions, and blue boxes are data. Blue text denotes actions or data for the MSTM-TE procedure.</figcaption></figure>

a relative time indicator *t* is adopted such that the first observation of each TC occurs at *t* = 1, and *L*<sub>n</sub> differs for each event. Further, applications in Section 4 below will consider the case *D* = 2, for variables *H*<sub>S</sub> and *U*<sub>10</sub>. Temporal maxima *x*<sub>n,m</sub><sup>d</sup>, and spatio-temporal maxima (STM) *s*<sub>n</sub><sup>d</sup> for each dimension *d*, are defined by

$$x_{n,m}^{d} = \max_{t=1,...,L_n} \tilde{x}_{n,m}^{d}(t) \text{ for } n=1,2,...,N, m=1,2,...,M \text{ and } d=1,...,D \quad (1)$$

and

$$s_n^d = \max_{m=1,...,M} x_{n,m}^d \text{ for } n=1,2,...,N \text{ and } d=1,...,D. \quad (2)$$

That is, a temporal maximum is the maximum observed value of a variable at a location for a given TC, and the STM is the maximum observed value of a variable over all locations for a given TC. For the TC indexed by *n*, the set of *D* STMs {*s*<sub>n</sub><sup>d</sup>}<sub>d=1</sub><sup>D</sup> is referred to collectively as a multivariate spatio-temporal maximum (MSTM) vector. Furthermore, we define the exposure (E) at location *m* relative to STM *s*<sub>n</sub><sup>d</sup> as

$$e_{n,m}^{d} = \max_{t=1,...,L_n} \tilde{e}_{n,m}^{d}(t) \text{ for } n=1,...,N, m=1,2,...,M \text{ and } d=1,...,D \quad (3)$$

where $\tilde{e}_{n,m}^d(t)$ is the temporal exposure (TE) defined by

$$\tilde{e}_{n,m}^{d}(t) = \frac{\tilde{x}_{n,m}^{d}(t)}{s_{n}^{d}} \text{ for } n=1,...,N, m=1,2,...,M, \quad t=1,2,...,L_n \text{ and } d=1,...,D. \quad (4)$$

For TC *n*, the set of Es {*e*<sub>n,m</sub><sup>d</sup>}<sub>d=1</sub><sup>D</sup> or TEs {*e*<sub>n,m</sub><sup>d</sup>(*t*)}<sub>d=1,t=1</sub><sup>D,L<sub>n</sub></sup> over the *D* variables are referred to as exposure sets. Exposure sets will also show dependence between constituent variables in general.

Using these definitions, we isolate MSTM and exposure sets for each of *N* TCs in the sample. Illustrations of MSTM and temporal exposure sets for the two TCs illustrated in Fig. 1 are shown and described in Fig. 3. STMs of the *D* different variables do not necessarily occur contemporaneously within a TC. Further, since our interest lies in characterising the behaviour of extreme TCs, we include only cyclones with STM values above some threshold, the choice of which is discussed in Section 3.6.

### 3.3. A model for MSTM

We require a model for the MSTM vector *S* with components *S*<sub>*d*</sub> (*d* = 1, 2, ..., *D*), when at least one of its components is large. We use the sample of *N* MSTM vectors {*s*<sub>*n*</sub><sup>*d*</sup>}<sub>*d*=1</sub><sup>*D*</sup> (*n* = 1, 2, ..., *N*) to estimate
<figure><img src="image_6.png" /><figcaption>Fig. 3. Examples of TC events to aid in understanding MSTM-E and MSTM-TE. Each row of panels correspond to a cyclone event, the track of which is given in Fig. 1, with the upper (lower) panel corresponding to the cyclone track in red (orange) in Fig. 1. Left to right: the corresponding MSTM (panel 1), TE at location 2 (see Fig. 1, panel 2), and the corresponding observed timeseries of ($H_s, U_{10}$) in time (panel 3) and jointly (panel 4). Note that the maximum TE with respect to $H_s$ and $U_{10}$ for both cyclones is strictly <1, indicating that location 2 is not the location of the STM for any combination of metocean variable and cyclone. Note further that TE (panel 2) multiplied by MSTM (panel 1) gives observed time-series (panel 3).</figcaption></figure>

this. The conditional extremes model of Heffernan and Tawn (2004) is an appropriate model to adopt, since it is theoretically well-founded, relatively easy to estimate, and has been used in applications in various fields. Estimation of the conditional extremes model is performed in two stages, involving (a) estimation of marginal models for each $S_d$ independently, followed by transformation to standard marginal Laplace scale, (b) estimation of $D \times (D-1)$ pairwise conditional extremes models for $S_d'|S_d$ ($d, d' = 1, 2, \dots, D, d' \neq d$). These stages are described in Sections 3.3.1 and 3.3.2 respectively below.

### 3.3.1. Modelling margins

The marginal distribution of each of the D components of MSTM is estimated independently. For component $S_d$ ($d = 1, 2, \dots, D$) of MSTM, we assume that sample values {$s_n^d$}$_{n=1}^N$ are independently and identically distributed. Moreover, we assume that the conditional distribution of excesses of $S_d$ above some threshold $u_d$ is the generalised Pareto distribution in Eq. (5), where $u_d$ is the marginal threshold referred to in Fig. 2.

$$ \mathbb{P}(S_d \le x | S_d > u_d) = \begin{cases} 1 - \left(1 + \frac{\xi_d(x-u_d)}{\sigma_d}\right)^{-1/\xi_d} & \text{for } \xi_d \neq 0, \\ 1 - \exp\left(-\frac{x-u_d}{\sigma_d}\right) & \text{for } \xi_d = 0 \end{cases} \quad (5) $$

for shape parameter $\xi_d \in \mathbb{R}$ and scale parameter $\sigma_d > 0$. When $\xi_d < 0$, the upper bound for $S_d$ is restricted to $u_d - \sigma_d/\xi_d$, otherwise $S_d \in \mathbb{R}$. The distribution $F_d^*$ of $S_d$ below threshold $u_d$ is estimated empirically from observed data. The full unconditional distribution $F_d$ of $S_d$ is thus

$$ F_d(x) = \begin{cases} 1 - \{F_d^*(u_d)\} \left\{1 + \frac{\xi_d(x-u_d)}{\sigma_d}\right\}^{-1/\xi_d} & \text{for } x \ge u_d, \\ F_d^*(x) & \text{for } x < u_d. \end{cases} \quad (6) $$

The estimated marginal distribution is then used to transform the sample of values for $S_d$ to the standard Laplace scale variables $S_d^L$ (with superscript L used to notate STM variables on Laplace scale), for each $d = 1, 2, \dots, D$, using

$$ S_d^L = \begin{cases} \log(2F_d(S_d)) & \text{for } F_d(S_d) < \frac{1}{2}, \\ -\log(2(1-F_d(S_d))) & \text{for } F_d(S_d) \ge \frac{1}{2}. \end{cases} \quad (7) $$

### 3.3.2. Modelling dependence

Suppose that $Y = (Y_1, Y_2, \dots, Y_D)$ is a D-dimensional random vector on standard Laplace margins. The conditional extremes model provides a means to describe the joint distribution of $D-1$ of these variables, conditional on the value of the remaining variable being large. For D-vector $v$, writing $v_{-d}$ as the $(D-1)$-vector consisting of all but component $d$, we can describe the conditional extremes model as follows. We write

$$ Z_{-d} := (Y_{-d} - a_{-d}Y_d) / Y_d^{b_{-d}} \quad (8) $$

where component-wise operations are assumed. That is, $Z_{-d}$ is a transformed form of $Y_{-d}$, using normalisation ($D-1$)-vectors $a_{-d} \in [-1, 1]^{(D-1)}$ and $b_{-d} \in (-\infty, 1]^{(D-1)}$. The conditional extremes model assumes that vectors $a_{-d}$ and $b_{-d}$ can be found such that the conditional distribution of $Z_{-d}$ given $Y_d = y$ converges to an unknown non-degenerate ($D-1$)-dimensional limiting distribution $G_{-d}$ for values $y > \psi_d$, as the value $\psi_d$ increases. This motivates the statistical model

$$ (Y_{d'} | (Y_d = y)) = a_{d'|d} y + y^{b_{d'|d}} Z_{d'|d} \quad \text{for } y > \psi_d, \\ \text{and } d, d' = 1, 2, \dots, D, \quad d' \neq d \quad (9) $$

for each pair of variables ($Y_d, Y_{d'}$), with parameters $a_{d'|d} \in [-1, 1]$ and $b_{d'|d} \in (-\infty, 1]$ to be estimated. Since the distribution of $Z_{d'|d}$ is not specified by theory, for model fitting purposes it is assumed Gaussian with mean $\mu_{d'|d} \in \mathbb{R}$ and standard deviation $\kappa_{d'|d} > 0$. Model fitting is thus performed over a partition of a sub-domain of the full D-dimensional Laplace-scale domain into sets $C_d$ ($d = 1, 2, \dots, D$), with set $C_d$ defined as the region of the space in which $Y_d > \psi_d$ and $\max\{Y_{-d} < Y_d; \text{ that is, in set } C_d, \text{ component } Y_d \text{ of Y exceeds } \psi_d, \text{ and } Y_d \text{ is also more extreme than any of the other components. Following model fitting, the set of empirical } (D-1)\text{-dimensional residuals } R_d \text{ from model fits for } Y_{-d} | (Y_d = y) \text{ is adopted as a sample from the joint distribution of } Z_{-d}, \text{ which can be sampled at random with replacement for subsequent simulations under the model.}$

## 3.4. Uncertainty quantification for marginal and conditional extreme value models

A bootstrap resampling scheme is used to quantify the uncertainty of estimated model parameters, for both marginal and conditional extreme value models. In the scheme, for a given model of interest,
$N_B$ different model variants are estimated, each based on a different bootstrap resample of the input data for modelling. Uncertainty in parameter estimation is then represented by the set of bootstrap parameter estimates. The bootstrap procedure is used to aid threshold selections in Section 3.5 in particular. We use $N_B = 1000$ for marginal models, and $N_B = 100$ for conditional extremes models.

In the current work, we choose not to propagate uncertainty through the chain of inferences of which MSTM-(T)E is comprised, but rather to adopt bootstrap mean model parameter estimates at each stage as the best models for subsequent inferences. Overall MSTM-(T)E model performance is quantified by estimating return values for long return periods from random subsamples of data corresponding to shorter periods of observation, as explained in Section 3.8.

### 3.5. Threshold selection

From a practical perspective, a critical aspect of marginal and dependence modelling is the reasonable choice of threshold. Choice of marginal thresholds {$u_d$}$_{d=1}^D$ influences the estimation of each marginal distribution, and hence also the conditional extremes model and the identification of Es and TEs. In this work, as illustrated in Section 4, we examine threshold stability plots for $\xi_d$ (e.g. Coles, 2001) (figs. S.3.1 and S.3.2 in the Supplementary Material) to make informed selections of marginal thresholds, using the bootstrap scheme of Section 3.4 for uncertainty quantification. The values of marginal thresholds are fixed prior to attempting dependence modelling.

The choice of dependency thresholds {$\psi_d$}$_{d=1}^D$ determines which events are considered extreme, affecting the inferred conditional extremes models, and whether exposure appears to be correlated with MSTM. In the notation of Section 3.3.2, for simplicity and in the absence of strong evidence against the reasonableness of the assumption, it was decided to adopt a common value $\psi_d = \psi$ ($d = 1, 2, ..., D$), typically selected to be larger than $\max_{d=1,2,...,D} F_d(u_d)$, so that the conditional extremes threshold can be viewed as more extreme than any marginal threshold. To assess the sensitivity of inference to the choice of {$\psi_d$}$_{d=1}^D$, we evaluate the convex hull of parameter pairs ($a_{d'}^B$, $b_{d'}^B$) under bootstrap resampling (specified with the superscript $^B$) of the Laplace-scale input data for conditional extremes modelling, for difference choices of $d, d'$ ($d, d' = 1, 2, ..., D, d' \neq d$), over a range of dependency thresholds. Illustrations of typical findings are given in Section 4.

### 3.6. Assessing the independence of MSTM and (T)E

The basic STM-E methodology introduced by Wada et al. (2018) assumes that STM and E are approximately independent. In the current context, the corresponding assumption is that occurrences of MSTM and (T)E can be considered approximately independent. That is, an occurrence of MSTM can be validly associated with any occurrence of (T)E. To apply the MSTM-(T)E models reasonably therefore, it is critical to assess whether this assumption is valid for the samples of MSTM and (T)E observed in the current work.

To achieve this, we use Kendall's rank correlation coefficient test, also known as Kendall's tau test. This is a non-parametric statistical hypothesis test to determine whether two variables are statistically dependent. Details of how the test is applied for MSTM-E and MSTM-TE are provided in Sections 3.6.1 and 3.6.2 respectively, and the results can be found in section S.4 of the Supplementary Material.

#### 3.6.1. MSTM-E

In the notation of Section 3.2, Kendall's rank correlation coefficient $\tau_m^d$ for component $d$ ($d = 1, 2, ..., D$) and location $m$ ($m = 1, 2, ..., M$) is given by

$$ \tau_m^d = \frac{1}{N(N-1)} \sum_{n' \neq n} \operatorname{sgn}(s_n^d - s_{n'}^d) \operatorname{sgn}(e_{n,m}^d - e_{n',m}^d) \quad (10) $$

where $N$ is the total number of TC events present, and $\operatorname{sgn}(x) = 1$ for $x > 0$, $= 0$ for $x = 0$ and $= -1$ for $x < 0$. If the value of $\tau_m^d$ is close to zero, then the rank correlation between STM and E is small. In practice, we state that STM and E are independent when the value of $\tau_m^d$ lies within the 2.5% and 97.5% percentiles of the null distribution of the test statistics under independence, for sample size $N$. We can then generate spatial maps for each component $d$ of the environment, indicating locations where the rank correlation of STM and E appears non-negligible. For the Guadeloupe application, results of Kendall's tau tests are summarised in Section 4. In particular, in Section 4.1 we discuss interesting findings regarding non-negligible levels of rank correlation between STM and E, associated with differences between TC characteristics to the west and east of Guadeloupe. Notice also that we only consider the marginal dependence between STM and E per TC component, and do not consider cross-dependence between STM and E for different TC components (although this could also be incorporated).

#### 3.6.2. MSTM-TE

In the notation of Section 3.2, the TE $\bar{\epsilon}_{n,m}^d(t)$ for component $d$ ($d = 1, 2, ..., D$) of TC $n$ ($n = 1, 2, ..., N$) at location $m$ ($m = 1, 2, ..., M$) is a time-series for $t \in [1, 2, ..., L_n]$, the length $L_n$ of which varies with TC. As such, the Kendall's tau test cannot be directly applied to STM and TE without some adaptation. Specifically, we need to estimate a scalar summary statistic, which be reasonably used to represent the overall characteristics of each of $N$ TEs at location $m$ and component $d$. To motivate this choice of summary statistic, Fig. 4 illustrates two TEs for an arbitrary choice of TC and location for the Guadeloupe application. We are interested in how the characteristics of TE change as STM increases, and therefore specify a metric which measures the relative dissimilarity of two trajectories of the temporal exposure, referred to as the trajectory distance.

For computational simplicity, we first impose a regular partition of the TE space $[0, 1]^D$, assuming $P$ cells per dimension $d$, resulting in a total of $P^D$ hypercube cells. For example, in the case of $D = 2$ with $P = 10$ illustrated in Fig. 4, the TE space is partitioned into $10^2 = 100$ cells with square shapes, indexed using $i_1$ (for $H_S$) and $i_2$ (for $U_{10}$). We then count the number of times that TE $n$ enters each of the cells, and store this as $\gamma_n^m(i_1, i_2)$. This quantity therefore summarises the distribution of TE $n$ in the $D$-dimensional TE space at location $m$. More generally, we estimate the trajectory distances

$$ \begin{gathered} \gamma_n^m(i_1, i_2, \dots, i_D) \text{ for } i_d = 1, 2, \dots, P, n = 1, 2, \dots, N, \\ m = 1, 2, \dots, M, d = 1, 2, \dots, D \end{gathered} \quad (11) $$

and use these to calculate the average trajectory distance between a TE $n$ and the remaining $N-1$ other TEs at location $m$. The average trajectory distance is calculated as

$$ \begin{gathered} \delta_n^m = \frac{1}{N} \sum_{n' \neq n} \sum_{1 \leq i_1, i_2, \dots, i_D} \left| \gamma_n^m(i_1, i_2, \dots, i_D) - \gamma_{n'}^m(i_1, i_2, \dots, i_D) \right| \\ \text{for } n = 1, 2, \dots, N, m = 1, 2, \dots, M. \end{gathered} \quad (12) $$

The average trajectory distance $\delta_n^m$ for TC $n$ at location $m$ therefore provides a univariate summary of multivariate TE characteristics, which can be used in Eq. (10) to calculate Kendall's tau statistic. In practice, we estimate the rank correlation coefficient between the average trajectory distance and each component of the MSTM per location. Results of applying this analysis for the Guadeloupe study are discussed in Section 4.4.

### 3.7. Simulation of tropical cyclones using fitted models

We use the procedure given in Algorithm 1 to generate realisations of MSTM from the set of estimated conditional models. In brief, the procedure requires us to first generate a realisation of MSTM using a combination of the marginal and conditional extreme value models. This realisation is then combined with a random drawing of a historical
<figure><img src="image_7.png" /><figcaption>Fig. 4. Left: example of a temporal exposure (TE, in red) for an arbitrary choice of TC n and location, in $e^{H_s} \cdot e^{U_{10}}$ space. The cells represent the values of $\gamma_n^m(i, j)$, with $i_1$ and $i_2$ indicating cell location in discretised $H_s$ and $U_{10}$ respectively. Right: example of TEs (in red and orange) at one location for two TCs n and n'. The cells now represent the values of $|\gamma_n^m(i, j) - \gamma_n^m(i, j)|$ used in Eq. (12). Note that the maximum of TEs (red and orange) for each of $H_s$ and $U_{10}$ correspond to exposure E in the corresponding MSTM-E model.</figcaption></figure>

E or TE to generate a realisation of a TC temporal maximum (from MSTM-E) or TC multivariate time-series (form MSTM-TE). These steps are repeated until $N^*$ realisations have been generated, sufficient to perform reliable structural risk assessment.

As mentioned in Section 3.3.2, the algorithm could also be used with modified marginal and conditional extreme models estimated using bootstrap-resampled data (see Section 3.4), as an alternative means to estimate uncertainties associated with model inferences.

## 3.8. Performance assessment

A number of different approaches can be used to assess the performance of inferences using MSTM-(T)E models. These are outlined here, and illustrated in Section 4 in application to the Guadeloupe study.

### 3.8.1. Environmental contours for MSTM and MSTM-E

For MSTM and MSTM-E models, realisations from the simulation procedure in Algorithm 1 can be used to estimate various forms of environmental contours corresponding to some long return period T years, provided that the number of realisations $N^*$ generated is large. For the Guadeloupe study, a comparison of T = 100 and T = 300-year environmental contours estimated from observed data (corresponding to a period $T_0$ = 1000 years) and from simulation (using MSTM-E models fitted to a subsample corresponding to $T_S$ = 200 years) is made. The contour estimate takes the form of the curve $(w_1(\theta), w_2(\theta))$ in $H_S - U_{10}$ space, for angle $\theta \in [0, \pi/2)$, where

$$ \mathbb{P}(H_S > w_1(\theta), U_{10} > w_2(\theta)) = \frac{1}{Tp_a p_e}. \qquad (13) $$

Here, in addition to the return period T of interest in years, $p_a$ is the annual rate of occurrence of TCs (estimated from the TC data to be 0.6). Further, parameter $p_e$ is the probability of a TC being incorporated in the sample used to estimate the contour. For contours estimated using data simulated under MSTM-E, this corresponds to the probability that the common Laplace-scale threshold $\psi$ for the conditional extremes model in Section 3 is exceeded. For contours estimated using the full original sample of data, $p_e = 1$, since all TCs are included.

In Section 4, we also assess the variability of estimates of the T = 1000 year environmental contour using MSTM-E. We achieve this by estimating MSTM-E models, and hence environmental contours, for random sub-samples (corresponding to a period $T_S$ = 200 years) of the full TC data (corresponding to a period of 1000 years). An arbitrary number of environmental contours can be estimated in this way, and

the distributional characteristics of their $(w_1(\theta), w_2(\theta))$ curve estimated as a function of $\theta \in [0, \pi/2)$. These MSTM-E contours can then be compared with a direct empirical estimate of the T = 1000-year contour obtained using the full TC data set.

### 3.8.2. Location-specific estimates for MSTM-E

It is also useful to assess estimates from MSTM-E for specific locations of interest, and to compare these with estimates obtained using location-specific estimation (LSE). The LSE estimate for location m is obtained by fitting marginal and conditional extremes models to the set of temporal maxima $\{x_{n,m}^d\}_{d=1}^D$. Since MSTM-E exploits information from all M locations in the spatial domain, we might expect estimates from MSTM-E for a specific location to be less variable than those from an LSE for that location; however, if MSTM-E modelling assumptions are not appropriate, it might be that MSTM-E estimates are also more biased that LSEs. Choice of marginal and dependence thresholds is of course also an important consideration for LSE. In the comparisons reported in Section 4, both marginal and dependence thresholds were set to the 70% percentile of the corresponding distributions.

In a similar fashion to environmental contours in Section 3.8.1, we can also estimate multiple marginal and conditional extremes models for random subsets of temporal maxima corresponding to periods of $T_S$ = 200-years at the location. In this way, as discussed in Section 4, we can quantify the variability of LSE estimates for locations of interest.

### 3.8.3. Location-specific joint distributions of TC variables for MSTM-TE

TC time-series simulated under MSTM-TE models are realisations of complex multi-dimensional objects. Assessment of the extent of agreement between historical TC time-series and those generate under MSTM-TE models is therefore also complex. In the current work, we choose to assess the performance of MSTM-TE in estimating location-specific joint distributions of observed TC variables ($H_S$ and $U_{10}$ in the case of Guadeloupe) from simulation under MSTM-TE, with empirical estimates from the full TC data set. In the comparisons in Section 4.4, we are careful to ensure a fair comparison, by only including extreme TCs (i.e. exceeding dependence thresholds with non-exceedance probability $\mathcal{L}(\psi)$ for conditional extremes in the empirical estimates from the full TC data set, where $\mathcal{L}$ is the cumulative distribution function of the standard Laplace distribution). We then compare the two estimates as contour plots exploiting Gaussian kernels for smoothing, using Scott's rule (Scott, 2015) for bandwidth selection.
**Algorithm 1:** Simulating realisations of TC temporal maxima and time-series using MSTM-(T)E

**Input:** Marginal model parameters {$u_d, \xi_d, \sigma_d$}$_{d=1}^D$ for extremes of each of *D* variables

**Input:** Conditional extremes model parameters and residuals {$\psi_d, a_{-d}, b_{-d}, R_d$}$_{d=1}^D$

**Input:** Empirically estimated rates of occurrence {$\rho_d$}$_{d=1}^D$ of MSTM over partition {$C_d$}$_{d=1}^D$ of MSTM space

**Input:** For MSTM-E, exposure set of Es {$e_{n,m}^d$}$_{n=1,m=1,d=1}^{N,M,D}$

**Input:** For MSTM-TE, exposure set of TEs {$\tilde{e}_{n,m}^d(t)$}$_{n=1,m=1,d=1,t=1}^{N,M,D,T_n}$

**STEP 1:** Simulate realisation $n^*$ of MSTM, {$s_{n^*}^d$}$_{d=1}^D$

for partition $C_{d'}$, $d' = 1, 2, ..., D$ of environmental space do

for conditioning variate $d = 1, 2, ..., D$ do

Simulate a realisation of the Laplace-scale conditioning variable *d* given rate of occurrence $\rho_d$ and dependence threshold $\psi_d$

Simulate a realisation of the conditioned variables using the conditional extremes model in Eq. (8)

Reject any realisations which generate conditioned values greater than the conditioning value

Transform realisation from Laplace to physical scale using estimated marginal cumulative distribution functions in Eq. (6)

end

end

**STEP 2:** Sample realisation $n^*$ of E or TE

For MSTM-E, sample a historical E at random, and call this {$e_{n^*,m}^d$}$_{m=1,d=1}^{M,D}$

For MSTM-TE, sample a historical TE at random, can call this {$\tilde{e}_{n^*,m}^d(t)$}$_{m=1,d=1,t=1}^{M,D,T_n^*}$

Note that the same historical (T)E is sampled for all locations and variables, preserving dependence in space and between variables

**STEP 3:** Combine MSTM and (T)E realisations (see Section 3.2)

For MSTM-E, generate a realisation of TC temporal maxima $x_{n^*,m}^d = s_{n^*,m}^d \times e_{n^*,m}^d$ for each component *d* (*d* = 1, 2, ..., *D*) at each location *m* (*m* = 1, 2, ..., *M*)

For MSTM-TE, generate a realisation of TC time-series $\tilde{x}_{n^*,m}^d(t) = s_{n^*,m}^d \times \tilde{e}_{n^*,m}^d(t)$ for each component *d* at each location *m*, for *t* = 1, 2, ..., $T_n^*$

**STEP 4:** Repeat Step 1–3 until a sufficient number $N^*$ of realisations is generated

**Output:** For MSTM-E, a set of $N^*$ realisations of TC temporal maxima {$x_{n^*,m}^d$}$_{n^*=1,m=1,d=1}^{N^*,M,D}$

**Output:** For MSTM-TE, a set of $N^*$ realisations of TC time-series {$\tilde{x}_{n^*,m}^d(t)$}$_{n^*=1,m=1,d=1,t=1}^{N^*,M,D,T_n^*}$

# 4. Application of MSTM-(T)E to Guadeloupe region

Now we apply the MSTM-(T)E methodology introduced in Section 3 to the wind and wave data introduced in Section 2 using case studies that focus on joint extremal behaviour at the five locations illustrated in Fig. 1. We estimate MSTM, MSTM-E, and MSTM-TE using subsamples of the original data set, each corresponding to a period of 200 years of observation. We then assess the performance of the estimated models in predicting return values, using the full 1000-year synthetic data set which serves as ground truth when the estimated return period is shorter than 1000 years. By repeatedly subsampling 200 years worth of data at random for estimation, we are able to quantify the bias and variance characteristics of MSTM-(T)E estimates.

The layout of the section is as follows. In Section 4.2, we characterise the joint extreme value structure for the $H_S-U_{10}$ MSTM, a common element of the MSTM-(T)E methodologies (see Fig. 2). Specifically, we discuss the estimation of marginal extreme value models for the components $H_S$ and $U_{10}$ of MSTM, and conditional extreme value models to describe the extremal dependence between $H_S$ and $U_{10}$. We emphasise the importance of sensible threshold selection for both marginal and conditional extreme value models. Environmental contours for the $H_S-U_{10}$ MSTM corresponding to a return period of 100 year are also illustrated. In Section 4.3, we discuss the estimation of the MSTM-E model, and test for dependence between MSTM and exposure. Environmental contours corresponding to 100-year and 300-year return values are estimated for each of the 5 locations given in Fig. 1. We also compare contours with those from location-specific estimation (LSE, described in Section 3.8.2). Finally, in Section 4.4, we discuss the application of MSTM-TE, comparing model-based estimates of joint evolutions of time-series for $H_S$ and $U_{10}$ with those from location-specific estimation. First however, in Section 4.1, we motivate and explain the partitioning of the spatial domain around Guadeloupe into west and east clusters, found necessary to account for different physical behaviours of extreme events in the two regions.

## 4.1. Initial analysis and partitioning of the Guadeloupe region

During preliminary data analysis, statistical testing indicated that it was not reasonable to assume that MSTM and exposure sets can be considered independent for the Guadeloupe region. Specifically, the test provided evidence for non-negligible levels of dependence between MSTM and exposure, associated with differences between TC characteristics to the west and east of Guadeloupe. To mitigate this finding, the set of TC events was partitioned into “clusters” of “east” TCs and “west” TCs based on the location of occurrence of the maximum of $H_S$ for the TC event, relative to longitude 61.5W. This partitions the set of 685 TCs into 467 TCs for the east cluster and 218 for the west cluster. MSTM-(T)E models are then developed independently for the east and west clusters. Notably, for each cluster, dependence between MSTM and exposure sets is shown to be low, as discussed in Section 4.2. The choice of longitude 61.5W as partition boundary is supported by evidence of bimodality in the distribution of MSTM components for the full Leeward region (see fig. S.2.3 in the Supplementary Material), in particular for $U_{10}$. This suggests the existence of mixed populations of TCs. We note however that this bimodal structure is not as clear when we consider TCs in the Guadeloupe region alone. There is also evidence of wave blocking and shadowing of the islands on the distribution of exposure. The method of clustering by STM location was ultimately motivated by the results of Kendall’s tau test, presented in Section 4.3.

Note that partitioning of TCs into east and west clusters is performed based on the location of occurrence of the maximum of $H_S$ for each TC event. This does not mean that TC data for a given cluster is limited spatially to the east of west subregions of the Guadeloupe region. Indeed, data from either east or west clusters typically provide a means to estimate exposure for the full Guadeloupe region, not just for the subregion corresponding to the cluster itself. As a result, we can estimate MSTM-(T)E models for the full Guadeloupe region using TC data for either east or west clusters. The performance of these models is assessed in Sections 4.3 and 4.4.

## 4.2. Estimation of MSTM

Here we assess critical features of the estimation of the $H_S-U_{10}$ MSTM model, the methodology of which is described in Section 3.3. We estimate parameters for extreme marginal and dependence models, based on results of diagnostic tests to ensure sensible marginal and dependence threshold selection. This is followed by estimation of environmental contours for MSTM (as opposed to MSTM-(T)E). Model fitting is illustrated for a typical random sample of TC data for the
<caption>**Table 1**</caption>

Bootstrap mean parameter estimates of marginal GP models, over $N_B = 1000$ bootstrap resamples, for the east and west clusters, with parameters *u* (threshold), *σ* (scale) and *ξ* (shape). Corresponding bootstrap standard deviations are given in parentheses.

<table><caption>**Table 2**</caption><thead><tr><th>Cluster</th><th>Variable</th><th>μ̂<sup>β</sup></th><th>δ̂<sup>β</sup></th><th>ξ̂<sup>β</sup></th></tr></thead><tbody><tr><td rowspan="2">East</td><td>H<sub>S</sub></td><td>10.327(0.026)</td><td>6.389(0.634)</td><td>-0.387(0.062)</td></tr><tr><td>U<sub>10</sub></td><td>33.222(0.125)</td><td>15.226(1.650)</td><td>-0.523(0.058)</td></tr><tr><td rowspan="2">West</td><td>H<sub>S</sub></td><td>6.995(0.052)</td><td>4.796(0.683)</td><td>-0.368(0.112)</td></tr><tr><td>U<sub>10</sub></td><td>26.140(0.097)</td><td>13.499(2.062)</td><td>-0.393(0.118)</td></tr></tbody></table>

Bootstrap mean parameter estimates of the conditional extremes model parameters, over $N_B = 1000$ bootstrap resamples, for the east and west clusters, with linear parameter *a*, power parameter *b* and residual mean parameter *μ*. Bootstrap standard deviations are shown in parentheses. Bootstrap mean parameter estimates are used for simulation. Standard deviation in parentheses.

<table><thead><tr><th>Cluster</th><th>Model form</th><th>δ<sup>β</sup></th><th>b<sup>β</sup></th><th>μ̂<sup>β</sup></th></tr></thead><tbody><tr><td rowspan="2">East</td><td>U<sub>10</sub><sup>L</sup>|H<sub>S</sub><sup>L</sup></td><td>0.182(0.150)</td><td>0.184(0.309)</td><td>1.044(0.411)</td></tr><tr><td>H<sub>S</sub><sup>L</sup>|U<sub>10</sub><sup>L</sup></td><td>0.024(0.250)</td><td>0.256(0.325)</td><td>1.511(0.826)</td></tr><tr><td rowspan="2">West</td><td>U<sub>10</sub><sup>L</sup>|H<sub>S</sub><sup>L</sup></td><td>0.542(0.246)</td><td>0.647(0.216)</td><td>0.513(0.311)</td></tr><tr><td>H<sub>S</sub><sup>L</sup>|U<sub>10</sub><sup>L</sup></td><td>0.601(0.219)</td><td>0.445(0.225)</td><td>0.506(0.307)</td></tr></tbody></table>

east cluster of the Guadeloupe region, corresponding to a period of observation of 200 years. We use the bootstrap scheme of Section 3.4 to quantify model uncertainty, and report estimated parameters in terms of bootstrap means and standard deviations. Then, assuming bootstrap mean model parameters for simulation, multiple model fits to differ- ent 200-year random samples are used to estimate the environmental contour for $H_S-U_{10}$ MSTM and its uncertainties. We first examine the stability of estimated models for $H_S-U_{10}$ MSTM against marginal and dependence threshold choice, and then report parameter estimates for both marginal and dependence models.

**Marginal extreme value models.** First, threshold selection for marginal extreme value models for both $H_S$ and $U_{10}$ is considered, for both west and east clusters. A non-exceedance probability of 0.60 was chosen as the marginal threshold for both clusters, corresponding to (7.00 m, 26.36 m/s) for the west cluster and (10.35 m, 33.30 m/s) for the east cluster. Details are given in fig. S.3.1 and fig. S.3.2 of section S.3.1 of the Supplementary Material. With this threshold, estimated GP parameters are summarised in Table 1. The results suggest short-tailed distributions ($\xi \ll 0$), for both $H_S$ and $U_{10}$ in both regions.

**Conditional extreme value models.** Next, we choose an appropriate threshold for conditional extremes models by assessing the stability of bootstrap estimates for model parameters *a* and *b* with respect to differ- ent combinations of marginal and dependency thresholds. Results are illustrated in fig. S.3.5, fig. S.3.6 fig. S.3.7, and fig. S.3.8 of section S.3.1 of the Supplementary Material. Dependence threshold non-exceedance probabilities of 0.70 and 0.90 were selected for west and east clusters respectively. This choice was driven in part by the lower rate of occurrence of TCs in the west cluster. For the selected dependency thresholds, samples of parameter estimates from bootstrap analysis are illustrated in fig. S.5.1 and fig. S.5.2 of section S.3.2. Note that bootstrap mean parameter estimates are used for subsequent estimation of the model residual distribution *Z* (see Eq. (8) in Section 3.3.2), and hence simulation under fitted MSTM-(T)E models. In the notation of Section 3.3.2, we also examine scatter plots of $Z_{-d}$ against large $Y_d$ to confirm no clear evidence of strong dependence; see fig. S.3.9, fig. S.3.10, and fig. S.3.11 in section S.3.2. Finally, estimated conditional model parameters are given in Table 2. It is interesting to note that estimates for linear (or slope) parameter *a* are considerably larger in the west sector that the east, but that this effect is mitigated by larger values of residual mean *μ* for the east. Estimates for power parameter *b* are larger in the west sector also.

The fitted conditional extreme value models for the east cluster are illustrated in Fig. 5, on Laplace scale (left) and physical scale (right).

In the right hand panel, the apparent “rays” or “streaks” in realisations under fitted models (shown as yellow and blue points) are the result of sampling randomly from a small residual set $\mathcal{R}$ (see Section 3.3.2), itself the result of the relative small sample of extreme TCs available for analysis. Corresponding figures for the west cluster are given in fig. S.6.1 in section S.6.

We next compare the environmental contour corresponding to a $T = 100$-year return period, estimated with bootstrap means under the MSTM model, to that estimated empirically from the full original $T_o = 1000$-year simulated TC data set. To quantify uncertainty in the location of the environmental contour in the $H_S-U_{10}$ environment space, MSTM models are estimated independently for a total of 1000 randomly-selected TC data sets, each corresponding to a period of observation $T_S = 200$ years, from the full TC data set. We then follow the procedure described in Section 3.8.1 to estimate the mean contour and its 95% uncertainty interval. Resulting environmental contours for the east cluster are provided in Fig. 6 (with the corresponding estimate for the west cluster given in fig. S.6.2 of section S.6). For both clusters, the empirical environmental contour agrees well with the mean MSTM contour, suggesting that the bias in the MSTM estimate is small. Moreover, the empirical contour falls well within the uncertainty band for the MSTM contour.

### 4.3. Estimation of MSTM-E

Using the MSTM model with bootstrap mean parameters and corre- sponding sets $\mathcal{R}$ of residuals, we now estimate $H_S-U_{10}$ MSTM-E models which characterise temporal maxima of extreme TC events in space (only). The key elements of the MSTM-E model are MSTM and (temp- oral maximum exposure) E. We proceed as follows. First, following Section 3.6 we test the assumption of independence of MSTM and E for the $H_S$ and $U_{10}$ components of the environment. Then we focus on developing estimates for environmental contours corresponding to $T = 100$-year and 300-year return periods, at the five reference locations (see Fig. 1), and comparing these with estimates obtained empirically for the full $T_o = 1000$-year data set. The procedure for simulating under the MSTM-E model is outlined in Algorithm 1.

The statistical independence of MSTM and E is quantified using rank correlation following the procedure of Section 3.6.1. For each of $H_S$ and $U_{10}$ independently, and each of the west and east sectors, we identify the spatial locations at which the rank correlation between STM and E lies beyond the central 95% band under the assumption of independence. The percentage of locations with non-negligible rank correlation is summarised in Fig. 7 for dependence thresholds $\mu$ with non-exceedance probabilities $\in [0.6, 0.9]$. For the selected dependence thresholds, corresponding to non-exceedance probabilities of 0.7 in the west cluster and 0.9 in the east cluster, the percentage of locations with non-negligible rank correlation is below 5% for both $H_S$ and $U_{10}$ in the west cluster, and for $H_S$ in the east cluster, and just slightly above 5% for $U_{10}$ in the east cluster. Spatial maps of non-negligible rank correlations for the selected dependence thresholds are given in fig. S.4.1 and fig. S.4.2 of Supplementary Material section S.4. We conclude, given wise choice of dependence thresholds, that the assumption of independence of STM and E is not obviously violated for TCs in Guadeloupe region, provided that region is considered as the union of west and east clusters of TCs with somewhat different characteristics.

We use the fitted MSTM-E models to estimate environmental con- tours at the five reference locations corresponding to return periods of $T = 100$ and $T = 300$ years. To quantify uncertainty in the location of the environmental contour, we follow the sampling proce- dure described in Section 4.2 for MSTM, using 1000 MSTM-E models from randomly-sampled data sets, each corresponding to a period of observation of $T_S = 200$ years. Also as in Section 4.2 for MSTM, we compare the MSTM-E contour with that estimated empirically from the full $T_o = 1000$-year TC data set. For further comparison, we also
<figure><img src="image_8.png" /><figcaption>**Fig. 5.** Conditional extreme values models for MSTM in east cluster. Left: conditional models on standard Laplace scale, for $U_{10}^L|H_S^L$ (yellow) and $H_S^L|U_{10}^L$ (blue). Solid lines represent models with bootstrap mean parameter estimates for *a*, *b* and *μ*. Dashed lines represent 95% uncertainty bands for predictions, estimated using the 2.5% and 97.5% quantiles of the set empirical *R* of model residuals *Z*. Right: conditional model on physical scale. Coloured points represent realisations of MSTM sampled under the fitted models for $U_{10}|H_S$ (yellow) and $H_S|U_{10}$ (blue), together with the original sample (black).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 6.** Environmental contours for $H_S-U_{10}$ MSTM of in the east cluster, corresponding to a return period of $T = 100$ years. The contour estimated under multiple MSTM models given random samples of length $T_S = 200$ years is shown in red, as a mean contour and its central 95% uncertainty band. The contour shown in black is estimated empirically using the full $T_c = 1000$ years of simulated TC data. The full underlying TC dataset for the east cluster is also shown as black dots.</figcaption></figure>

estimate a location-specific estimate (LSE) of the contour, using data only from that location following the procedure of Section 3.8.2; it was found that dependency thresholds with non-exceedance probabilities of 0.6 for the west cluster and 0.7 in the west were suitable. Uncertainty in the LSE-based contour was assessed in the same manner as that for the MSTM-E-based estimate.

Estimated environmental contours with 100-year return period, for the five reference locations, based on TCs from the east cluster only, are shown in Fig. 8 (and in ig. S.6.3 of section S.6 for the west cluster). There is good agreement between the median MSTM-E estimate (solid blue line) and the empirical ground truth estimate (solid black line) using the full $T_0 = 1000$ years of east cluster TC data for the location directly. Moreover, the empirical estimate is contained within the central 95% uncertainty interval (dashed blue line) for the MSTM-E estimate at all locations. Compared with the LSE estimates (orange), MSTM-E provides lower estimated contour bias and variance. Corresponding estimated contours for a 300-year return period are shown in Fig. 9 (and in fig. S.6.4 of section S.6 for the west cluster). The general characteristics of the figure are similar to those of Fig. 8, except perhaps that the uncertainties are greater in Fig. 9 due to the fact that we are

extrapolating further into the joint tail of $H_S$ and $U_{10}$. The LSE (orange) is particularly uncertain at the 300-year level. Again, the empirical ground truth estimate (solid black line) is contained within the central 95% uncertainty interval (dashed blue line) for the MSTM-E estimate at all locations, and the median MSTM-E (solid orange line) is in good general agreement with the empirical estimate. Note that Figs. 8 and 9 of the main text for the east cluster make use of TC data for the east cluster exclusively, whereas fig. S.6.3 and fig. S.6.4 of section S.6 make use of TC data for the west cluster. Hence, for example, empirical ground truth estimates of the 100-year environmental contour for reference location 3 are expected to be different based on east and west cluster data in Fig. 8 and fig. S.6.3 respectively.

As a further quantification of the performance of MSTM-E and LSE relative to the empirical ground truth estimate, we compare the location on each contour where the marginal quantile values are equal. This is achieved by transforming the contour to standard Laplace $H_S^L-U_{10}^L$ scale, using the appropriate marginal models, and then measuring the distance from the origin of coordinates to the point of intersection of the contour along the line $U_{10}^L = H_S^L$. This distance is referred to as the “contour distance” below, for each of discussion. The comparison is illustrated in Fig. 10 as box-whisker representations of the distributions of contour distances for MSTM-E (blue) and LSE (orange). The empirical ground truth estimate is given as a thick black horizontal line. Again, it can be seen that the bias of MSTM-E estimates is generally lower than that of LSE, and that the variance of MSTM-E estimates is considerably lower than that of LSE. We note that numerous other comparisons of this form are possible, and that the current comparison of contour distances is included as an illustration. Results of Fig. 10 are summarised in table S.1–table S.4 in section S.9 of the Supplementary Material.

The MSTM-E analysis requires the practitioner to make a number of choices of tuning parameter prior to analysis. Clearly, the user must ensure that inferences made are not overly sensitive to these choices. In section S.7, we explore the effect of reducing the size of the Guadeloupe region on the estimation of environmental contours. Our cautious conclusion is that a reasonable choice of region does not materially alter the performance of the MSTM-E model.

### 4.4. Estimation of MSTM-TE

Using the MSTM model with bootstrap mean parameters and corresponding sets $R$ of residuals, we now estimate $H_S-U_{10}$ MSTM-TE models which characterise the temporal evolution of extreme TC events in space and time. The key elements of the MSTM-TE model are
<figure><img src="image_10.png" /><figcaption>Fig. 7. Percentage of locations with non-negligible rank correlations between STM and E for $H_s$ (left) and $U_{10}$ (right), in the west cluster (triangle) and the east cluster (circle), as a function of the non-exceedance probability of the dependence threshold.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 8. Environmental contours for return period $T = 100$ years at each of reference locations 1–5, using east cluster data. In each panel, black points represent temporal maxima for TCs at the location. Coloured lines correspond to contours estimated empirically (black), using MSTM-E (blue), and using LSE (orange) respectively. Solid lines give median contours over random selections of 200 years of training data (see Section 3.8.1), and dashed lines corresponding central 95% uncertainty intervals. Empirical ground truth estimates are obtained from the full 1000 years of east cluster TC data for the location of interest.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 9. Environmental contours for return period $T = 300$ years at each of reference locations 1–5, using east cluster data. For full description, see Fig. 8.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Fig. 10.** Box-whisker plots comparing contour distances for MSTM-E (blue) and LSE (orange) with empirical ground truth estimates (thick black), for contour return periods $T = 100$ years (left) and 300 years (right) using east cluster data. The colour scheme follows Fig. 8.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Fig. 11.** Examples of TC events simulated using MSTM-TE at Location 2. Each line drawn represents the trajectory of a single event, with 3 events coloured (red, green, blue) to illustrate specific instances of the joint evolution of $H_S$ and $U_{10}$. The left and middle panels show marginal time-series plots of $H_S$ and $U_{10}$ per event, and the right hand panel shows a scatter plot of the joint time-series trajectory in $H_S - U_{10}$ space. Note that events were selected based on MSTM-E (rather than MSTM-TE) $H_S, U_{10}$ values having joint return periods in the interval of 95–105 years. Note further, in the right hand panel, that the boundary created by the set of event trajectories in $H_S - U_{10}$ space resembles the isocontour lines in Fig. 8.</figcaption></figure>

MSTM and temporal exposure TE. The TE for a TC is simply the set of time-series for $H_S$ and $U_{10}$ for all locations of interest on the spatial domain, corresponding to the time period of the TC. As can be seen from Eq. (3) in Section 3, TE is the precursor for calculation of (temporal maximum) exposure E. It is again critical to establish the approximate independence of MSTM and TE for the $H_S$ and $U_{10}$ components of the environment, in order that the MSTM-TE model can be used with confidence. This was achieved using the procedure explained in Section 3.6.2, the results illustrated in fig. S.8.1 of section S.8 of the Supplementary Material. It was found that p-values for Kendall's tau tests were ≥0.05 for less than 5% of the locations for both $H_S$ and $U_{10}$. Hence the assumption of approximate independence is deemed to be valid for each of the east and west clusters of TCs.

Examples of the TC events simulated under the MSTM-TE model at Location 2 are illustrated in Fig. 11. These examples of TC events have MSTM-E (not MSTM-TE) $H_S, U_{10}$ values with joint return period in the interval of 95–105 years. The interval is introduced to increase the population of obtaining MSTM-E samples. One hundred time series are presented, and 3 events are highlighted in colour to better illustrate the relationship between time series of $H_S$ and $U_{10}$ (right panel). The boundary created by the set of event trajectories in $H_S - U_{10}$ space resembles the isocontour lines in Fig. 8.

The MSTM-TE model combines a multivariate extreme value model for MSTM with multivariate spatial time-series model for temporal evolution of a TC over a spatial domain. Assessing the performance of a MSTM-TE model is therefore also multi-faceted. In the current work, we choose to assess the performance of the MSTM-TE in simulating time-series of sea-state $H_S$ and $U_{10}$ with realistic joint distributions at

each of the five reference locations. As explained in Section 3.8.3, we choose to do this by comparing gridded estimates of joint probability density functions of sea-state $H_S$ and $U_{10}$ estimated under the MSTM-TE model, with the ground truth estimate obtained empirically from the full $T_0 = 1000$ years of simulated TC data.

Fig. 12 illustrates the comparison for reference location 1. Comparing the left hand and centre plots, regions of high probability density are similar, at least visually. Inspection of estimates for other locations in figs. S.8.2 to S.8.5 of section S.8 confirms a similar level of agreement for the other reference locations. The right-hand panel compares contours derived from the two densities, and again suggests reasonable agreement, particularly at shorter return periods. Contours of low probability density in the empirical estimate are noisy due to small sample size; corresponding estimates from the MSTM-TE model are smoother. Similar findings have been observed for multiple other sets of subsampled cyclones. We conclude that the MSTM-TE model estimated using $T_S = 200$ years of data produces a good representation of the empirical density estimated using $T_0 = 1000$ years of data.

## 5. Discussion and conclusion

In this article, we propose a straightforward statistical model to characterise the spatio-temporal evolution of extreme tropical cyclones (TCs). The approach can be thought of as an extension of the STM-E methodology of Wada et al. (2018). STM-E assumes it is reasonable to describe the spatial statistical properties of a TC in terms of those of its spatio-temporal maxima (STM) and the spatial exposure (E) of the TC relative to its STM. Here, we extend STM-E to incorporate joint
<figure><img src="image_15.png" /><figcaption>**Fig. 12.** Joint density of sea-state $H_s$ and $U_{10}$ within an extreme tropical cyclone, at reference location 1. The left-hand panel gives the ground truth estimate obtained empirically from the full $T_0$ = 1000 years of east cluster TC data for the location. The centre panel shows the estimate obtained under multiple MSTM-TE models, each generated from (randomly-selected) $T_s$ = 200 years of east cluster TC data. The right-hand panel provides a contour plot comparing the empirical ground truth estimate (solid lines) with the MSTM-TE estimates (dashed lines).</figcaption></figure>

extremes of multiple environmental variables (like significant wave height $H_s$ and wind speed $U_{10}$) and a history matching procedure to capture time-series evolution of environmental variables within TC. Joint extremes of STMs of environmental variables are characterised using the conditional extremes model of Heffernan and Tawn (2004). Two variants of the methodology are considered: MSTM-E describes the spatial distribution of temporal maxima of TCs, ignoring time-series evolution; and MSTM-TE describes the full spatio-temporal evolution.

We demonstrate the usefulness of the MSTM-(T)E models in application to a data set of synthetic TCs in the Guadeloupe region of the Leeward Islands in the Caribbean. A set of simulated TC data corresponding to 1000 years of observations was used as ground truth to assess the performance of MSTM-(T)E models based on subsets of the full data, corresponding typically to a period of observation of 200 years.

The motivation for STM-E and MSTM-(T)E models is the provision of a pragmatic but principled approach to offshore design from small samples of metocean data. A critical assumption of all the models is that the spatial exposure *E* of a TC is independent of its STM. The models can therefore only be applied reasonably when this assumption can be justified; here, we use an approach based on Kendall's rank correlation test to quantify this. For the Guadeloupe region, we find it is necessary to develop separate models for TCs with maximum severities to the east and west of Guadeloupe.

The key results of this work are (a) demonstration that MSTM-E produces estimates of environmental design contours at a location which have low bias and variance compared to estimates made using data for just the location of interest, and (b) demonstration that MSTM-TE is able to generate spatial time-series of multiple metocean variables, the joint statistical characteristics of which reflect the ground truth.

Marginal generalised Pareto models for extremes of STMs for $H_s$ and $U_{10}$ from the synthetic TCs both indicate a negative shape parameter, and hence a finite upper end point to the marginal distribution. This is particularly clear in the data for $U_{10}$, and is reflected in simulations under MSTM-(T)E models. Conditional extremes models for $H_s$ given large values of $U_{10}$ (and vice versa) indicate positive dependence between the variables (with model slope parameter *a* estimated to be >0). However, we also find *a* < 1 in all cases, suggesting asymptotic independence between $H_s$ and $U_{10}$; that is, the very largest values of one metocean variable do not typically coincide with the very largest values of the other. This may align with other studies into the effect of *extended fetch* in hurricanes (Young, 2003) or the saturation of drag coefficient under high winds (Donelan, 2018; Powell et al., 2003; Takagaki et al., 2016), all of which suggest is a non-linear relationship between extreme wind and wave, especially in cyclonic conditions. More discussion on the conditional model parameters can be found in section S.5 of the Supplementary Material. It is also interesting

that different marginal threshold choices were found appropriate for inference on east and west cluster TCs, but this can be explained at least in part by the different rate of occurrence of TCs to the east and west of Guadeloupe.

There are many aspects of the MSTM-(T)E that require further study, but we think that we have already established that the methodologies are useful. For the Guadeloupe region, the prior partitioning of TCs to those with maximum severities to the east and west of Guadeloupe suggests the possibility of adopting a mixture model for MSTM-(T)E, within which the effects of TCs from both the east and west clusters of TCs would be captured. The specific partitioning employed, using longitude 61.5W as boundary, is motivated by physical intuition. However, we do not claim that this partitioning is optimal, but it is simple. The conditional extremes model has proved a useful approach to quantifying the extremal dependence between metocean variables. It would be interesting to consider to what extent the characteristics of multivariate extreme events simulated under conditional extremes models preserve known or anticipated physical constraints.

MSTM-TE generates coupled time-series of multiple metocean variables across a spatial domain, the statistical characteristics of which are complex. To date, we have demonstrated that the approach is capable of capturing the joint distribution of sea-state $H_s$ and $U_{10}$ for an extreme tropical TC at a location. There are many other diagnostic analyses which might supplement this, in order to more fully assess the relative characteristics of MSTM-TE estimates of the ground truth. Specifically, it would be interesting to summarise the full multivariate spatio-temporal MSTM-TE output in terms of one or two summary "structure variables", such as an approximate environmental $H_s-U_{10}$ load on an offshore structure (Masoomi et al., 2019). Practically, better quantification of bias and the variance of estimates for marginal and joint extremes corresponding to longer return periods from limited data is necessary to assess the usefulness of MSTM-TE. We can generate multiple MSTM-TE models from random 200-year subsamples of the full TC data, as we did in MSTM and MSTM-E validation. However, the full 1000 years of synthetic TC data is insufficient to characterise the ground truth in this case. Given the limitation of 1000 years of ground truth data, one possible study would be to assess the usefulness of TC models based on samples corresponding to much shorter periods of time (e.g. 20 years) in estimating extreme sets with return periods of around 100 years.

### CRediT authorship contribution statement

**Kosuke Sando:** Writing – original draft, Methodology, Investigation, Formal analysis. **Ryota Wada:** Writing – review & editing, Writing – original draft, Validation, Supervision, Project administration, Methodology, Investigation, Funding acquisition, Formal analysis.
**Jérémy Rohmer:** Writing – review & editing, Investigation, Data curation. **Philip Jonathan:** Writing – review & editing, Validation, Supervision, Methodology, Investigation, Conceptualization.

Declaration of competing interest

The authors declare that they have no known competing finan-
cial interests or personal relationships that could have appeared to
influence the work reported in this paper.

Data availability

The authors do not have permission to share data.

Acknowledgements

We acknowledge the contribution of Rodrigo Pedreros, Andrea Fil-
ippini and Sophie Lecacheux from BRGM (French Geological Survey)
for providing the data set of TCs passing through Guadeloupe used
in this research. The second co-author acknowledges funding of The
Japan Society of Naval Architects and Ocean Engineers (JASNAOE).
The third co-author acknowledges funding of the Carib-Coast Inter-
reg project (https://www.interreg-caraibes.com/carib-coast, grant no.
2014TC16RFTN008).

**Appendix A. Supplementary data**

Supplementary material related to this article can be found online
at https://doi.org/10.1016/j.oceaneng.2024.118365.

References

Anderson, D., Rueda, A., Cagigal, L., Antolinez, J., Mendez, F., Ruggiero, P., 2019. Time-varying emulator for short and long-term analysis of coastal flood hazard potential. J. Geophys. Res.: Oceans 124 (12), 9209–9234.

Arduin, F., Rogers, E., Babanin, A.V., Filipot, J.-F., Magne, R., Roland, A., Van Der Westhuysen, A., Queffeulou, P., Lefevre, J.-M., Aouf, L., et al., 2010. Semiem- pirical dissipation source functions for ocean waves. Part I: Definition, calibration, and validation. J. Phys. Oceanogr. 40 (9), 1917–1941.

Bloemendaal, N., Haigh, I.D., de Moel, H., Muis, S., Haarsma, R.J., Aerts, J.C., 2020. Generation of a global synthetic tropical cyclone hazard dataset using STORM. Sci. Data 7 (1), 1–12.

Brown, B.G., Katz, R.W., Murphy, A.H., 1984. Time series models to simulate and forecast wind speed and wind power. J. Appl. Meteorol. Climatol. 23 (8), 1184–1195.

Caires, S., Sterl, A., 2005. 100-Year return value estimates for ocean wind speed and significant wave height from the ERA-40 data. J. Clim. 18 (7), 1032–1048.

Cardone, V.J., Cox, A.T., 2013. Tropical cyclone marine surface wind modeling: the shape of the radial wind profile matters. In: Advances in Hurricane Engineering: Learning from Our Past. ASCE, pp. 941–959.

Chavez-Demoulin, V., Davison, A., 2012. Modelling time series extremes. Revstat 10, 109–133.

Coles, S., 2001. An Introduction to Statistical Modeling of Extreme Values. In: Springer Series in Statistics, Springer London, London.

Davison, A.C., Padoan, S.A., Ribatet, M., 2012. Statistical modelling of spatial extremes. Statist. Sci. 27, 161–186.

Donelan, M.A., 2018. On the decrease of the oceanic drag coefficient in high winds. J. Geophys. Res.: Oceans 123 (2), 1485–1501.

Dong, S., Wang, N., Liu, W., Soares, C.G., 2013. Bivariate maximum entropy distribution of significant wave height and peak period. Ocean Eng. 59, 86–99.

Dufoy, A., 2021. Multivariate extreme value theory: practice and limits. In: Bousquet, N., Bernardara, P. (Eds.), Extreme Value Theory with Applications to Natural Hazards: From Statistical Theory to Industrial Practice. Springer International Publishing, pp. 191–228.

Ferreira, J., Soares, C.G., 2002. Modelling bivariate distributions of significant wave height and mean wave period. Appl. Ocean Res. 24 (1), 31–45.

Hansen, H.F., Randell, D., Zeeberg, A.R., Jonathan, P., 2020. Directional-seasonal extreme value analysis of north sea storm conditions. Ocean Eng. 195, 106665.

Haver, S., 1987. On the joint distribution of heights and periods of sea waves. Ocean Eng. 14 (5), 359–376.

Heffernan, J.E., Tawn, J.A., 2004. A conditional approach for multivariate extreme values (with discussion). J. R. Stat. Soc. Ser. B Stat. Methodol. 66 (3), 497–546.

Heideman, J.C., Mitchell, D.A., 2009. Grid point pooling in extreme value analysis of hurricane hindcast data. J. Waterw. Port Coast. Ocean Eng. 135 (2), 31–37.

Holland, G.J., 1980. An analytic model of the wind and pressure profiles in hurricanes. Mon. Weather Rev. 108 (8), 1212–1218.

Huser, R., Opitz, T., Wadsworth, J., 2024. Modeling of spatial extremes in envi-
ronmental data science: time to move away from max-stable processes. arXiv:
2401.17430.

Huser, R., Wadsworth, J.L., 2022. Advances in statistical modeling of spatial extremes. WIREs Comput. Stat. 14, e1537.

Jonathan, P., Ewans, K., 2007. Uncertainties in extreme wave height estimates for hurricane-dominated regions. J. Offshore Mech. Arct. Eng. 129 (4), 300–305. http://dx.doi.org/10.1151/1.2746401.

Jonathan, P., Ewans, K., 2013. Statistical modelling of extreme ocean environments for marine design: a review. Ocean Eng. 62, 91–109.

Landsea, C.W., Franklin, J.L., 2013. Atlantic hurricane database uncertainty and presentation of a new database format. Mon. Weather Rev. 141 (10), 3576–3592.

Mackay, E., de Hauteclocque, G., Vanem, E., Jonathan, P., 2021. The effect of serial correlation in environmental conditions on estimates of extreme events. Ocean Eng. 242, 110092.

Masoomi, H., van de Lindt, J.W., Ameri, M.R., Do, T.Q., Webb, B.M., 2019. Combined wind-wave-surge hurricane-induced damage prediction for buildings. J. Struct. Eng. 145 (1), 04018227.

Powell, M.D., Vickery, P.J., Reinhold, T.A., 2003. Reduced drag coefficient for high wind speeds in tropical cyclones. Nature 422 (6929), 279–283.

Ross, E., Randell, D., Ewans, K., Feld, G., 2017. Efficient estimation of return value distributions from non-stationary marginal extreme value models using Bayesian inference. Ocean Eng. 142, 315–328.

Scott, D.W., 2015. Multivariate Density Estimation: Theory, Practice, and Visualization. John Wiley & Sons.

Takagaki, N., Komori, S., Suzuki, N., Iwano, K., Kurose, R., 2016. Mechanism of drag coefficient saturation at strong wind speeds. Geophys. Res. Lett. 43 (18), 9829–9835.

Tawn, J., Shooter, R., Towe, R., Lamb, R., 2018. Modelling spatial extreme events with environmental applications. Spatial Stat. 28, 39–58.

Teena, N., Sanil Kumar, V., Sudheesh, K., Sajeev, R., 2012. Statistical analysis on extreme wave height. Nat. Hazards 64, 223–236.

Tendijck, S., Jonathan, P., Randell, D., Tawn, J., 2023a. Temporal evolution of the extreme excursions of multivariate kth order Markov processes with application to oceanographic data. Environmetrics e2834.

Tendijck, S., Tawn, J.A., Jonathan, P., 2023b. Extremal characteristics of conditional models. Extremes 26, 139–156.

Tolman, H., 2014. The WAVEWATCH III Development Group: User Manual and System Documentation of WAVEWATCH III Version 4.18. Technical Note, Environmental Modeling Center, National Centers for Environmental Prediction, National Weather Service, National Oceanic and Atmospheric Administration, US Department of Commerce, College Park, MD.

Vickery, P., Skerlj, P., Twisdale, L., 2000. Simulation of hurricane risk in the US using empirical track model. J. Struct. Eng. 126 (10), 1222–1237.

Wada, R., Rohmer, J., Krien, Y., Jonathan, P., 2022. Statistical estimation of spatial wave extremes for tropical cyclones from small data samples: validation of the STM-E approach using long-term synthetic cyclone data for the Caribbean sea. Nat. Hazards Earth Syst. Sci. 22 (2), 431–444.

Wada, R., Waseda, T., 2020. Assessment of data-inherited uncertainty in extreme wave analysis. J. Offshore Mech. Arct. Eng. 142 (2), 021204.

Wada, R., Waseda, T., Jonathan, P., 2016. Extreme value estimation using the likelihood-weighted method. Ocean Eng. 124, 241–251.

Wada, R., Waseda, T., Jonathan, P., 2018. A simple spatial model for extreme tropical cyclone seas. Ocean Eng. 169, 315–325.

Walshaw, D., 2023. Multivariate extreme value theory. In: Wiley StatsRef: Statistics Reference Online. John Wiley and Sons, Ltd.

Young, I.R., 2003. A review of the sea state generated by hurricanes. Mar. Struct. 16 (3), 201–218.