# Estimating Metocean Environments Associated with Extreme Structural Response

Matthew Speers  David Randell  Jonathan Tawn  Philip Jonathan  Department of
Mathematics and Statistics, Lancaster University LA1 4YF, United Kingdom.
Shell Global Solutions International BV, Amsterdam, The Netherlands.  Shell
Research Limited, London, United Kingdom.

###### Abstract

Extreme value analysis (EVA) uses data to estimate long-term extreme
environmental conditions for variables such as significant wave height and
period, for the design of marine structures. Together with models for the
short-term evolution of the ocean environment and for wave-structure
interaction, EVA provides a basis for full probabilistic design analysis.
Environmental contours provide an alternate approach to estimating structural
integrity, without requiring structural knowledge. These contour methods also
exploit statistical models, including EVA, but avoid the need for structural
modelling by making what are believed to be conservative assumptions about the
shape of the structural failure boundary in the environment space. These
assumptions, however, may not always be appropriate, or may lead to
unnecessary wasted resources from over design. We introduce a methodology for
full probabilistic analysis to estimate the joint probability density of the
environment, conditional on the occurrence of an extreme structural response,
for simple structures. We use this conditional density of the environment as a
basis to assess the performance of different environmental contour methods. We
demonstrate the difficulty of estimating the contour boundary in the
environment space for typical data samples, as well as the dependence of the
performance of the environmental contour on the structure being considered.

###### keywords:

Structural design, extreme, full probabilistic analysis, contour, IFORM,
conditional simulation, importance sampling, significant wave height, wave
steepness.

###### PACS:

0000 , 1111

###### MSC:

0000 , 1111

††journal: Ocean Engineering

##  1 Introduction

Ocean engineers use different approaches to quantify extreme conditions for
design and reassessment of offshore and coastal structures. The natural full
probabilistic approach (henceforth, the “forward” approach) is to construct a
sequence of statistical models to characterise the extreme multivariate ocean
environment, as well as the interaction between that environment and the
structure (e.g., Towe et al. 2021). This approach considers the response of
the structure to be a stochastic function of the environment, summarised by
underlying sea state statistics such as significant wave height and period,
contrary to previous work (e.g., Coles and Tawn 1994) where a deterministic
relationship is assumed. The forward approach thus seeks a multivariate
distribution $F_{\mathbf{X}}$ for environmental variables $\mathbf{X}$, as
well as a distribution $F_{R_{L}|\mathbf{X}}$ which characterises the maximum
stochastic response $R_{L}$ induced on the structure by the environment
$\mathbf{X}$ over the period of a sea state of duration $L$ (e.g., 3) hours. A
key property of this method is that uncertainty from the estimation of each
distribution can be naturally quantified and propagated through the sequence
of models. Structural risk assessment centres on the estimation of the
probability of structural failure using the distribution $F_{R_{S}}$ for the
maximum response $R_{S}$ in a random storm, or alternatively the distribution
$F_{R_{A}}$ for the maximum response $R_{A}$ per annum, evaluated by
marginalisation of the distribution $F_{R_{L}|\mathbf{X}}$ over the
environment space. In addition, $F_{R_{L}|\mathbf{X}}$ may be used to find the
probability of the event of structural failure $\mathcal{F}_{L}$ within a sea
state (of duration $L$) with environmental conditions $\mathbf{X}$. The first
objective of the current work is to demonstrate the performance of an
efficient forward model to estimate the tail of the distributions
$F_{R_{L}|\mathbf{X}}$ and $F_{R_{S}}$, corresponding to return periods $P$ of
the order of $10^{3}$ years.

A combination of better models for the extreme ocean environment, techniques
to reduce the computational complexity of the forward approach, and improved
computational resources, have made forward estimation of structural failure
probability more routinely achievable. These include the development of
practically-useful statistical models for non-stationary, or covariate
effected, margins (e.g., Chavez-Demoulin and Davison 2005; Randell et al.
2016; Youngman 2019) as well as the conditional multivariate extremes model
developed by Heffernan and Tawn (2004), conditional simulation of extreme time
series (of waves and wave kinematics) proposed by Taylor et al. (1997), and
efficient importance sampling from distributions (e.g., Gelman et al. 2013).
We demonstrate the direct estimation of the distribution of base shear,
approximated using the loading equation of Morison et al. (1950) on simple
structures, following the procedure of Tromans and Vanderschuren (1995). The
environmental and structural models utilised here are sufficiently complex to
illustrate key methodological steps, whilst being simple enough to avoid
unnecessary complexity. We will also estimate the joint conditional density of
environmental variables (henceforth referred to as $\text{CDE}_{P}$) given the
occurrence of an extreme $P$-year Morison load on example structures.

Historically however, forward estimation of the probability of failure has
proved computationally intractable or prohibitively expensive. Instead,
metocean design has tended to focus on a dominant variable (such as
significant wave height) at a location, placing less emphasis on other
associated environmental variables (e.g., Feld et al. 2015); the metocean
engineer’s challenge is then to estimate marginal return values of the
dominant variable, and perhaps follow an engineering recipe to specify design
values for associated variables. The development of methods of structural
reliability, associated with Madsen et al. (1986), made it clearer that good
models for the joint distribution of environmental variables were necessary.
For this reason, approaches to structural reliability which make use of more
than one environmental variable became popular, such as environmental design
contours derived from parametric hierarchical models of the environmental
variables $\mathbf{X}$ (e.g., Haver 1987). A recent review of statistical
methodologies for metocean design is provided by Vanem et al. (2022).

In one respect, contour methods are advantageous over the forward approach in
that they characterise the environment only and so just require estimation of
the joint environmental distribution $F_{\mathbf{X}}$. Therefore, combined
with appropriate assumptions for the nature of the ocean-structure
interaction, environmental contours can be used in principle to assess any
structure in that environment. The decoupling of environment and structure is
achieved by making what are thought to be conservative assumptions about the
environment-structure interaction, leading to what are believed to be
conservative estimates of structural reliability. Many different methods exist
to estimate environmental design contours (e.g., Ross et al. 2020;
Haselsteiner et al. 2021; Mackay and Haselsteiner 2021; Hafver et al. 2022;
Mackay and de Hauteclocque 2023). These methods, whilst each attempting to
construct the $P$-year level contour, can produce quite different estimates to
one another. In addition, the environmental contour produced by a particular
contour estimation method will vary with the estimate of $F_{\mathbf{X}}$ used
to obtain it. In this work, we consider the IFORM environmental contour method
(e.g., Haver and Winterstein 2009), recommended by both NORSOK N-003 (2017)
and DNVGL-RP-C205 (2017) standards. We will assess the relative predictive
performance of different hierarchical models for $F_{\mathbf{X}}$ and their
respective IFORM environmental contours. The second objective of this work is
to compare estimates for the $\text{CDE}_{P}$ of metocean variables with
different estimators of IFORM environmental contours for the same return
period $P$, for a set of simple structures.

The layout of the article is as follows. In Section 2, we seek to motivate our
analysis using a sample of data for storm peak significant wave height and
second spectral moment wave period from a location in the central North Sea.
Section 3 describes the methodologies combined to achieve the forward model
for estimation of the distributions $F_{R_{S}}$ and $F_{R_{A}}$. The approach
used to estimate environmental design contours, and the various parametric
forms considered for the hierarchical estimation of $F_{\mathbf{X}}$ are given
in Section 4. In Section 5, we present estimates of $\text{CDE}_{P}$ for three
variants of a simple stick structure, and use these to quantify the
performance of different IFORM environmental contours. In Section 6, we
discuss the implications of our results, and make recommendations for
structural design practice. We provide an online Supplementary Material (SM)
with a fuller description of aspects of the procedures above.

##  2 Motivating metocean dataset

We motivate the analysis using hindcast data for sea state significant wave
height and second spectral moment wave period for a location in the central
North Sea. The data consist of 124671 observations for the period January 1979
to September 2013, calculated for consecutive 3-hour sea states. Intervals
corresponding to storm events are isolated from the hindcast data, using the
approach of Ewans and Jonathan (2008), resulting in a total of 2462 values for
storm peak significant wave height ($H_{S}$) and corresponding wave period
($T_{2}$), for an average of $73$ storm events per annum. Figure 1(a) shows
the storm peak data $(H_{S},T_{2})$. Despite the fact that we expect these
variables not to be identically distributed due to environmental covariates
(e.g., direction and season), for the purposes of the current work we assume
these to be independently and identically distributed.

(a) Storm peak variables $(H_{S},T_{2})$.

(b) Storm peak variables $(H_{S},S_{2})$.

Figure 1: Data of storm peak sea state variables from a location in the
central North Sea. Storm peaks are extracted from hindcast data of consecutive
$3$-hour sea states, using the method of Ewans and Jonathan (2008), in terms
of significant wave height $H_{S}$ [m], second spectral moment wave period
$T_{2}$ [s] and wave steepness $S_{2}$.

We choose to use storm peak wave steepness $S_{2}$ in favour of wave period
$T_{2}$ in the analysis below. Values of $S_{2}$ are calculated via
$S_{2}=2\pi H_{S}/(gT_{2}^{2})$, where $g$ is the acceleration due to gravity.
Figure 1(b) shows the resulting storm peak data $\mathbf{X}=(H_{S},S_{2})$.
Note throughout the paper, that we restrict the notation $H_{S}$, $S_{2}$ and
$T_{2}$ to refer to storm peak quantities, and further that all numerical
values are quoted in SI units. Additionally, all physical properties are taken
to be real-valued unless stated otherwise.

The choice of $(H_{S},S_{2})$ over $(H_{S},T_{2})$ is motivated by the fact
that extreme value models are generally developed for the joint upper tail of
variables. Extreme environmental loads are typically generated by large values
of $H_{S}$, large values of $S_{2}$ but non-extreme values of $T_{2}$; that
is, joint extremes of $(H_{S},S_{2})$ in the upper tail induce the most
extreme environments and structural responses. Therefore, it is appealing to
structure the analysis in terms of $(H_{S},S_{2})$ (e.g., Myrhaug 2018).

##  3 Environment and response modelling

###  3.1 Outline of the forward model and $\text{CDE}_{P}$

####  3.1.1 Storm peak characteristics and intra-storm evolution

We illustrate the methodologies combined to form the forward approach for
direct estimation of the distribution of extreme structural response
(specifically, base shear), in order to estimate the conditional distribution
of the environment $\text{CDE}_{P}$. For generality, and alignment with the
work of others (e.g., Towe et al. 2021), initially we present a form of the
forward model which incorporates full intra-storm evolution and the effects of
covariates. We subsequently restrict the model, and focus our analysis on
storm peak variables as introduced in Section 2, so as to emphasise the key
methodological steps only.

We assume access to metocean data for storm peak variables
$\mathbf{X}^{\text{sp}}$ (e.g., $H_{S}$ and $S_{2}$ from Section 2). Since the
characteristics of storm peak events generally vary with respect to covariates
(e.g., Randell et al. 2015) we also assume access to storm peak covariates
$\boldsymbol{\Theta}^{\text{sp}}$, which encode this information. The response
of an offshore structure to environmental loading occurs continuously, and in
particular for the full duration of the storm event, consisting of a series of
sea states (of duration $L$ hours) indexed by
$s\in\mathcal{S}_{T}=\\{1,2,\ldots,T\\}$, for unknown storm length $T$. To
estimate the distribution of the maximum structural response $R_{S}$ in a
storm, we need to consider (a) the variability of the duration $T$ of a random
storm, (b) the full evolution of sea state variables
$\\{\mathbf{X}_{s}\\}_{s\in\mathcal{S}_{T}}$ over a given storm (themselves
not identically distributed given the sea state covariates
$\\{\boldsymbol{\Theta}_{s}\\}_{s\in\mathcal{S}_{T}}$) and (c) the resulting
maximum response $R_{S}|\\{(\mathbf{X}_{s},\boldsymbol{\Theta}_{s})\\}_{s\in
S_{T}}$ induced over sea states within the storm.

As sea states are dependent over time, therefore $R_{L,s}$ and
$R_{L,s^{\prime}}$, for any $s,s^{\prime}\in\mathcal{S}_{T}$ where $s\neq
s^{\prime}$, will also be dependent, where $R_{L,s}$ is the maximum response
in a sea state of length $L$ at time index $s$. However, this dependence is
purely due to the sea state characteristics
$\\{(\mathbf{X}_{s},\boldsymbol{\Theta}_{s})\\}_{s\in\mathcal{S}_{T}}$
evolving over time. This is because a load $R_{L,s}$ from the model of Morison
et al. (1950) is produced by an individual random wave event, and the
statistical properties of wave events within a sea state at time index $s$ are
determined by the characteristics $(\mathbf{X}_{s},\boldsymbol{\Theta}_{s})$,
yet the interval of time over which consecutive waves are correlated is
considerably shorter than the length $L$ of a sea state. Therefore, the
conditional dependence between the variables $R_{L,s}$ and $R_{L,s^{\prime}}$,
given $(\mathbf{X}_{s},\boldsymbol{\Theta}_{s})$ and
$(\mathbf{X}_{s^{\prime}},\boldsymbol{\Theta}_{s^{\prime}})$, is negligible.
We exploit this reasoning to facilitate step (c), where we assume that the
random variables $\\{R_{L,{s}}\\}_{s\in\mathcal{S}_{T}}$ are conditionally
independent given sea states
$\\{(\mathbf{X}_{s},\boldsymbol{\Theta}_{s})\\}_{s\in\mathcal{S}_{T}}$, which
leads to the following simplification, for response $r>0$

$\displaystyle F_{R_{S}|\\{(\mathbf{X}_{s},\boldsymbol{\Theta}_{s})\\}_{s\in
S_{T}}}(r|\\{(\mathbf{x}_{s},\boldsymbol{\theta}_{s})\\}_{s\in
S_{\tau}})=\prod_{s\in\mathcal{S}_{T}}F_{R_{L}|(\mathbf{X}_{s},\boldsymbol{\Theta}_{s})}(r|\mathbf{x}_{s},\boldsymbol{\theta}_{s}),$
(1)

where we omit the subscript $s$ in writing
$F_{R_{L}|(\mathbf{X}_{s},\boldsymbol{\Theta}_{s})}$ since the dependence of
$R_{L}$ on $s$ is contained through the vales of
$(\mathbf{X}_{s},\boldsymbol{\Theta}_{s})$.

####  3.1.2 Distribution of maximum response per storm and per annum

The forward model estimates the cumulative distribution function $F_{R_{S}}$
of the maximum structural response $R_{S}$ in a storm, and subsequently the
distribution of maximum response $R_{A}$ per annum. Estimation of $F_{R_{S}}$
requires (a) modelling of multivariate storm peak variables
$\mathbf{X}^{\text{sp}}$ given storm peak covariates
$\boldsymbol{\Theta}^{\text{sp}}$, (b) characterisation of the conditional
time-varying within-storm evolution of sea state characteristics
$\\{(\mathbf{X}_{s},\boldsymbol{\Theta}_{s})\\}_{s\in\mathcal{S}_{T}}$ given
storm peak characteristics
$\mathbf{X}^{\text{sp}},\boldsymbol{\Theta}^{\text{sp}}$, and (c) the
estimation of the maximum response $R_{L,s}$ given sea state characteristics
$(\mathbf{X}_{s},\boldsymbol{\Theta}_{s})$ for an $L$-hour sea state. We must
also consider the variability in the duration $T$ of a storm event. See Towe
et al. (2021) for previous discussions of similar models. Given knowledge of
the above, and exploiting assumption (1) made in Section 3.1, the distribution
of $R_{S}$ can be written as

$\displaystyle F_{R_{S}}(r)$
$\displaystyle=\int_{\boldsymbol{\theta}^{\text{sp}}}\int_{\mathbf{x}^{\text{sp}}}\int_{(\left\\{(\mathbf{x}_{s},\boldsymbol{\theta}_{s}\right)\\}_{s\in\mathcal{S}_{\tau}},\tau)}\prod_{s\in\mathcal{S}_{\tau}}F_{R_{L}|(\mathbf{X}_{s},\boldsymbol{\Theta}_{s})}(r|\mathbf{x}_{s},\boldsymbol{\theta}_{s})$
(2) $\displaystyle\times
f_{(\left\\{(\mathbf{X}_{s},\boldsymbol{\Theta}_{s})\right\\}_{s\in\mathcal{S}_{T}},T)|\mathbf{X}^{\text{sp}},\boldsymbol{\Theta}^{\text{sp}}}\left(\left\\{(\mathbf{x}_{s},\boldsymbol{\theta}_{s})\right\\}_{s\in\mathcal{S}_{\tau}},\tau\mid\mathbf{x}^{\text{sp}},\boldsymbol{\theta}^{\text{sp}}\right)$
(3) $\displaystyle\times
f_{\mathbf{X}^{\text{sp}}|\boldsymbol{\Theta}^{\text{sp}}}(\mathbf{x}^{\text{sp}}|\boldsymbol{\theta}^{\text{sp}})\times
f_{\mathbf{\Theta}^{\text{sp}}}(\boldsymbol{\theta}^{\text{sp}})\,\,\text{d}(\\{(\mathbf{x}_{s},\boldsymbol{\theta}_{s})\\}_{s\in\mathcal{S}_{\tau}},\tau)\,\,\text{d}\mathbf{x}^{\text{sp}}\,\,\text{d}\boldsymbol{\theta}^{\text{sp}},$
(4)

for $r>0$, where $f_{\boldsymbol{\Theta}^{\text{sp}}}$ is the joint
probability density of storm peak covariates,
$f_{\mathbf{X}^{\text{sp}}|\boldsymbol{\Theta}^{\text{sp}}}$ is the joint
probability density of storm peak variables given storm peak covariates, and
$f_{(\left\\{(\mathbf{X}_{s},\boldsymbol{\Theta}_{s})\right\\}_{s\in\mathcal{S}_{T}},T)|\mathbf{X}^{\text{sp}},\boldsymbol{\Theta}^{\text{sp}}}$
is the joint probability density for the full time-series evolution of within-
storm sea state characteristics (and storm duration) given storm peak
characteristics.

If we assume that the number of occurrences of storm events in a year is
Poisson-distributed with expectation $\lambda$ per annum, we can use
$F_{R_{S}}$ to estimate the corresponding cumulative distribution function
$F_{R_{A}}$ of the maximum response $R_{A}$ in a year, i.e.,

$F_{R_{A}}(r)=\sum_{m=0}^{\infty}[F_{R_{S}}(r)]^{m}\frac{\lambda^{m}e^{-\lambda}}{m!}=\exp\left[-\lambda\left(1-F_{R_{S}}(r)\right)\right].$
(5)

for $r>0$. From this expression, we may define return values for maximum
response corresponding to a return period of $P$ years as
$r_{P}=F_{R_{A}}^{-1}(1-1/P)$.

####  3.1.3 Reduced forward model

The objective of the current work is to compare estimates for $\text{CDE}_{P}$
made using the forward approach, summarised by (4), with environmental
contours at the $P$-year level. This comparison is useful even in the absence
of covariate effects, and whilst neglecting intra-storm evolution. Therefore,
to minimise computational complexity, we now assume that the effects of
covariates $\boldsymbol{\Theta}^{\text{sp}}$ and $\boldsymbol{\Theta}_{s}$ can
be ignored, and that the maximum response in a storm always occurs in the
storm peak sea state, so that intra-storm evolution can also be ignored. As a
result, integral (4) for $F_{R_{S}}$ reduces to

$F_{R_{S}}(r)=\int_{\mathbf{x}^{\text{sp}}}F_{R_{L}|\mathbf{X}^{\text{sp}}}(r|\mathbf{x}^{\text{sp}})f_{\mathbf{X}^{\text{sp}}}(\mathbf{x}^{\text{sp}})\text{d}\mathbf{x}^{\text{sp}},$
(6)

for $r>0$, where $F_{R_{L|\mathbf{X}^{\text{sp}}}}$ is the cumulative
distribution function of the maximum response over an $L$-hour sea state given
storm peak variables $\mathbf{X}^{\text{sp}}$. For brevity, we henceforth omit
the storm peak superscript and write

$F_{R_{S}}(r)=\int_{\mathbf{x}}F_{R_{L}|\mathbf{X}}(r|\mathbf{x})f_{\mathbf{X}}(\mathbf{x})\text{d}\mathbf{x},$
(7)

for $r>0$, where $f_{\mathbf{X}}$ is now the joint density of storm peak
variables $\mathbf{X}=(X_{1},\ldots,X_{p})$, and $F_{R_{L}|\mathbf{X}}$ is the
distribution of maximum response over an $L$-hour storm peak sea state with
variables $\mathbf{X}$. This brings us back to the setting of Section 2, where
we introduced observations of storm peak sea state data $\mathbf{X}$ with
$p=2$. The key steps in evaluating the reduced forward model in (7) therefore
become the estimation of $f_{\mathbf{X}}$ and $F_{R_{L}|\mathbf{X}}$. The
first of these inferences is achieved using the conditional extremes model of
Heffernan and Tawn (2004), as described in Section 3.2. The second inference
involves conditional simulation of environmental time-series following Taylor
et al. (1997), and importance sampling from the tail of $F_{R_{L}|\mathbf{X}}$
following Towe et al. (2021), as described in Section 3.3.

####  3.1.4 Probability of structural failure

When designing offshore structures, it is often desirable to determine the
probability of the event of structural failure $\mathcal{F}_{L}$ for a given
$L$-hour sea state with variables $\mathbf{X}=\mathbf{x}$. Given an estimate
for the distribution of $R_{L}|(\mathbf{X}=\mathbf{x})$, we can evaluate this
probability using the expression

$\Pr\left(\mathcal{F}_{L}|\mathbf{X}=\mathbf{x}\right)=\int_{r_{L}}\Pr\left(\mathcal{F}_{L}|R_{L}=r_{L}\right)f_{R_{L}|\mathbf{X}}(r_{L}|\mathbf{x})\,\text{d}r_{L},$
(8)

where the failure probability $\Pr\left(\mathcal{F}_{L}|R_{L}=r_{L}\right)$
depends on the nature of the structure, as illustrated in Figure 6 in Section
5.

####  3.1.5 Conditional density of the environment

The conditional density of the environment $f_{\mathbf{X}|R_{L}}(\cdot|r_{P})$
describes the joint density of the environmental variables $\mathbf{X}$,
conditional on the appearance of a $P$-year maximum response $r_{P}$ within a
sea state of length $L$ hours. Given estimates for $f_{\mathbf{X}}$,
$F_{R_{L}|\mathbf{X}}$, $F_{R_{L}}$ and $r_{P}$, $\text{CDE}_{P}$ can
therefore be evaluated using Bayes’ rule

$f_{\mathbf{X}|R_{L}}(\mathbf{x}|r_{P})=\frac{f_{R_{L}|\mathbf{X}}(r_{P}|\mathbf{x})f_{\mathbf{X}}(\mathbf{x})}{f_{R_{L}}(r_{P})},$
(9)

for $\mathbf{x}\in\mathbb{R}^{p}$, where $f_{R_{L}|\mathbf{X}}$ and
$f_{R_{L}}$ are the densities corresponding to the distributions
$F_{R_{L}|\mathbf{X}}$ and $F_{R_{L}}$ respectively. Examples of
$\text{CDE}_{P}$ for the central North Sea application are given in Section
5.3.

###  3.2 Joint modelling of storm peak conditions

####  3.2.1 Outline of the conditional extremes model

The upper extremes of the marginal and joint distributions of the
environmental variables $\mathbf{X}=(X_{1},\ldots,X_{p})$, corresponding to
the storm peak sea state, are described using the conditional extremes model
of Heffernan and Tawn (2004). For our illustrative example, in Section 2,
$p=2$ but we present the methodology here and for environmental contours for
dimension $p$, to cover more general cases. This asymptotically justified
flexible framework allows for the characterisation of joint tail behaviour
from a sample of independently identically distributed observations of
$\mathbf{X}$, without the need for making a subjective choice for a particular
form of extremal dependence model (copula) between variables. This method has
been applied extensively to the modelling of oceanographic data (e.g.,
Jonathan et al. 2014; Towe et al. 2019; Shooter et al. 2021; Tendijck et al.
2023).

The conditional extremes method uses univariate extreme value techniques to
characterise the distribution of each variable individually, with the joint
structure specified for variables on standard (typically Laplace) marginal
scales (e.g., Keef et al. 2013). Estimation of the conditional extremes model
is thus performed in two stages: (a) marginal extreme value modelling of each
variable $X_{j}$ ($j=1,\ldots,p$) in turn, followed by the marginal
transformation $X_{j}\mapsto Y_{j}$ of each variable to standard Laplace scale
$Y_{j}$ and (b) estimation of the conditional extremes model for the set of
Laplace-scale variables $\mathbf{Y}=(Y_{1},\ldots,Y_{p})$. Subsequently, we
estimate $f_{\mathbf{X}}$ within an extreme joint tail region using the fitted
conditional extremes model. The above steps are discussed below.

####  3.2.2 Marginal modelling and marginal transformation to Laplace scale

We adopt the approach of Davison and Smith (1990) for marginal modelling of
storm peak variable $X_{j}$ for $j=1,\ldots,p$. We fit a generalised Pareto
distribution (GPD) to exceedances of high threshold $u_{j}$, and model
threshold non-exceedances empirically. Our marginal model $F_{X_{j}}$ for the
cumulative distribution function of $X_{j}$ can thus be written

$F_{X_{j}}(x)=\begin{cases}\tilde{F}_{X_{j}}(x)&x\leq u_{j}\\\
\tilde{F}_{X_{j}}(u_{j})+\\{1-\tilde{F}_{X_{j}}(u_{j})\\}F_{\text{GPD},j}(x;u_{j},\sigma_{j},\xi_{j})&x>u_{j},\end{cases}$
(10)

where $\tilde{F}_{X_{j}}$ is the empirical distribution of $X_{j}$ and

$F_{\text{GPD},j}(x;u_{j},\sigma_{j},\xi_{j})=1-\left(1+\frac{\xi_{j}(x-u_{j})}{\sigma_{j}}\right)_{+}^{-1/\xi_{j}},$

for $x>u_{j}$, with scale and shape parameters $\sigma_{j}>0$ and
$\xi_{j}\in\mathbb{R}$ and where $y_{+}=\max(y,0)$ for $y\in\mathbb{R}$. The
values of $u_{j}$ ($j=1,\ldots,p$) are selected using the univariate extreme
threshold selection methods summarised by Coles (2001), see SM. Parameters
$\sigma_{j}$ and $\xi_{j}$ are jointly estimated using standard maximum
likelihood techniques. The probability integral transform

$Y_{j}=\begin{cases}\log\left\\{2F_{X_{j}}\left(X_{j}\right)\right\\}&\text{
for }X_{j}<F_{X_{j}}^{-1}(0.5)\\\
-\log\left\\{2\left[1-F_{X_{j}}\left(X_{j}\right)\right]\right\\}&\text{ for
}X_{j}>F_{X_{j}}^{-1}(0.5),\end{cases}$ (11)

is then applied to each variable in turn, obtaining the standard Laplace-
scaled equivalents $\mathbf{Y}$ of $\mathbf{X}$, such that, for
$j=1,\ldots,p$, the distribution of $Y_{j}$ is

$F_{Y_{j}}(y)=\begin{cases}\frac{1}{2}\exp\left({y}\right)&y\leq 0\\\
1-\frac{1}{2}\exp\left(-{y}\right)&y>0.\end{cases}$

####  3.2.3 Joint dependence modelling

Having transformed environmental variables $\mathbf{X}$ to standard Laplace-
scale equivalents $\mathbf{Y}$, we now apply the model of Heffernan and Tawn
(2004) to estimate the joint distribution $F_{\mathbf{Y}}$ in the upper tail
region. It is shown by Heffernan and Tawn (2004) and Keef et al. (2013) that,
for $j=1,\ldots,p$, there exist unique values of parameters
$\boldsymbol{\alpha}_{|j}\in[-1,1]^{p-1}$,
$\boldsymbol{\beta}_{|j}\in(-\infty,1]^{p-1}$, satisfying the constraints of
Keef et al. (2013), and $\mathbf{z}_{|j}\in\mathbb{R}^{p-1}$, $y>0$, such that

$\lim_{v_{j}\rightarrow\infty}\Pr\left(\frac{\mathbf{Y}_{-j}-\boldsymbol{\alpha}_{|j}Y_{j}}{Y_{j}^{\boldsymbol{\beta}_{|j}}}<\mathbf{z}_{|j},Y_{j}-v_{j}>y|Y_{j}>v_{j}\right)=e^{-y}G_{|j}(\mathbf{z}_{|j}),$
(12)

where the $(p{-}1)$-vector $\mathbf{Y}_{-j}$ denotes the $p$-vector
$\mathbf{Y}$ with $j$th element $Y_{j}$ removed, $G_{|j}$ is a
$(p{-}1)$-dimensional distribution function with non-degenerate marginals, and
componentwise operations are assumed. Property (12) can be leveraged by
assuming a non-linear regression of $\mathbf{Y}_{-j}$ onto $Y_{j}$ holds for
all values of $\mathbf{Y}$ within the region
$\\{\mathbf{Y}\in\mathbb{R}^{p}:Y_{j}>v_{j}\\}$, for some suitably large
finite threshold $v_{j}>0$. For conditioning variable $Y_{j}$
($j=1,\ldots,p$), the form of this regression is

$(\mathbf{Y}_{-j}|\\{Y_{j}=y\\})=\boldsymbol{\alpha}_{|j}y+y^{\boldsymbol{\beta}_{|j}}\mathbf{Z}_{|j},\quad
y>v_{j},$ (13)

where $\mathbf{Z}_{|j}\sim G_{|j}$ is a $(p{-}1)$-dimensional residual random
variable that is independent of $Y_{j}$ given $Y_{j}>v_{j}$. We estimate
parameter vectors $\boldsymbol{\alpha}_{|j}$ and $\boldsymbol{\beta}_{|j}$
using standard maximum likelihood techniques, assuming for model fitting only
that $G_{|j}$ corresponds to independent Gaussian distributions with unknown
means and variances. The distribution $G_{|j}$ is modelled via the kernel
density estimate of the observed values of the $(p{-}1)$-dimensional residual

$\mathbf{Z}_{|j}=\frac{\mathbf{Y}_{-j}-\boldsymbol{\alpha}_{|j}Y_{j}}{Y_{j}^{\boldsymbol{\beta}_{|j}}},\quad\text{for
}\,Y_{j}>v_{j},$ (14)

as in Winter and Tawn (2017).

####  3.2.4 Simulation under the conditional model and estimation of the
environment joint density

Inferences in $\mathbb{R}^{p}$ using the fitted conditional extremes model are
typically made by careful combination of Laplace-scale simulations in each of
the upper tail regions $\\{\mathbf{Y}\in\mathbb{R}^{p}:Y_{j}>v_{j}\\}$, for
$j=1,\ldots,p$, together with empirical estimation in the remaining region
$\\{\mathbf{Y}\in\mathbb{R}^{p}:Y_{j}\leq v_{j}\,\forall j\\}$, as described
in Heffernan and Tawn (2004), to give a set of size $N_{\text{sim}}$
realisations from the estimate of the joint density
${f}_{\mathbf{\mathbf{Y}}}$. The $p$ fitted marginal models (10) can then be
used componentwise to transform this sample of $\mathbf{Y}$ with Laplace-scale
marginals to a sample of $\mathbf{X}$ on the original physical scale.

We can use the simulated sample to estimate the probability of $\mathbf{X}$
being in sub-regions of $\mathbb{R}^{p}$. Specifically, if $D$ is the set of
feasible $\mathbf{X}$ values such that
$\Pr\left(\mathbf{X}\in\mathbb{R}^{p}\setminus D\right)=0$, then we consider a
partition $(D_{1},\ldots,D_{M})$ of $D$. Then, if $N_{\text{sim},i}$ is the
number of realisations in set $D_{i}$, we can estimate $\Pr\left(\mathbf{X}\in
D_{i}\right)$, for any $i=1,\ldots,M$, as ${\Pr}({\mathbf{X}\in
D_{i}})={N_{\text{sim},i}}/N_{\text{sim}}$. To obtain an estimate
$f_{\mathbf{X}}(\mathbf{x})$ for any $\mathbf{x}\in D$ we exploit the property
that, if $\mathbf{x}\in D_{i}$ and $|D_{i}|$ is sufficiently small that
${f}_{\mathbf{X}}$ is reasonably constant for all $\mathbf{x}\in D_{i}$, then

${\Pr}({\mathbf{X}\in D_{i}})=\int_{\mathbf{x}^{\prime}\in
D_{i}}{f}_{\mathbf{X}}(\mathbf{x}^{\prime})\text{d}\mathbf{x}^{\prime}\approx|D_{i}|{f}_{\mathbf{X}}(\mathbf{x}),$
(15)

yielding the estimate $\hat{f}_{\mathbf{X}}(\mathbf{x})=N_{\text{sim,
i}}/(N_{\text{sim}}|D_{i}|)$ for all $\mathbf{x}\in D_{i}$. We can achieve the
required conditions for the approximation to be reliable by taking $M$ to be
sufficiently large and selecting all $D_{i}$ such that $|D_{i}|\propto
M^{-1}$.

###  3.3 Estimation of maximum response in a storm peak sea state given storm
peak variables

####  3.3.1 Outline of estimation of $F_{R_{L}|\mathbf{X}}$

To derive properties of $R_{L}$ we first need to model the behaviour of the
maximum response $R_{I}$ to an individual wave in sea state $\mathbf{X}$.
However, to estimate the distribution $F_{R_{I}|\mathbf{X}}$ of the maximum
response due to the action of an individual wave in sea state $\mathbf{X}$, we
first evaluate the distribution $F_{R_{I}|(\mathbf{X},C)}$ of the maximum
response $R_{I}$ to an individual wave in the sea state $\mathbf{X}$ with
crest elevation $C$. This is achieved by simulation of wave fields under sea
state conditions $\mathbf{X}$ with known crest elevation $C$, followed by
propagation of the resulting stochastic wave fields through to the structural
response model; see Section 3.3.2 for details. Then, integrating out $C$, we
have

$F_{R_{I}|\mathbf{X}}(r|\mathbf{x})=\int_{\mathbb{R}^{+}}F_{R_{I}|\mathbf{X},C}(r|\mathbf{x},c)f_{C|\mathbf{X}}(c|\mathbf{x})\text{d}c,$
(16)

for $r>0$, where $f_{C|\mathbf{X}}$ is the density of crest elevation in the
sea state $\mathbf{X}$, where we assume that crests are Rayleigh-distributed,
with density

$f_{C|\mathbf{X}}(c|\mathbf{x})=\frac{16c}{h(\mathbf{x})^{2}}\exp\left(-8\frac{c^{2}}{h(\mathbf{x})^{2}}\right).$
(17)

for $c>0$, with sea state $\mathbf{X}=\mathbf{x}$ with significant wave height
$h(\mathbf{x})=x_{1}$. Computationally efficient estimation of
$F_{R_{I}|\mathbf{X}}$ following (16) is achieved using importance sampling;
see Section 3.3.3. Finally, we obtain $F_{R_{L}|\mathbf{X}}$ from
$F_{R_{I}|\mathbf{X}}$ by assuming that (a) there are a fixed number of waves
per $L$-hour sea state $\mathbf{x}$, given by
$Q_{L}(\mathbf{x})=60^{2}L/t_{2}(\mathbf{x})$ where $t_{2}(\mathbf{x})$ is the
second spectral moment wave period for the sea state, and (b) individual-wave
maximum responses (i.e., the $R_{I}$) in a given sea state are independent of
each other. Assumption (a) approximates the stochastic number of waves per sea
state with an ‘average’ value, and (b) holds since individual base shears
calculated in Section 3.3.2 are observed for fractions of a second,
significantly less than the typical wave period; therefore, there is no
correlation between responses induced by different waves for a known sea
state. Combining these assumptions gives

$F_{R_{L}|\mathbf{X}}(r|\mathbf{x})=\\{F_{R_{I}|\mathbf{X}}(r|\mathbf{x})\\}^{Q_{L}(\mathbf{x})},$
(18)

for $r>0$, following the definition for the distribution of the maximum of
independent random variables.

####  3.3.2 Simulation of maximum response to the action of an individual
wave, given sea state variables and crest elevation

We estimate the distribution of $R_{I}|(\mathbf{X},C)$ in two stages: (a)
simulation of realisations of wave fields under sea state conditions
$\mathbf{X}$ with known crest elevation $C$, followed by (b) propagation of
the resulting wave fields through a suitable structural response model. The
details of each stage are outlined below.

The model of Taylor et al. (1997) allows for conditional simulation of a wave
field given the occurrence of a turning point of surface elevation in time,
with specified crest elevation $C=c>0$ at the structure location at time
$t=0$, for a given sea state $\mathbf{X}=\mathbf{x}$ with wave spectrum
$S(\cdot;\mathbf{x})$, using linear wave theory. The JONSWAP spectrum of
Hasselmann et al. (1973) is chosen as the form of $S(\cdot;\mathbf{x})$ due to
its applicability to the North Sea wave conditions (Holthuijsen, 2010), see SM
for further details. Taylor et al. (1997) provides expressions for linear
crest elevation $E(t;\mathbf{x},c)$, horizontal velocity
$U(t;z,\mathbf{x},c)$, and horizontal acceleration
$\dot{U}(t;z,\mathbf{x},c)$, at time $t\in\mathbb{R}$ and vertical position
$z\in\mathbb{R}$, relative to the mean water level, each conditioned on the
wave process (a) attaining a turning point of $E$ at time $t=0$, with (b)
$E(t=0;\mathbf{x},c)=c$, both at the location of the structure. The forms of
$E(t;\mathbf{x},c)$ and $U(t;z,\mathbf{x},c)$ are

$\displaystyle E(t;\mathbf{x},c)$
$\displaystyle=\sum_{n=1}^{N}\left\\{(A_{n}+\mathcal{Q}\sigma^{2}_{n})\cos\left(\omega_{n}t\right)+(B_{n}+\mathcal{R}\sigma^{2}_{n}\omega_{n})\sin\left(\omega_{n}t\right)\right\\},$
(19)

and

$U(t;z,\mathbf{x},c)=\sum_{n=1}^{N}\omega_{n}\frac{\cosh\left(k_{n}(d+z)\right)}{\sinh(k_{n}d)}\left\\{(A_{n}+\mathcal{Q}\sigma^{2}_{n})\cos(\omega_{n}t)+(B_{n}+\mathcal{R}\sigma^{2}_{n}\omega_{n})\sin(\omega_{n}t)\right\\},$
(20)

for $z<E(t;\mathbf{x},c)$ and zero otherwise, for a regular grid of angular
frequencies $\omega_{1},\ldots,\omega_{N}>0$ with spacing $\delta\omega>0$ and
$N\in\mathbb{N}$ specified below, where $A_{n},B_{n}\in\mathbb{R}$
($n=1,\ldots,N)$ and $\mathcal{Q},\mathcal{R}\in\mathbb{R}$ are random
coefficients, $d$ is water depth and $k_{n}$ ($n=1,\ldots,N$) are wave numbers
given implicitly by $\omega_{n}^{2}=gk_{n}\tanh(k_{n}d)$. The equation for
$\dot{U}(t;z,\mathbf{x},c)$ can be found by differentiation of
$U(t;z,\mathbf{x},c)$ with respect to $t$. Coefficients $A_{n},B_{n}$
$(n=1,\ldots,N)$ are a series of independently and identically distributed
$N(0,\sigma^{2}_{n})$ random variables with variance
$\sigma^{2}_{n}=S(\omega_{n};\mathbf{x})\delta\omega$, the integrated spectral
density in the frequency band
$(\omega_{n}-\delta\omega/2,\omega_{n}+\delta\omega/2)$ of the discretised
wave spectrum. The random coefficients $\mathcal{Q}$ and $\mathcal{R}$ are
defined as

$\mathcal{Q}=\frac{1}{\sum_{n}\sigma^{2}_{n}}\left(c-\sum_{n=1}^{N}A_{n}\right)\,\,\text{and}\,\,\mathcal{R}=\frac{1}{\sum_{n}\omega_{n}^{2}\sigma^{2}_{n}}\left(-\sum^{N}_{n=1}\omega_{n}B_{n}\right).$

Next, we estimate the total base shear response of the structure to the
simulated conditional wave field. We assume the wave-structure interaction to
be quantified by the equation of Morison et al. (1950), which estimates drag
and inertial loads applied by the ocean environment on a stick structure.
These loads are calculated from the wave velocity and acceleration fields
respectively. Under the assumptions of linear wave theory, these fields can be
derived entirely from knowledge of the wave spectrum. For the applications
described in this work, we assume that the values of sea state $(H_{S},S_{2})$
are sufficient to define the wave spectrum; hence, for a given structure, the
2-dimensional storm peak representation
$\mathbf{X}=(X_{1},X_{2})=(H_{S},S_{2})$ from Section 2 is sufficient to
describe the extreme ocean environment and its associated structural response.

Under our simplifying assumptions, waves are assumed to be unidirectional,
propagating in a single direction towards the vertical cylindrical structure
with nominal small diameter. Waves are assumed to pass through the structure,
whilst also exerting force, without being obstructed, and the effects of
current and wind are ignored. This wave field model provides a basis to
approximate the induced load on a jacket structure. The Morison loading
equation estimates the base shear $M(t;z,\mathbf{x},c)$ induced on a cylinder
by the wave at time $t$ and vertical position $z$, and is given by

$M(t;z,\mathbf{x},c)=\rho c_{m}(z)V\dot{U}(t;z,\mathbf{x},c)+\frac{1}{2}\rho
c_{d}(z)AU(t;z,\mathbf{x},c)|U(t;z,\mathbf{x},c)|,$ (21)

where $c_{m}(z),c_{d}(z)>0$ are inertia and drag coefficients, $\rho=1024$
(recall that SI units are used throughout) is the density of water, $V$ is the
volume of the body and $A$ is the area of the structure perpendicular to the
wave propagation. We assume a cylindrical structure with diameter 1 and height
of 150 situated within water of depth $d=100$. Since the probability of a
crest elevation greater than 50 is near zero for all relevant sea states, this
structure scenario amounts to a cylinder of infinite height. In order to
approximate models of different structure types, $c_{m}(z)$ and $c_{d}(z)$ can
be made to vary with $z$, as discussed in Section 5. To evaluate the total
base shear $B_{S}(t;\mathbf{x},c)$ on the structure at time $t$, we integrate
$M(t;z,\mathbf{x},c)$ to give

$B_{S}(t;\mathbf{x},c)=\int_{-d}^{E(t;\mathbf{x},c)}M(t;z,\mathbf{x},c)\text{d}z,$
(22)

the total Morison load induced up the water column at the structure location.

The response $R_{I}|(\mathbf{X}=\mathbf{x},C=c)$ may be obtained by
considering the portion of the time series
$\\{B_{S}(t;\mathbf{x},c)\\}_{t\in\mathbb{R}}$ that corresponds to the central
wave conditioned to attain $E(t=0;\mathbf{x},c)=c$; that is, the period of
time $t\in\mathcal{T}_{0}\subset\mathbb{R}$, with $0\in\mathcal{T}_{0}$, for
which the wave surrounding the conditioning crest of elevation $c$ at time
$t=0$ acts on the structure. We define

$R_{I}|(\mathbf{X}=\mathbf{x},C=c)=\max_{t\in\mathcal{T}_{0}}{B_{S}(t;\mathbf{x},c)}.$
(23)

We obtain realisations of $R_{I}|(\mathbf{X}=\mathbf{x},C=c)$ from a time
series of the base shear response (22) evaluated using Morison loads (21), in
turn calculated from wave fields simulated according to expressions (19) and
(20). The interval of time over which wave fields are simulated corresponds to
a period of $120$ seconds, sufficiently large to ensure reliable performance
of the FFT algorithm (Cooley and Tukey, 1965), meaning here
$\mathcal{T}_{0}\subset[-60,60]$. Realisations of conditional crest elevation
and wave kinematics are simulated for a regular grid
$\\{(t_{i},z_{j})\\}_{i=1,j=1}^{n_{t},n_{z}}$ of values $t\in[-60,60]$ and
$z\in[-100,150]$. We set $N=n_{t}$ in expressions (19) and (20), which is
necessary to evaluate the wave field equations using the FFT algorithm; see SM
for details. The values of $n_{t}=480$ and $n_{z}=50$ are chosen, sufficient
to ensure reasonable response approximation. The simulated kinematics are then
propagated through the Morison equation, providing a realisation of
$\\{M(t_{i};z_{j},\mathbf{x},c)\\}_{i=1,j=1}^{n_{t},n_{z}}$. Numerical
evaluation of integral (22) with respect to $z$ yields a realisation of the
time series $\\{B_{S}(t_{i};\mathbf{x},c)\\}_{i=1}^{n_{t}}$. A realisation of
the maximum individual wave response $R_{I}|(\mathbf{X}=\mathbf{x},C=c)$ is
then obtained by applying (23), using the set
$t\in\\{t_{i}\\}_{i=1}^{n_{t}}\cap\mathcal{T}_{0}$ as an approximation to
$\mathcal{T}_{0}$. Given conditioning crests $\\{c_{i}\\}_{i=1}^{k}$, the
above procedure can be used to map $c_{i}\mapsto r_{i}$, for $i=1,\ldots,k$,
obtaining a set of maximum responses $\\{r_{i}\\}_{i=1}^{k}$, for a given sea
state $\mathbf{x}$.

####  3.3.3 Importance sampling of simulated maximum responses

The procedure discussed in Section 3.3.2 is used to obtain realisations
$\\{r_{i}\\}_{i=1}^{k}$ of ${R_{I}|(\mathbf{X},C)}$, for a set of $k$
conditioning crests $\\{c_{i}\\}_{i=1}^{k}$ and specified values of storm peak
variables $\mathbf{X}$. These are then used in integral (16) to estimate the
distribution $F_{R_{I}|\mathbf{X}}$ of the maximum response to an individual
wave in sea state $\mathbf{X}$, assuming a Rayleigh distribution (17) for
$C|\mathbf{X}$. However, evaluation of integral (16) via Monte Carlo methods
sampling from the Rayleigh density is inefficient in targeting the tail of the
response distribution $F_{R_{I}|\mathbf{X}}$. Given our interest in the
extreme structural response on the structure, we therefore employ the
importance sampling approach described by Towe et al. (2021), writing integral
(16) as

$F_{R_{I}|\mathbf{X}}(r|\mathbf{x})=\int_{r}F_{R_{I}|\mathbf{X},C}(r|\mathbf{x},c)\frac{f_{C|\mathbf{X}}(c|\mathbf{x})}{g^{(\epsilon)}_{C|\mathbf{X}}(c|\mathbf{x})}g^{(\epsilon)}_{C|\mathbf{X}}(c|\mathbf{x})\text{d}c,$
(24)

for $r>0$, where $g^{(\epsilon)}_{C|(\mathbf{X}=\mathbf{x})}$ is the density
of the $\text{Uniform}[0,\epsilon h(\mathbf{x})]$ distribution, for
significant wave height $h(\mathbf{x})$ and some $\epsilon>0$. For a fixed
number of conditional wave simulations, sampling of the conditioning crest $c$
from $g^{(\epsilon)}_{C|\mathbf{X}}$ ensures greater coverage of the feasible
range of large crest elevations and of the induced maximum response than is
achieved when sampling $c$ from $f_{C|\mathbf{X}}$. Therefore, the upper tail
of the distribution $F_{R_{I}|\mathbf{X}}$ is more efficiently estimated using
the sampling distribution $g^{(\epsilon)}_{C|\mathbf{X}}$. The value of
$\epsilon$ in (24) is selected so that $g^{(\epsilon)}_{C|\mathbf{X}}$
provides adequate coverage of the domain of $f_{C|\mathbf{X}}$, i.e., the
exceedance probability

$\Pr\\{C>\epsilon h(\mathbf{x})|\mathbf{X}=\mathbf{x}\\}=\int_{c>\epsilon
h(\mathbf{x})}f_{C|\mathbf{X}}(c|\mathbf{x})\text{d}c=\exp\left(-8\epsilon^{2}\right),$
(25)

is sufficiently close to zero. We set $\epsilon=2$, which gives a value of
probability (25) in the order of $10^{-14}$.

Integral (24) is then estimated as follows. The $k$ conditional crests
$\\{c_{i}\\}_{i=1}^{k}$ are sampled from the uniform proposal density
$g^{(\epsilon)}_{C|\mathbf{X}}$. Corresponding realisations
$\\{r_{i}\\}_{i=1}^{k}$ of single-wave maximum responses are obtained using
the procedure described in Section 3.3.2. The distribution
$F_{R_{I}|\mathbf{X}}$ is then estimated as

$\hat{F}_{R_{I}|\mathbf{X}}(r|\mathbf{x})=\frac{\sum_{i=1}^{k}\mathbbm{1}_{\\{r_{i}\leq
r\\}}\iota^{(\epsilon)}(c_{i}|\mathbf{x})}{\sum^{k}_{i=1}\iota^{(\epsilon)}(c_{i}|\mathbf{x})},$
(26)

for $r>0$, where
$\iota^{(\epsilon)}(c|\mathbf{x})=f_{C|\mathbf{X}}(c|\mathbf{x})/{g^{(\epsilon)}_{C|\mathbf{X}}(c|\mathbf{x})}$
is the importance sampling ratio and $\mathbbm{1}_{\\{r_{i}\leq r\\}}=1$ if
$r_{i}\leq r$, zero otherwise, for $i=1,\ldots,k$. Estimate (26) is an
empirical cumulative distribution function of the simulated responses,
weighted to remove bias introduced from sampling crests from
$g^{(\epsilon)}_{C|\mathbf{X}}$ rather than $f_{C|\mathbf{X}}$. We use
estimate (26) to evaluate the distribution of maximum response per $L$-hour
sea state using relation (18). Given an estimate for $f_{\mathbf{X}}$ obtained
as in Section 3.2, we may then calculate the marginal maximum response
distribution using integral (7).

##  4 Environmental contours

###  4.1 Overview of environmental contours

Environmental contours provide a method of determining extremal conditions
which are in some way related to an extreme structural response. These
contours often make assumptions about the interaction between environment and
response, usually regarding the shape of some failure boundary in the
environment space such that environmental conditions beyond the boundary will
result in structural failure. For instance, IFORM (e.g., Winterstein et al.
1993) contours assume a convex form for this boundary, whereas ISORM (e.g.,
Chai and Leira 2018) assumes it to be concave. These assumptions may or may
not be valid depending on the specific features of the structure type in
question. Here, we outline the methodology of the IFORM contour (Section 4.2)
and the fitting approach we employ to estimate it for our example dataset
(Section 4.3).

###  4.2 IFORM design contours

Section 3.1.4 details how, given an estimate for the distribution of
$R_{L}|\mathbf{X}$ from Section 3, in principle we can evaluate the
probability of structural failure
$\Pr\left(\mathcal{F}_{L}|\mathbf{X}=\mathbf{x}\right)$ for a sea state (of
duration $L$) with variables $\mathbf{X}=\mathbf{x}$. IFORM offers an approach
to structural design which avoids direct calculation of
$\Pr\left(\mathcal{F}_{L}|\mathbf{X}=\mathbf{x}\right)$, by attempting to make
conservative assumptions. For an ocean environment represented by a set of
random variables $\mathbf{X}$ transformed to independent standard Gaussian
random variables $\mathbf{U}$, IFORM assumes that
$\Pr\left(\mathcal{F}_{L}|\mathbf{U}=\mathbf{u}\right)$ is deterministic for
all $\mathbf{u}$, taking values $\\{0,1\\}$, contrary to the failure
probability discussed in Section 3.1.4 which takes any value in $[0,1]$.
Writing the region of environmental space corresponding to failure
$\mathcal{F}_{L}$ as $\mathcal{F}$, IFORM assumes that the boundary
$\partial\mathcal{F}$ of $\mathcal{F}$ is linear, and lies tangential to a
contour of constant transformed environmental density $f_{\mathbf{U}}$, making
the assumed location of $\partial\mathcal{F}$ dependent on the joint Gaussian
distribution $f_{\mathbf{U}}$. The assumption of a failure boundary of this
type is typically conservative, in that estimates for
$\Pr(\mathbf{U}\in\mathcal{F})$ using it have positive bias.

The transformation of $\mathbf{X}\mapsto\mathbf{U}$ is achieved via the method
of Rosenblatt (1952), which proceeds as follows. For storm peak variables
$\mathbf{X}=(X_{1},\ldots,X_{p})$, suppose we can estimate the nested
conditional distributions
$F_{X_{1}},F_{X_{2}|X_{1}},F_{X_{3}|(X_{1},X_{2})},\ldots,F_{X_{p}|(X_{1},\ldots,X_{p-1})}$.
Estimation of these distributions is non-trivial as it involves estimating a
sequence of conditional dependence models with results dependent on the
sequencing of the $p$ environmental variables; see Section 4.3 for an example
approach. The Rosenblatt transformation maps a realisation
$\mathbf{x}=(x_{1},\ldots,x_{p})$ of $\mathbf{X}$ to the realisation
$\mathbf{u}=(u_{1},\ldots,u_{p})$ of $\mathbf{U}$ via
$u_{1}=\Phi^{-1}\\{F_{X_{1}}(x_{1})\\}$ and
$u_{j}=\Phi^{-1}\left\\{F_{X_{j}|(X_{1},\ldots,X_{j-1})}(x_{j}|x_{1},\ldots,x_{j-1})\right\\}$,
for $j=2,\ldots,p$, where $\Phi$ is the standard Gaussian cumulative
distribution function.

In $\mathbf{U}$-space, contours of constant probability joint density
correspond to boundaries of hyperspheres centred at the origin. In particular,
the $P$-year IFORM contour in $\mathbf{U}$-space is the boundary of a
hypersphere with radius

$\beta_{I}=\Phi^{-1}\left(1-\frac{1}{N_{\text{an}}P}\right),$ (27)

where $N_{\text{an}}$ is the average number of independent storm peak
observations per annum. That is, the probability of a point lying outside the
set enclosed by the $P$-year IFORM contour is $1/(N_{\text{an}}P)$. Consider
failure region $\mathcal{F}_{\beta_{I}}$ with boundary
$\delta\mathcal{F}_{\beta_{I}}$ tangential to the hypersphere centred at the
origin with radius $\beta_{I}$. For any angle $\boldsymbol{\psi}$, in
spherical polar coordinates for a point on the hypersphere, we have

$\Pr\left(\mathbf{U}\in\mathcal{F}_{\beta_{I}}\right)=\Pr\left(\sum_{i=1}^{p}w_{i}(\boldsymbol{\psi})U_{i}>\beta_{I}\right),$
(28)

where $\sum^{p}_{i=1}w_{i}^{2}(\boldsymbol{\psi})=1$, due to the linearity of
$\delta\mathcal{F}_{\beta_{I}}$ and it being tangential to the hypersphere at
radius $\beta_{I}$. So, as $W=\sum_{i=1}^{p}w_{i}(\boldsymbol{\psi)}U_{i}\sim
N(0,1)$ given the independence of $(U_{1},\ldots,U_{p})$, it follows that
$\Pr\left(\mathbf{U}\in\mathcal{F}_{\beta_{I}}\right)=\Pr\left(W>\beta_{I}\right)$,
which directly gives expression (27) for $\beta_{I}$.

Once the environmental contour has been estimated in $\mathbf{U}$-space, it
can be represented in the original $\mathbf{X}$-space via the inverse
Rosenblatt transformation $x_{1}=F^{-1}_{X_{1}}\\{\Phi(u_{1})\\}$ and
$x_{j}=F^{-1}_{X_{j}|(X_{1},\ldots,X_{j-1})}\\{\Phi(u_{j})|x_{1},\ldots,x_{j-1}\\}$,
for $j=2,\ldots,p$. Unlike the hypersphere-shaped contours in the
$\mathbf{U}$-space, these contours are not guaranteed to be convex (see
Section 5.2). The procedure for construction of a $P$-year IFORM contour in
terms of environmental variable $\mathbf{X}=(X_{1},\ldots,X_{p})$ is
summarised in Algorithm 1.

Algorithm 1 IFORM contour calculation for $\mathbf{X}=(X_{1},\ldots,X_{p})$

Input Return period $P$; Average number of independent storm peaks per annum
$N_{\text{an}}$; Estimates of distributions
$F_{X_{1}},F_{X_{2}|X_{1}},F_{X_{3}|(X_{1},X_{2})},\ldots,F_{X_{p}|(X_{1},\ldots,X_{p-1})}$.
Output $P$-year IFORM contour.

1:Define $\beta_{I}=\Phi^{-1}\left(1-\frac{1}{N_{\text{an}}P}\right)$ where
$\Phi$ is the standard normal cdf.

2:Obtain a set of $k$ equally spaced points
$\\{(u_{1(j)},\ldots,u_{p(j)})\\}^{k}_{j=1}$ on the hypersphere given by
$u_{1}^{2}+\ldots+u_{p}^{2}=\beta_{I}^{2}$, using a regular grid of values of
$\boldsymbol{\psi}$ in spherical polar coordinates.

3:for $j$ in $1,\ldots,k$ do

4: Compute $x_{1(j)}=F_{X_{1}}^{-1}\\{\Phi(u_{1(j)})\\}$.

5: Compute $x_{2(j)}=F_{X_{2}|X_{1}}^{-1}\\{\Phi(u_{2(j)})|x_{1(j)}\\}$.

6: ⋮

7: Compute
$x_{p(j)}=F^{-1}_{X_{p}|(X_{1},\ldots,X_{p-1})}\\{\Phi(u_{p(j)})|(x_{1(j)},\ldots,x_{p-1(j)})\\}$.
return $\\{(x_{1(j)},\ldots,x_{p(j)})\\}_{j=1}^{k}$ points along the
$p$-dimensional IFORM contour.

###  4.3 Joint parametric models for storm peak variables

Construction of the IFORM contour for environmental variables
$\mathbf{X}=(X_{1},X_{2})=(H_{S},S_{2})$ via Algorithm 1 requires estimates
for the marginal distribution $F_{H_{S}}$ and conditional distribution
$F_{S_{2}|H_{S}}$ (while $F_{S_{2}}$ and $F_{H_{S}|S_{2}}$ could equally be
used, we use the former to reflect past approaches). Mirroring the
hierarchical approach of Winterstein et al. (1993), we select the GPD tail
model (10) for marginal $H_{S}$, and evaluate a range of parametric forms for
the distribution of $S_{2}|H_{S}$, selecting the most appropriate based on an
assessment of predictive performance. We estimate the model $S_{2}|H_{S}$ as
follows.

We allow for two sources of flexibility: (a) the conditional distributional
form for $S_{2}|H_{S}$, and (b) the nature of the parametric form for how the
parameters of the distribution vary as a function of $H_{S}$. The
distributions we consider in modelling step (a) are the Lognormal($\mu_{L}$,
$\sigma_{L}$), as in Winterstein et al. (1993), Gamma($\alpha$, $\beta$),
Weibull($\lambda$, $k$) and the Generalised Extreme Value distribution
GEV($\mu_{G}$, $\sigma_{G}$, $\xi$); see SM. We also consider conditional
distributions fitted to the transformed negative steepness
$S^{-}_{2}=\min(S_{2})-S_{2}$, enabling the right-hand tail of the
distribution to be fitted to small $S_{2}$; see SM for a summary of all model
combinations. In step (b), we impose linear, quadratic and exponential forms
for the $S_{2}$ distribution parameters as functions of $H_{S}$. To assess the
performance of candidate models for $S_{2}|H_{S}$ and $S_{2}^{-}|H_{S}$, we
use a cross validation approach in which we evaluate the predictive likelihood
of the models, focusing on the performance in the tail region (for large
$H_{S}$). This proceeds as follows.

The sample $\\{\mathbf{x}_{i}\\}_{i=1}^{N_{\text{all}}}$ (see Section 2) is
partitioned into a ‘body’ where $H_{S}\leq v$ and a ‘tail’ where $H_{S}>v$ for
$v>0$, denoted $\\{\mathbf{x}_{i}^{B}\\}_{i=1}^{N_{B}}$ and
$\\{\mathbf{x}_{i}^{T}\\}_{i=1}^{N_{T}}$ respectively, with
$N_{B}+N_{T}=N_{\text{all}}$. The tail portion
$\\{\mathbf{x}_{i}^{T}\\}_{i=1}^{N_{T}}$ is itself partitioned into $K$
subsets $\\{\mathcal{S}_{j}\\}_{j=1}^{K}$ each of sizes $\lfloor
N_{T}/K\rfloor$ or $\lceil N_{T}/K\rceil$. A $K$-fold cross validation is
performed using training set
$\\{\mathbf{x}_{i}\\}_{i=1}^{N}\setminus\mathcal{S}_{j}$ and test set
$\mathcal{S}_{j}$ at fold $j\in\\{1,\ldots,K\\}$. That is, we always utilise
the entirety of the body of the data within the training set alongside all but
a single fold of the tail. Excluding a subset of the tail from the training
data in this way is appropriate for estimation of the conditional distribution
for $S_{2}|H_{S}$, but leads to biased estimation of the marginal distribution
of $H_{S}$, so it is important to note estimation of the marginal model (10)
(see Section 3.2.2) is carried out using the entire dataset. The predictive
likelihood is then only calculated on extreme data points, and so measures the
fit of each model to the extremes of the data.

We repeat the above process for $K=5,10$ over values of $v=0,0.8,0.9$ to
determine the sensitivity of model performance to the choice of extreme
threshold. Setting $v=0$ recovers a standard cross validation approach for
assessing fit to all of data, which we also include to ensure the best
performing models fit the body of the data well. In addition, we also evaluate
the model fit using AIC. The AIC and cross validation scores are each
standardised by dividing by the number of observations for which we evaluate
the optimised likelihood (negated when considering AIC, which is standardised
over both terms). This standardisation results in a set of loosely comparable
scores for each model across different threshold choices, and these scores for
each model are averaged over values of both $K$ and $v$ to obtain a single
score, referred to as the aggregate score (AS); see SM for details. The model
with the largest AS is deemed to be the best predictive model for the data.

##  5 Results

###  5.1 Estimating the joint density of storm peaks

We employ the forward methodology of Sections 3.2 and 3.3 to estimate the
environmental density $f_{\mathbf{X}}$ and response distribution
$F_{R_{L}|\mathbf{X}}$ for our motivating dataset of
$\mathbf{X}=(H_{S},S_{2})$ introduced in Section 2. In turn, these are used to
evaluate the distributions $F_{R_{S}}$ and $F_{R_{A}}$ (as in Sections 3.1.2
and 3.1.3) and subsequently $\text{CDE}_{P}$ (as in Section 3.1.5). The
$\text{CDE}_{P}$ is then compared with various IFORM contour estimates using
the methods of Section 4.

We model the joint environment $\mathbf{X}$ using the method of Heffernan and
Tawn (2004) discussed in Section 3.2, fitted to data in the region
$A_{\nu}=\\{\mathbf{X}\in\mathbb{R}^{2}:H_{S}>\nu\\}$ for conditioning
threshold $\nu$. We take $\nu=\tilde{F}_{H_{S}}^{-1}(0.95)$. Inspection of
plots for the variability of conditional extremes model parameters with
respect to threshold indicated this choice of threshold to be within the
interval for which these parameters are invariant.

Using this fitted model, the density $f_{\mathbf{X}}$ is estimated in
$A_{\nu}$ as described in Section 3.2.4. The density in the complement
$\mathbb{R}^{2}\setminus A_{\nu}$ is modelled empirically. Since our interest
lies in environments with large $H_{S}$ and associated structural responses,
we are not concerned with (a) smooth estimation of this lower portion of the
density and (b) the density in the region corresponding to large $S_{2}$ but
small $H_{S}$.

Figure 2: Storm peak density for $H_{S}$ [m] and $S_{2}$ in $A_{\nu}$ on log
scale, estimated from the fitted conditional extremes model. Observations of
storm peak $\mathbf{X}=(H_{S},S_{2})$ are shown as black dots.

In the 2-dimensional case $\mathbf{X}=(H_{S},S_{2})$, the parameters of model
(13) reduce to $\alpha\in[-1,1]$ and $\beta\in(-\infty,1]$ when conditioning
on $H_{S}$. We find corresponding estimates $\hat{\alpha}=0.378$ and
$\hat{\beta}=0.533$, as well as an estimate $\hat{\xi}=-0.02$ for the $H_{S}$
marginal shape parameter $\xi$ in (10), the latter indicating a near
exponential upper tail for $H_{S}$. The effect of $\hat{\alpha}$ is seen in
Figure 2, as a positive trend in the resulting density estimate. This density
also appears to agree with the shape of the data in the extreme region
$\mathcal{A}_{\nu}$.

###  5.2 Selection of model form for the conditional distribution

For IFORM, we consider four distributional forms for each of $S_{2}|H_{S}$ and
$S^{-}_{2}|H_{S}$ (summarised in Section 4.3) each with two distribution
parameters modelled as functions of $H_{S}$. We do not model the GEV shape
$\xi$ as a function of $H_{S}$ as its estimation using maximum likelihood is
difficult for finite samples; instead, we assume it to be an unknown constant.
The variation of each of these eight distribution parameters with $H_{S}$ is
represented by one of three parametric forms (linear, quadratic and
exponential), giving a total of 72 combined candidate models for $S_{2}|H_{S}$
and $S^{-}_{2}|H_{S}$. These models are ranked using the AS, introduced in
Section 4.3, yielding results given in full in the SM. Table 1 summarises
these results, showing the optimal model for each of $S_{2}|H_{s}$ and
$S_{2}^{-}|H_{S}$, with unique $a,b,c\in\mathbb{R}$ for each distributional
parameter, constrained such that their respective domains are not violated.
Standard errors are found as the sample standard deviation of the AS,
evaluated over thirty replicates of each cross validation setting.

Label | Distribution (with optimal functional form of parameters) | AS
---|---|---
$\mathcal{C}_{P}^{1}$ | $S^{-}_{2}|\\{H_{S}=h\\}\sim\text{GEV}(\mu_{G},\sigma_{G},\xi)$ | $\mu_{G}=a+b\exp(ch)$ | $\sigma_{G}=a+b\exp(ch)$ | 3.999 (0.002)
$\mathcal{C}_{P}^{2}$ | $S_{2}|\\{H_{S}=h\\}\sim\text{Weibull}(k,\lambda)$ | $k=a+bc$ | $\lambda=a(h+b)^{2}+c$ | 3.983 (0.003)
$\mathcal{C}_{P}^{3}$ | $S^{-}_{2}|\\{H_{S}=h\\}\sim\text{Lognormal}(\mu_{L},\sigma_{L})$ | $\mu_{L}=a+bh$ | $\sigma_{L}=a+bh$ | 3.963 (0.001)
$\mathcal{C}_{P}^{4}$ | $S^{-}_{2}|\\{H_{S}=h\\}\sim\text{Gamma}(\alpha,\beta)$ | $\alpha=a(h+b)^{2}+c$ | $\beta=a+bh$ | 3.963 (0.006)
$\mathcal{C}_{P}^{5}$ | $S^{-}_{2}|\\{H_{S}=h\\}\sim\text{Weibull}(k,\lambda)$ | $k=a+b\exp(ch)$ | $\lambda=a+b\exp(ch)$ | 3.897 (0.002)
$\mathcal{C}_{P}^{6}$ | $S_{2}|\\{H_{S}=h\\}\sim\text{Gamma}(\alpha,\beta)$ | $\alpha=a+bh$ | $\beta=a+bh$ | 3.873 (0.002)
$\mathcal{C}_{P}^{7}$ | $S_{2}|\\{H_{S}=h\\}\sim\text{Lognormal}(\mu_{L},\sigma_{L})$ | $\mu_{L}=a(h+b)^{2}+c$ | $\sigma_{L}=a+b\exp(ch)$ | 3.824 (0.005)
$\mathcal{C}_{P}^{8}$ | $S_{2}|\\{H_{S}=h\\}\sim\text{GEV}(\mu_{G},\sigma_{G},\xi)$ | $\mu_{G}=a+b\exp(ch)$ | $\sigma_{G}=a(h+b)^{2}+c$ | 3.533 (0.000)

Table 1: AS (with standard error) for the optimal forms of each distribution
considered for $S_{2}|H_{S}$ and $S^{-}_{2}|H_{S}$. Large values of AS
indicate good performance. The three best performing models’ scores are
indicated in bold. Complete results are given in the SM.

The models in Table 1 are used with Algorithm 1 to construct the IFORM
contours in Figure 3. Each contour corresponds to a return period of $P=1000$
years. Contours are labelled $\mathcal{C}_{P}^{1}$ to $\mathcal{C}_{P}^{8}$
and ordered according to their AS, with the best fitting models having the
lowest labelling. All of the contour estimates provide plausible descriptions
of the shape of the sample, but from an engineering design perspective, we
note clear differences in the shape and position of the contours for larger
$H_{S}$. Even so, the three highest ranking models generate contours which
agree to a reasonable degree in all regions. These three contours also appear
visually to be the best descriptions of shape of the data. In comparison, the
other contours do not agree in the region of large $H_{S}$, and fail to
capture the shape of the main body of the data. We therefore select the
highest three ranking contours as the best representations of IFORM to compare
to $\text{CDE}_{P}$ in Section 5.3.

Figure 3: IFORM contours in $H_{S}$ [m] and $S_{2}$ constructed from the
models in Table 1, corresponding to an exceedence probability of
$p=10^{-3}/73$, or a return period of $P=1000$ years for data with
$N_{\text{an}}=73$ observations per annum. Contours $\mathcal{C}_{P}^{1}$ to
$\mathcal{C}_{P}^{8}$ are listed and coloured in order of decreasing
performance AS, from red to purple, with contour $\mathcal{C}^{i}_{P}$
labelled C$i$ for $i=1\ldots,8$.

###  5.3 Estimating the conditional density of associated environmental
variables

####  5.3.1 Estimation of $\text{CDE}_{P}$ for example structure models

We evaluate $\text{CDE}_{P}$ for three examples of the stick structure model
(Section 3.3.2), denoted A, B & C. These structures assume different values
for drag and inertia coefficients $c_{d}(z)$ and $c_{m}(z)$ along their height
$z$, as shown in Table 2. Structure A represents the simplest stick structure,
with homogeneous drag and inertia coefficients along its entire height, i.e.,
$c_{d}(z)=c_{m}(z)=1$ for all $z\in[-100,150]$. We mimic wave-in-deck loads
for structure B, with $c_{d}(z)=c_{m}(z)=100$ increased for a portion of the
structure above sea level. A portion of structure C near the sea bed incurs
increased load.

|  | Value of $c_{d}(z)=c_{m}(z)$
---|---|---
|  | $5<z\leq 15$ | $-95<z\leq-85$ | Elsewhere
Structure | A | 1 | 1 | 1
B | 100 | 1 | 1
| C | 1 | 100 | 1
Table 2: The drag and inertia coefficients of (21) with $c_{d}(z)=c_{m}(z)$
for all $z$, for varying $z$ associated with stick structures A, B & C.

Figure 4 shows the corresponding estimates for $\text{CDE}_{P}$ for $P=1000$
years. The shape and position of the conditional density varies between
structures, due to their differing loading characteristics. For structure C in
particular, the conditional density extends to larger $H_{S}$, and over a
wider interval of $S_{2}$; we comment further on this feature in the
discussion of environmental contours in Figure 6.

Figure 4: Estimated $\text{CDE}_{P}$ for the three structure models given in
Table 2 (left to right A, B & C), for a period of $P=1000$ years. These are
evaluated for the example dataset given in Section 2, using the methods of
Section 3.

####  5.3.2 Comparison between $\text{CDE}_{P}$ and IFORM contours

Provided we have accurate models for the series of $p$ nested conditional
distributions $F_{X_{1}},\ldots,F_{X_{p}|(X_{1},\ldots,X_{p-1})}$, and
provided that the assumptions underlying IFORM are valid, the $P$-year IFORM
contour gives design points at which evaluation of a response model will
provide conservative estimates of the $P$-year response, as indicated by
Winterstein et al. (1993). That is, it aims to provides environmental
conditions at least as extreme as those which induce the $P$-year response for
specific structures. It is natural therefore to consider assessing IFORM
contour performance using $\text{CDE}_{P}$, since the latter provides
asymptotically-justified estimates of the environmental conditions
corresponding to the $P$-year response, obtained from application from the
forward model of Section 3.

Assuming that $\text{CDE}_{P}$ provides a valid estimate of the environment
density conditional on the $P$-year response, we reason that a well-estimated
IFORM contour $\mathcal{C}_{P}$ should intersect with $\text{CDE}_{P}$, and
that the discrepancy between $\mathcal{C}_{P}$ and $\text{CDE}_{P}$ points to
inadequacy of the IFORM methodology, potentially due to a mis-specified
parametric environmental model or invalid IFORM assumptions regarding the
characteristics of the wave-structure interaction for the application at hand.
For instance, see Figure 5, which shows the estimated $\text{CDE}_{P}$ for
structures A, B and C, for $P=1000$. Contour $\mathcal{C}_{P}^{2}$ in the left
panel overlaps with darker points of $\text{CDE}_{P}$ so provides estimates of
the $P$-year response which agree with those obtained from the forward model,
meaning it performs well in this region of the contour. The centre panel of
Figure 5 shows $\mathcal{C}_{P}^{1}$ taking $(H_{S},S_{2})$ values more
extreme than the darker points of $\text{CDE}_{P}$, indicating that the
contour will provide estimates for the $P$-year response more extreme than
those obtained from the forward model, i.e., the conservative outcome
intended. Conversely, the right panel shows $\mathcal{C}_{P}^{2}$ taking
$(H_{S},S_{2})$ values less extreme than the darker points of
$\text{CDE}_{P}$, suggesting estimates for the $P$-year response obtained from
points along this contour will be smaller than those obtained from the forward
model and hence they fail to give conservative estimates, which is contrary to
their claimed properties. More generally, the majority of the presented IFORM
contours exhibit over-conservatism in regions of small $S_{2}$, i.e., they lie
in regions of the environment space with more severe $H_{S}$ than in the
regions with non-zero $\text{CDE}_{P}$, and this could lead to substantial
over-design if this region of the environmental space was important.

We observe from the above analysis of Figure 5 that the overlap between the
region of non-zero $\text{CDE}_{P}$ and the region $\mathcal{A}_{P}$ bounded
below by the IFORM contour $\mathcal{C}_{P}$ relates to the level of
conservatism, and therefore is a good scalar metric for measuring performance
of IFORM contours relative to the more strongly justified forward approach. If
this overlap is zero, then the contour lies in a region less extreme than the
region of non-zero $\text{CDE}_{P}$, indicating non-conservatism. Conversely,
if the overlap includes all of the non-zero region of $\text{CDE}_{P}$, the
contour appears to be conservative. To formalise this observation, we define
the following metric, which quantifies the level of this overlap. Consider

$\zeta(P,\mathcal{A}_{P})=2\int_{\mathbf{x}}\mathbbm{1}_{\mathcal{A}_{P}}(\mathbf{x})f_{\mathbf{X}|R_{L}}(\mathbf{x}|r_{P})\text{d}\mathbf{x}-1,$
(29)

for $f_{\mathbf{X}|R_{L}}$ as in Section 3.1.5, with
$\mathbbm{1}_{\mathcal{A}_{P}}(\mathbf{x})=1$ if
$\mathbf{x}\in\mathcal{A}_{P}$ and
$\mathbbm{1}_{\mathcal{A}_{P}}(\mathbf{x})=0$ otherwise, and where $r_{P}$ is
the $P$-year response of the structure. The metric $\zeta(P,\mathcal{C}_{P})$
takes values in [-1, 1]. Here, $\zeta(P,\mathcal{C}_{P})>0$ indicates
conservatism of $\mathcal{C}_{P}$ due to overestimation of the $P$-year
response, $\zeta(P,\mathcal{C}_{P})<0$ indicates underestimation of the
$P$-year response and so non-conservatism, and
$\zeta(P,\mathcal{C}_{P})\approx 0$ indicates accurate estimates for the
$P$-year response and so optimal overall structural design given by IFORM.

To justify our findings regarding $\zeta(P,\mathcal{C}_{P})$ in a general
case, we apply the following heuristic argument. First, we observe from Figure
4 that $\text{CDE}_{P}$ exhibits approximate marginal symmetry with respect to
both $H_{S}$ and $S_{2}$ (i.e., little marginal skew). We also assume
$\mathcal{A}_{P_{1}}\subset\mathcal{A}_{{P_{2}}}$ for any $P_{2}>P_{1}\geq 1$,
and that a contour $\mathcal{C}_{P}$ varies smoothly with $H_{S}$ and $S_{2}$
for a given $P$. Then, (a) if $\zeta(P,\mathcal{C}_{P})$ $\approx$ 0, the
integral over $\text{CDE}_{P}$ within $\mathcal{A}_{P}$ is approximately equal
to the integral over $\text{CDE}_{P}$ within $\mathcal{A}^{c}_{P}$. We
interpret this as indicating that points on $\mathcal{C}_{P}$ coincide with
high values of $\text{CDE}_{P}$, as in the cases seen in Figure 5, hence
structural responses corresponding to points on $\mathcal{C}_{P}$ will be of
similar magnitudes to those from the high density regions of $\text{CDE}_{P}$.
However, (b) if $\zeta(P,\mathcal{C}_{P})>0$, $\mathcal{A}_{P}$ contains the
environmental region where $\text{CDE}_{P}$ is non-zero, hence points on
$\mathcal{C}_{P}$ produce structural responses beyond the $P$-year level,
resulting in conservative design using IFORM contour $\mathcal{C}_{P}$.
Further, (c) if $\zeta(P,\mathcal{C}_{P})<0$, the intersection between
$\mathcal{A}_{P}$ and the non-zero $\text{CDE}_{P}$ region of the environment
space is negligible. Given our assumptions, this arises only when
$\mathcal{A}_{P}$ occupies a region of the environmental space less extreme
that that with non-zero $\text{CDE}_{P}$. Under these circumstances, points on
$\mathcal{C}_{P}$ will produce structural responses corresponding to return
periods less that $P$, resulting in a lack of conservatism in design.

Estimates for $\zeta(P,\mathcal{C}_{P})$ in Table 3 support this
interpretation, relative to our estimates in Figure 5. For structure A, only
contour $\mathcal{C}_{P}^{2}$ appears to pass through the highest density
region of $\text{CDE}_{P}$. The other contours lie beyond the highest density
region of $\text{CDE}_{P}$, corresponding to $\zeta(P,\mathcal{C}_{P})>0$.
Observations for structure B are similar, since $\text{CDE}_{P}$ does not vary
considerably between structures A and B. For structure C, relative to A and B,
the highest density regions in $\text{CDE}_{P}$ occur at higher $H_{S}$ but
lower $S_{2}$. Hence, despite no change in the locations of contour
$\mathcal{C}_{P}^{j}$ ($j=1,2,3$), now only contour $\mathcal{C}_{P}^{3}$ is
conservative. Contours $\mathcal{C}_{P}^{1}$ and $\mathcal{C}_{P}^{2}$ are
clearly non-conservative, as confirmed in Table 3.

Contour | $\zeta(P,\mathcal{C}_{P})$
---|---
Structure A | Structure B | Structure C
$\mathcal{C}_{P}^{1}$ | 0.353 | 0.475 | -0.259
$\mathcal{C}_{P}^{2}$ | -0.125 | -0.016 | -0.539
$\mathcal{C}_{P}^{3}$ | 0.176 | 0.200 | 0.042
Table 3: Estimates of $\zeta(P,\mathcal{C}_{P})$ from (29) for IFORM contours
corresponding to the three best fitting models for $H_{S}$, $S_{2}$ from Table
1, for each of structures A, B and C from Table 2. Figure 5: Density of storm
peak sea state parameters $(H_{S},S_{2})$ conditioned on observing the
1000-year marginal response on stick-type structures. Overlaid are the three
highest scoring IFORM fits by AS, with contour $\mathcal{C}^{i}_{P}$ labelled
C$i$ for $i=1\ldots,3$. The colouring of the contours (red, orange, grey)
indicates the order of ranking, in terms of decreasing predictive performance.
Figure 6: Estimated log probability of exceeding the $P$-year response for
$P=1000$ within a 3-hour sea state as a function of $H_{S}$ and $S_{2}$, for
structures A, B and C (left to right), obtained using the foward method of
Section 3. Overlaid are IFORM contours for the three best-fitting
environmental models, with contour $\mathcal{C}^{i}_{P}$ labelled C$i$ for
$i=1\ldots,3$. The original sample of storm peak $(H_{S},S_{2})$ is shown as
grey dots.

Figure 6 shows the estimated log probability of exceeding the $P$-year
response 3-hour sea state, given by
$\log\\{1-F_{R_{L}|\mathbf{X}}(r_{p}|\mathbf{x})\\}$ where
$F_{R_{L}|\mathbf{X}}$ is estimated as in Section 3, for $P=1000$ and for all
$\mathbf{X}=(H_{S},S_{2})$ within a subset $[0,25]\times[0,0.08]$ of the
environment space $\mathbb{R}^{2}$. IFORM contours using the highest ranking
environmental model are overlaid. The panels show that the probability of
exceeding $P$-year response for given environmental conditions varies with
structure type, resulting in differing locations of a ‘frontier’ where this
probability becomes non-zero (i.e., where the log probability exceeds roughly
-30), corresponding to the shaded regions on the plots. For structure A, the
frontier lies along constant $S_{2}$, indicating that $H_{S}$ alone effects
the base shear induced on the structure. For structures B and C, we see convex
and concave curvature of the frontier respectively, indicating that these
structures are more susceptible to high- and low-steepness conditions
respectively. These differences in the underlying wave-structure interaction
are reflected in the positions of $\text{CDE}_{P}$ in Figures 4 and 5, but not
in the locations of the IFORM contours which remain the same across all
structures. This change in frontier location is therefore the cause of
different contour performances across structures, seen in Table 3 and Figure
5. For instance, contours $\mathcal{C}_{P}^{1}$ and $\mathcal{C}_{P}^{2}$ give
non-conservative estimates for the $P$-year response on structure C because
they fail to account for the concave curvature of the frontier see in the
third panel of Figure 6, and so not react to the increased severity of
response seen at lower values of $S_{2}$.

##  6 Discussion

Offshore structures are subjected to extreme environmental conditions (in the
current context, significant wave height $H_{S}$ and wave steepness $S_{2}$),
making structural risk assessment a critical step in the design process. Ideal
methods combine models for environmental extremes with direct estimation of
the environment-structure interaction from fluid loading. Historically,
however, these approaches have been computationally infeasible; instead, met-
ocean engineers use approximate methods relying on characterisation of the
joint environment, in terms of contours, with accompanying assumptions about
the nature of the environment-structure interaction (see Ross et al. 2020;
Haselsteiner et al. 2021 for a recent review).

We demonstrate an efficient approach to estimating the distribution of base
shear on three stick structures, in a central North Sea environment.
Environmental modelling is carried out using a combination of asymptotically
justified models for univariate (Davison and Smith, 1990) and conditional
extremes (Heffernan and Tawn, 2004). Simulation of short term conditional wave
fields then exploits the conditional method of Taylor et al. (1997) with
efficient sampling (Towe et al., 2023) to estimate the distribution of Morison
base shear (see Morison et al. 1950). Findings are summarised in terms of the
conditional density of the environment ($\text{CDE}_{P}$), the conditional
joint density of environmental variables given occurrence of the $P$-year
response. We use $\text{CDE}_{P}$ to evaluate the usefulness of corresponding
$P$-year IFORM contours (see e.g., Winterstein et al. 1993), estimated for
environmental models for $H_{S}$ and $S_{2}|H_{S}$ with good predictive
performance. Results highlight two deficiencies of the IFORM method, and hence
of design contours in general. First, although design contours are intended to
be conservative by construction, they are unable to reflect the type of
structure under consideration and so may perform well for one structure but
not for another. Second, identification of a good environmental model to
construct the contour is challenging, and there is sufficient variability in
the location of the contour due to the choice of environmental model to effect
its performance when compared to $\text{CDE}_{P}$; current practice tends to
ignore this source of uncertainty. In addition, IFORM represents one of a
number of possible variants of environmental design contours (see Ross et al.
2020; Haselsteiner et al. 2021; Mackay and Haselsteiner 2021; Mackay and de
Hauteclocque 2023); it is often not clear which contour is most appropriate in
a given application.

We employ techniques for selecting the extreme thresholds in the models of
Davison and Smith (1990) and Heffernan and Tawn (2004) that can lead to
subjective choices for each. Recent work by Murphy et al. (2023) provides an
automated approach to threshold selection that eliminates this subjectivity in
the univariate case. Additional methods for handling the choice of
conditioning threshold for the model of Heffernan and Tawn (2004) exist, such
as testing for the independence between exceedance and residual values
(discussed by Jonathan et al. 2012) or bootstrap sampling to quantify the
uncertainty in parameter estimates due to threshold choice (see e.g., Jonathan
et al. 2010). Future analysis could therefore be improved by developing
methods to automate the choice of conditional model threshold which
incorporate the aforementioned techniques, to be used alongside the univariate
method of Murphy et al. (2023).

In this work, simple stick structures are considered, with a response
dependent on only two environmental variables ($H_{S}$ and $S_{2}$). We
believe that this framework provides sufficiently realistic wave surface and
kinematic models for our structure types. Despite this, future analysis might
benefit from the inclusion of more complex structure models and wave-structure
interactions. In reality, there are factors such as wind and current,
alongside additional directional and seasonal covariate effects, not
accommodated in this work. Structure models which include effects such as
local loading and wave breaking could also be utilised, alongside improved
models for the environment itself. For example, linear wave theory may be
extended by transforming linear wave characteristics to their respective non-
linear equivalents following the approach outlined in Swan (2020) and Gibson
(2020). Use of more complex structure and environment models, however, incurs
higher computational cost, and so an approach for efficient estimation of
$\text{CDE}_{P}$ that avoids the need for numerical simulation from fluid
loading models may be desirable. Moreover, adoption of methods of full
probabilistic structural design must be undertaken with care, to ensure
rational evidence-based evolution of design procedures. For example, Standard
Norge (2022) identifies that a number of the features of the methodology of
the LOADS joint industry project (Swan, 2020; Gibson, 2020) are either not
compatible with the NORSOC standard, or not yet sufficiently well stress-
tested for adoption within the standard. From this perspective, contours
retain the advantages of being less computationally costly to employ, whilst
also only requiring knowledge of the joint environment.

The current analysis illustrates the benefits of full probabilistic structural
analysis relative to environmental contours, and highlights the importance of
principled environmental model selection for the application of environmental
contours. Wherever possible, we recommend the application of full
probabilistic structural design.

## Acknowledgements

The work was completed while Matthew Speers was part of the EPSRC funded
STOR-i centre for doctoral training (grant no. EP/S022252/1), with part-
funding from the ARC TIDE Industrial Transformational Research Hub at the
University of Western Australia. The authors wish to acknowledge the support
of colleagues at Lancaster University and Shell.

## References

  * Airy (1845) Airy, G.B., 1845.  Tides and waves.  B. Fellowes.
  * Chai and Leira (2018) Chai, W., Leira, B.J., 2018.  Environmental contours based on inverse SORM.  Marine Structures 60, 34–51.
  * Chavez-Demoulin and Davison (2005) Chavez-Demoulin, V., Davison, A.C., 2005.  Generalized additive modelling of sample extremes.  Journal of the Royal Statistical Society Series C: Applied Statistics 54, 207–222.
  * Coles (2001) Coles, S.G., 2001.  An Introduction to Statistical Modeling of Extreme Values. volume 208.  Springer.
  * Coles and Tawn (1994) Coles, S.G., Tawn, J.A., 1994.  Statistical methods for multivariate extremes: an application to structural design (with discussion).  Journal of the Royal Statistical Society: Series C (Applied Statistics) 43, 1–31.
  * Cooley and Tukey (1965) Cooley, J.W., Tukey, J.W., 1965.  An algorithm for the machine calculation of complex Fourier series.  Mathematics of Computation 19, 297–301.
  * Davison and Smith (1990) Davison, A.C., Smith, R.L., 1990.  Models for exceedances over high thresholds (with discussion).  Journal of the Royal Statistical Society Series B 52, 393–425.
  * DNVGL-RP-C205 (2017) DNVGL-RP-C205, 2017.  Environmental conditions and environmental loads.  Det Norske Veritas group, Norway.
  * Ewans and Jonathan (2008) Ewans, K., Jonathan, P., 2008.  The effect of directionality on northern North Sea extreme wave design criteria.  Journal of Offshore Mechanics and Arctic Engineering 130, 041604.
  * Feld et al. (2015) Feld, G., Randell, D., Wu, Y., Ewans, K., Jonathan, P., 2015.  Estimation of storm peak and intra-storm directional-seasonal design conditions in the North Sea.  Journal of Offshore Arctic Engineering 137, 021102.
  * Gelman et al. (2013) Gelman, A., Carlin, J.B., Stern, H.S., Dunson, D.B., Vehtari, A., Rubin, D.B., 2013.  Bayesian Data Analysis.  Chapman and Hall/CRC.
  * Gibson (2020) Gibson, R., 2020.  Extreme Environmental Loading of Fixed Offshore Structures: Summary Report, Component 2.  https://www.hse.gov.uk/offshore/assets/docs/summary-report-component2.pdf.
  * Hafver et al. (2022) Hafver, A., Agrell, C., Vanem, E., 2022.  Environmental contours as Voronoi cells.  Extremes 25, 451–486.
  * Haselsteiner et al. (2021) Haselsteiner, A.F., Coe, R.G., Manuel, L., Chai, W., Leira, B., Clarindo, G., Guedes Soares, C., Hannesdóttir, Á., Dimitrov, N., Sander, A., Ohlendorf, J.H., Thoben, K.D., de Hauteclocque, G., Mackay, E., Jonathan, P., Qiao, C., Myers, A., Rode, A., Hildebrandt, A., Schmidt, B., Vanem, E., Huseby, A.B., 2021.  A benchmarking exercise for environmental contours.  Ocean Engineering 236, 109504.
  * Hasselmann et al. (1973) Hasselmann, K., Barnett, T.P., Bouws, E., Carlson, H., Cartwright, D.E., Enke, K., Ewing, J., Gienapp, A., Hasselmann, D., Kruseman, P., Meerburg, A., Muller, P., Olbers, D., Richter, K., Sell, W., Walden, H., 1973.  Measurements of wind-wave growth and swell decay during the Joint North Sea Wave Project (JONSWAP).  Ergänzungsheft zur Deutschen Hydrographischen Zeitschrift, Reihe A, Nr. 12.
  * Haver (1987) Haver, S., 1987.  On the joint distribution of heights and periods of sea waves.  Ocean Engineering 14, 359–376.
  * Haver and Winterstein (2009) Haver, S., Winterstein, S., 2009.  Environmental contour lines: a method for estimating long term extremes by a short term analysis.  Transactions - Society of Naval Architects and Marine Engineers 116, 116–127.
  * Heffernan and Tawn (2004) Heffernan, J.E., Tawn, J.A., 2004.  A conditional approach for multivariate extreme values (with discussion).  Journal of the Royal Statistical Society Series B 66, 497–546.
  * Holthuijsen (2010) Holthuijsen, L.H., 2010.  Waves in Oceanic and Coastal Waters.  Cambridge University Press.
  * Jonathan et al. (2012) Jonathan, P., Ewans, K., Flynn, J., 2012.  Joint modelling of vertical profiles of large ocean currents.  Ocean Engineering 42, 195–204.
  * Jonathan et al. (2014) Jonathan, P., Ewans, K., Randell, D., 2014.  Non-stationary conditional extremes of northern North Sea storm characteristics.  Environmetrics 25, 172–188.
  * Jonathan et al. (2010) Jonathan, P., Flynn, J., Ewans, K., 2010.  Joint modelling of wave spectral parameters for extreme sea states.  Ocean Engineering 37, 1070–1080.
  * Keef et al. (2013) Keef, C., Papastathopoulos, I., Tawn, J.A., 2013.  Estimation of the conditional distribution of a multivariate variable given that one of its components is large: Additional constraints for the Heffernan and Tawn model.  Journal of Multivariate Analysis 115, 396–404.
  * Mackay and Haselsteiner (2021) Mackay, E., Haselsteiner, A.F., 2021.  Marginal and total exceedance probabilities of environmental contours.  Marine Structures 75, 102863.
  * Mackay and de Hauteclocque (2023) Mackay, E., de Hauteclocque, G., 2023.  Model-free environmental contours in higher dimensions.  Ocean Engineering 273, 113959.
  * Madsen et al. (1986) Madsen, H.O., Krenk, S., Lind, N.C., 1986.  Methods of Structural Safety.  Englewood Cliffs: Prentice-Hall.
  * Morison et al. (1950) Morison, J.R., O’Brien, M.P., Johnson, J.W., Schaaf, S.A., 1950.  The force exerted by surface waves on piles.  Journal of Petroleum Technology 189, 149–154.
  * Murphy et al. (2023) Murphy, C., Tawn, J.A., Varty, Z., 2023.  Automated threshold selection and associated inference uncertainty for univariate extremes.  arXiv preprint arXiv:2310.17999.
  * Myrhaug (2018) Myrhaug, D., 2018.  Some probabilistic properties of deep water wave steepness.  Oceanologia 60, 187–192.
  * NORSOK N-003 (2017) NORSOK N-003, 2017.  NORSOK Standard N-003:2017: Actions and action effects.  NORSOK, Norway.
  * Nussbaumer (1982) Nussbaumer, H.J., 1982.  The fast Fourier transform.  Springer.
  * Randell et al. (2015) Randell, D., Feld, G., Ewans, K., Jonathan, P., 2015.  Distributions of return values for ocean wave characteristics in the South China Sea using directional–seasonal extreme value analysis.  Environmetrics 26, 442–450.
  * Randell et al. (2016) Randell, D., Turnbull, K., Ewans, K., Jonathan, P., 2016.  Bayesian inference for non-stationary marginal extremes.  Environmentrics 27, 439–450.
  * Rosenblatt (1952) Rosenblatt, M., 1952.  Remarks on a multivariate transformation.  The Annals of Mathematical Statistics 23, 470–472.
  * Ross et al. (2020) Ross, E., Astrup, O.C., Bitner-Gregersen, E., Bunn, N., Feld, G., Gouldby, B., Huseby, A., Liu, Y., Randell, D., Vanem, E., 2020.  On environmental contours for marine and coastal design.  Ocean Engineering 195, 106194.
  * Shooter et al. (2021) Shooter, R., Tawn, J.A., Ross, E., Jonathan, P., 2021.  Basin-wide spatial conditional extremes for severe ocean storms.  Extremes 24, 241–265.
  * Standard Norge (2022) Standard Norge, 2022.  Shall NORSOK N-0031 and NORSOK N-0062 be updated as a result of findings in LOADS JIP? Conclusions from the evaluation committee.  https://standard.no/globalassets/fagomrader-sektorer/petroleum/loads-jip-and-norsok-n_003.pdf.
  * Swan (2020) Swan, C., 2020.  Extreme Environmental Loading of Fixed Offshore Structures: Summary Report, Component 1.  https://www.hse.gov.uk/offshore/assets/docs/summary-report-component1.pdf.
  * Taylor et al. (1997) Taylor, P.H., Jonathan, P., Harland, L.A., 1997.  Time domain simulation of jack-up dynamics with the extremes of a Gaussian process.  Journal of Vibration and Acoustics 119, 624–628.
  * Tendijck et al. (2023) Tendijck, S., Eastoe F., E., Tawn, J.A., Randell, D., Jonathan, P., 2023.  Modeling the extremes of bivariate mixture distributions with application to oceanographic data.  Journal of the American Statistical Association 118, 1373–1384.
  * Towe et al. (2023) Towe, R., Randell, D., Kensler, J., Feld, G., Jonathan, P., 2023.  Estimation of associated values from conditional extreme value models.  Ocean Engineering 272, 113808.
  * Towe et al. (2021) Towe, R., Zanini, E., Randell, D., Feld, G., Jonathan, P., 2021.  Efficient estimation of distributional properties of extreme seas from a hierarchical description applied to calculation of un-manning and other weather-related operational windows.  Ocean Engineering 238, 109642.
  * Towe et al. (2019) Towe, R.P., Tawn, J.A., Lamb, R., Sherlock, C.G., 2019.  Model-based inference of conditional extreme value distributions with hydrological applications.  Environmetrics 30, e2575.
  * Tromans and Vanderschuren (1995) Tromans, P.S., Vanderschuren, L., 1995.  Response based design conditions in the North Sea: application of a new method, in: Offshore Technology Conference, OTC. pp. OTC–7683.
  * Vanem et al. (2022) Vanem, E., Zhu, T., Babanin, A., 2022.  Statistical modelling of the ocean environment - a review of recent developments in theory and applications.  Marine Structures 86, 103297.
  * Winter and Tawn (2017) Winter, H.C., Tawn, J.A., 2017.  kth-order Markov extremal models for assessing heatwave risks.  Extremes 20, 393–415.
  * Winterstein et al. (1993) Winterstein, S.R., Ude, T.C., Cornell, C.A., Bjerager, P., Haver, S., 1993.  Environmental parameters for extreme response: Inverse form with omission factors, in: Proceedings of the ICOSSAR-93, Innsbruck, Austria, pp. 551–557.
  * Youngman (2019) Youngman, B.D., 2019.  Generalized additive models for exceedances of high thresholds with an application to return level estimation for U.S. wind gusts.  Journal of the American Statistical Association 114, 1865–1879.

Supplementary Material to ‘Estimating Metocean Environments Associated with
Extreme Structural Response’

## Overview

Here we present further details on aspects of the methodology covered in the
main text. Section S1 describes the diagnostic techniques used to select the
marginal exceedance thresholds for the generalised Pareto distribution (GPD)
tail model (4) of Section 3.2.2. In Section S2.1, we provide further details
on the origin of the wave kinematics equations given in Section 3.3.2, as well
as information on the efficient simulation of wave fields from this model.
Section S3 gives the specific parametric form of the models considered for
$S_{2}|H_{S}$ and $S_{2}^{-}|H_{S}$ in Section 4.3, as well as the complete
results summarised by Table 1 in Section 5.2.

##  S1 Univariate threshold selection for $H_{S}$ modelling

The suitable threshold for the modelling of $H_{S}$ and $S_{2}$ via (4) in
Section 3.2.2 is selected using standard extreme value diagnostics, such as
stability plots and mean-residual-life plots (see Coles 2001 for examples).
These methods are summarised below.

First, we consider the choice of threshold $u_{H_{S}}$ for $H_{S}$. Figure S1
shows the values of the GPD standardised scale ${\sigma}^{*}=\sigma_{u}-\xi u$
(for $\sigma_{u}$ estimated with threshold $u>0$) and shape parameter $\xi$
for choices of threshold non-exceedance probability
$p_{u_{H_{S}}}=\tilde{F}_{H_{S}}(u_{H_{S}})$, alongside their respective 95%
confidence intervals, obtained using block bootstrapping. The estimates of
each parameter appear stable beyond the 0.8th percentile. We thus select
$u_{H_{S}}=\tilde{F}_{H_{S}}^{-1}(0.8)$ as the exceedance threshold for the
marginal modelling of $H_{S}$. We select threshold $u_{S_{2}}$ for marginal
modelling of $S_{2}$ using the same approach. Figure S2 shows the equivalent
plots obtained when fitting to $S_{2}$. Again, we see stability of parameter
estimates for thresholds with exceedance probability past 0.8, and so we
select $u_{S_{2}}=\tilde{F}_{S_{2}}^{-1}(0.8)$ as the exceedance threshold for
the marginal modelling of $S_{2}$.

We verify these threshold choices using mean residual life plots (see Coles
2001 for details). The form of the GPD implies that a linear trend in mean
excesses with respect to threshold will occur above thresholds satisfying the
model conditions. Figure S3 shows that there is a linear trend in mean excess
point estimates when conditioning on thresholds above the 0.8th percentile for
both $H_{S}$ and $S_{2}$.

Figure S1: Stability plots for estimated values of the GPD scale $\sigma$ and
shape $\xi$ parameters when fitting to $H_{S}$ peak data, with respect to
varying levels of the conditioning threshold percentile. Point estimates
obtained via maximum likelihood estimation are marked in black. Corresponding
95% confidence intervals are shown in red. Figure S2: Stability plots for
estimated values of the GPD scale $\sigma$ and shape $\xi$ parameters when
fitting to $S_{S}$ peak data, with respect to varying levels of the
conditioning threshold percentile. Point estimates obtained via maximum
likelihood estimation are marked in black. Corresponding 95% confidence
intervals are shown in red. Figure S3: Mean residual life (mean excess) plots
for excesses of $H_{S}$ and $S_{2}$. Point estimates and 95% confidence
intervals of mean excess with given exceedance probability are shown. The GPD
implies a linear trend in mean excesses with respect to threshold. Thus,
exceedance probabilities for which this linear trend is observed provide valid
exceedance thresholds.

##  S2 Wave field simulation

###  S2.1 Airy wave theory

Section 3.3.2 introduces the wave model of Taylor et al. (1997) which is
itself derived from the work of Airy (1845), commonly referred to as linear
wave theory, which provides physically based models for the surface elevation
of the wave surface and associated kinematics. Airy (1845) models the
stochastic surface elevation $E$ at time $t\in\mathbb{R}$ and location
$x\in\mathbb{R}$ as

$E(t;x)=\sum_{n=1}^{N}\left\\{A_{n}\cos\left(\omega_{n}t-k_{n}x\right)+B_{n}\sin\left(\omega_{n}t-k_{n}x\right)\right\\},$
(S.30)

with contributing angular frequencies $\omega_{n}>0$, coefficients
$A_{n},B_{n}>0$ and wave numbers $k_{n}$ as defined in Section 3.3.2, for
$n=1,\ldots,N$. This is the non-conditioned form of expression (10), for an
arbitrary location $x\in\mathbb{R}$. This model ensures conservation of mass,
momentum satisfaction of necessary boundary conditions in a simple setting;
see Holthuijsen (2010) for details.

Airy (1845) also introduces the ‘velocity potential’ equation

$\phi(t;x,z)=\sum_{n=1}^{N}\left\\{\frac{\omega_{n}\cosh[k_{n}(d+z)]}{k_{n}\sinh(k_{n}d)}\left[-A_{n}\sin(\omega_{n}t-k_{n}x)+B_{n}\cos(\omega_{n}t-k_{n}x)\right]\right\\},$
(S.31)

for time $t\in\mathbb{R}$, horizontal position $x\in\mathbb{R}$ and vertical
position $z\in\mathbb{R}$ relative to mean surface level, with water depth
$d$. As in the main text, all physical quantities are given in SI units.
Expression (S.31) possesses the property

$U(t;x,z)=\frac{\partial\phi}{\partial x}\quad\text{and}\quad
V(t;x,z)=\frac{\partial\phi}{\partial z},$ (S.32)

where $U(t;x,z)$ and $V(t;x,z)$ are horizontal and vertical velocities
respectively. Hence, differentiation of (S.31) yields

$\displaystyle U(t;x,z)$
$\displaystyle=\sum_{n=1}^{N}\left[\left\\{A_{n}\cos(\omega_{n}t-k_{n}x)+B_{n}\sin(\omega_{n}t-k_{n}x)\right\\}\omega_{n}\frac{\cosh\left[k_{n}(d+z)\right]}{\sinh(k_{n}d)}\right],$
(S.33) $\displaystyle V(t;x,z)$
$\displaystyle=\sum_{n=1}^{N}\left[\left\\{-A_{n}\sin(\omega_{n}t-k_{n}x)+B_{n}\cos(\omega_{n}t-k_{n}x)\right\\}\omega_{n}\frac{\sinh\left[k_{n}(d+z)\right]}{\sinh(k_{n}d)}\right],$
(S.34) $\displaystyle\dot{U}(t;x,z)$
$\displaystyle=\sum_{n=1}^{N}\left[\left\\{-A_{n}\sin(\omega_{n}t-k_{n}x)+B_{n}\cos(\omega_{n}t-k_{n}x)\right\\}\omega_{n}^{2}\frac{\cosh\left[k_{n}(d+z)\right]}{\sinh(k_{n}d)}\right],$
(S.35) $\displaystyle\dot{V}(t;x,z)$
$\displaystyle=\sum_{n=1}^{N}\left[\left\\{-A_{n}\cos(\omega_{n}t-k_{n}x)-B_{n}\sin(\omega_{n}t-k_{n}x)\right\\}\omega_{n}^{2}\frac{\cosh\left[k_{n}(d+z)\right]}{\sinh(k_{n}d)}\right],$
(S.36)

for horizontal velocity $U$, horizontal acceleration $\dot{U}$, vertical
velocity $V$ and vertical acceleration $\dot{V}$. Note that in the main text
we only utilise models derived from the above expressions for the horizontal
velocity $U$ and acceleration $\dot{U}$ with $x=0$; the expressions for $V$
and $\dot{V}$ are included here for completeness.

###  S2.2 JONSWAP wave spectrum

In Section 3.3.2, we assume the wave spectrum $S(\omega;\mathbf{X})$ has a
JONSWAP parametric form, which is shown by Hasselmann et al. (1973) to be
suitably flexible to capture measure offshore spectral behaviour. The JONSWAP
spectrum has spectral density function

$S(\omega;\mathbf{X}=\mathbf{x})=\alpha\omega^{-r}\exp\left\\{-\frac{r}{4}\left(\frac{|\omega|}{\omega_{p}(\mathbf{x})}\right)^{-4}\right\\}\gamma\delta(\omega;\mathbf{x}),$

for $\omega>0$, where $\mathbf{X}=(H_{s},S_{2})$ and
$\omega_{p}(\mathbf{x})=2\pi/t(\mathbf{x})$, where $t(\mathbf{x})$ is the
observed value of the second spectral moment wave period $T_{2}=\sqrt{({2\pi
H_{S}})/({gS_{2}})}$ in sea state $\mathbf{X}=\mathbf{x}$, with

$\delta(\omega;\mathbf{x})=\exp\left\\{-\frac{1}{2\left(0.07+0.02\cdot
I\\{{\omega_{p}(\mathbf{x})}>|\omega|\\}\right)^{2}}\left(\frac{|\omega|}{\omega_{p}(\mathbf{x})}-1\right)^{2}\right\\},$

and constants $r,\alpha,\gamma>0$. The normalising constant $\alpha$ is chosen
so that

$4\cdot\left\\{\int_{-\infty}^{\infty}S(\omega;\mathbf{x})\text{d}\omega\right\\}^{\frac{1}{2}}=h(\mathbf{x}),$

where $h(\mathbf{x})$ is the observed value of significant wave height $H_{S}$
in sea state $\mathbf{X}=\mathbf{x}$.

###  S2.3 Efficient wave simulation using the fast Fourier transform

We apply the fast Fourier transform (FFT) algorithm (Nussbaumer, 1982) to
efficiently compute the linear wave behaviour given by expressions (10) and
(11) of Section 3.3.2, by formulating them as the discrete Fourier transforms
of appropriately chosen ‘link functions’. Specifically, for

$A^{\prime}_{n}=A_{n}+\mathcal{Q}\sigma^{2}_{n}\quad\text{and}\quad
B^{\prime}_{n}=B_{n}+\mathcal{R}\sigma^{2}_{n}\omega_{n},$

we write

$\displaystyle E(t_{j})$
$\displaystyle=\sum^{N}_{n=1}g^{(\eta)}_{n}\exp\left\\{{\frac{-2\pi
i}{N}jn}\right\\},\quad g^{(\eta)}_{n}=A^{\prime}_{n}+B^{\prime}_{n}i,$
$\displaystyle U(t_{j},z)$
$\displaystyle=\sum^{N}_{n=1}g^{(u)}_{n}\exp\left\\{{\frac{-2\pi
i}{N}jn}\right\\},\quad
g^{(u)}_{n}=(A^{\prime}_{n}+B^{\prime}_{n}i)\cdot\omega_{n}\cdot\frac{\cosh\\{k_{n}(z+d)\\}}{\sinh(k_{n}d)},$
$\displaystyle\dot{U}(t_{j},z)$
$\displaystyle=\sum^{N}_{n=1}g^{(\dot{u})}_{n}\exp\left\\{{\frac{-2\pi
i}{N}jn}\right\\},\quad
g^{(\dot{u})}_{n}=(B^{\prime}_{n}-A^{\prime}_{n}i)\cdot\omega_{n}^{2}\cdot\frac{\cosh\\{k_{n}(z+d)\\}}{\sinh(k_{n}d)},$

where $i^{2}=-1$ and $j=1,\ldots,N$. This allows us to utilise the fast
Fourier transform (Nussbaumer, 1982) to efficiently compute the linear wave
behaviour given by expressions (10) and (11) for $z\in\mathbb{R}$ at times
$\\{t_{j}\\}_{j=1}^{N}$. Due to the intractability of the above kinematics at
values of $z>0$, we employ simple linear kinematic stretching and evaluate the
above equations at

$z^{\prime}=\begin{cases}z&\text{for}\,\,z\leq 0,\\\
0&\text{for}\,\,z>0,\end{cases}$

rather than $z$.

##  S3 Storm peak hierarchical model selection

###  S3.1 Distributional forms for steepness modelling

Here we give the parametric forms of the four distributions considered for the
modelling of $S_{2}|H_{S}$ and $S^{-}_{2}|H_{S}$ in Section 4.3, described by
the most appropriate of their distribution and density functions. First, is
the Generalised Extreme Value (GEV) distribution, which has distribution
function

$F_{\text{GEV}}(x)=\exp\left\\{-\left[1+\xi\left(\frac{x-\mu_{G}}{\sigma_{G}}\right)\right]^{-1/\xi}\right\\},$
(S.38)

defined on the set $\\{x:1+\xi(x-\mu)/\sigma>0\\}$, with location
$\mu_{G}\in\mathbb{R}$, scale $\sigma_{G}>0$ and shape $\xi\in\mathbb{R}$.
Second, the Weibull distribution, which has distribution function

$F_{\text{Weibull}}(x)=1-e^{-(x/\lambda)^{k}},$ (S.39)

for $x\geq 0$, with scale $\lambda>0$ and shape $k>0$. Third, the Lognormal
distribution, with density function

$f_{\text{Lognormal}}(x)=\frac{1}{x\sigma_{L}\sqrt{2\pi}}\exp\left(-\frac{(\ln
x-\mu_{L})^{2}}{2\sigma_{L}^{2}}\right),$ (S.40)

for $x>0$, with location $\mu_{L}$ and log-scale $\sigma_{L}>0$. Fourth, the
gamma distribution, with density function

$f_{\text{Gamma}}(x)=\frac{\beta^{\alpha}}{\Gamma(\alpha)}x^{\alpha-1}e^{-\beta
x},$ (S.41)

where $\Gamma$ is the gamma function, for $x\geq 0$, with shape $\alpha>0$ and
rate $\beta>0$.

###  S3.2 Detailed cross validation results

Here we present the full details of the results summarised by Table 1 in
Section 5.2, summarised in Table S1 which shows the AIC, cross validation
scores and AS for all 72 candidate models for $S_{2}|H_{S}$ and
$S^{-}_{2}|H_{S}$. The distributional form and function of $H_{S}$ imposed on
distribution parameters are given, as well as individual scores for each of
the assessment criteria. Columns titled ‘$p_{v}$CV$K$’ contain $K$-fold cross
validation scores for an extreme threshold $v=\tilde{F}_{S_{2}}(p_{v})$. The
AS score is obtained by averaging of AIC and cross validation scores. Each
cross validation scenario is repeated for 30 replicates for each data set,
allowing the calculation of standard errors. The highest ranking model per
distribution per scoring method is marked in bold. In the event of a tie, we
opt for the model with the highest AIC.

Distribution | Parameter Forms | AIC | 0.0CV5 | 0.0CV10 | 0.8CV5 | 0.8CV10 | 0.9CV10 | 0.9CV10 | Aggregate Score (AS)
---|---|---|---|---|---|---|---|---|---
$S^{-}_{2}|H_{S}\sim\text{GEV}$ | Exp $\mu_{G}$, Exp $\sigma_{G}$ | 3.931 | 3.931 (0.002) | 3.931 (0.001) | 4.032 (0.014) | 4.037 (0.004) | 4.024 (0.002) | 4.11 (0.002) | 3.999 (0.002)
$S^{-}_{2}|H_{S}\sim\text{GEV}$ | Exp $\mu_{G}$, Lin $\sigma_{G}$ | 3.930 | 3.937 (0.004) | 3.936 (0.003) | 3.945 (0.027) | 3.959 (0.014) | 3.998 (0.004) | 4.084 (0.005) | 3.97 (0.005)
$S^{-}_{2}|H_{S}\sim\text{GEV}$ | Exp $\mu_{G}$, Qua $\sigma_{G}$ | 3.927 | 3.93 (0.001) | 3.93 (0.001) | 4.034 (0.007) | 4.041 (0.003) | 4.025 (0.005) | 4.106 (0.004) | 3.999 (0.002)
$S^{-}_{2}|H_{S}\sim\text{GEV}$ | Lin $\mu_{G}$, Exp $\sigma_{G}$ | 3.918 | 3.919 (0.001) | 3.919 (0.001) | 4.012 (0.006) | 4.014 (0.005) | 3.975 (0.001) | 4.061 (0.001) | 3.974 (0.004)
$S^{-}_{2}|H_{S}\sim\text{GEV}$ | Lin $\mu_{G}$, Lin $\sigma_{G}$ | 3.913 | 3.926 (0.004) | 3.925 (0.003) | 3.953 (0.017) | 3.951 (0.01) | 3.975 (0.005) | 4.062 (0.001) | 3.958 (0.003)
$S^{-}_{2}|H_{S}\sim\text{GEV}$ | Lin $\mu_{G}$, Qua $\sigma_{G}$ | 3.915 | 3.917 (0.001) | 3.916 (0.001) | 4.015 (0.007) | 4.016 (0.003) | 3.98 (0.004) | 4.06 (0.003) | 3.974 (0.001)
$S^{-}_{2}|H_{S}\sim\text{GEV}$ | Qua $\mu_{G}$, Exp $\sigma_{G}$ | 3.929 | 3.929 (0.004) | 3.93 (0.001) | 4.026 (0.012) | 4.029 (0.01) | 4.011 (0.013) | 4.101 (0.01) | 3.994 (0.003)
$S^{-}_{2}|H_{S}\sim\text{GEV}$ | Qua $\mu_{G}$, Lin $\sigma_{G}$ | 3.927 | 3.934 (0.004) | 3.933 (0.004) | 3.955 (0.016) | 3.97 (0.016) | 3.992 (0.012) | 4.086 (0.008) | 3.971 (0.004)
$S^{-}_{2}|H_{S}\sim\text{GEV}$ | Qua $\mu_{G}$, Qua $\sigma_{G}$ | 3.926 | 3.927 (0.003) | 3.927 (0.002) | 4.028 (0.012) | 4.043 (0.008) | 4.002 (0.018) | 4.105 (0.007) | 3.994 (0.004)
$S^{-}_{2}|H_{S}\sim\text{Weibull}$ | Exp $k$, Exp $\lambda$ | 3.745 | 3.744 (0.001) | 3.744 (0.001) | 3.99 (0.003) | 3.99 (0.002) | 3.992 (0.008) | 4.075 (0.007) | 3.897 (0.002)
$S^{-}_{2}|H_{S}\sim\text{Weibull}$ | Exp $k$, Lin $\lambda$ | 3.679 | 3.666 (0.017) | 3.672 (0.006) | 3.647 (0.258) | 3.843 (0.035) | 3.656 (0.238) | 3.809 (0.088) | 3.71 (0.047)
$S^{-}_{2}|H_{S}\sim\text{Weibull}$ | Exp $k$, Qua $\lambda$ | 3.711 | 3.726 (0.005) | 3.721 (0.003) | 3.955 (0.01) | 3.956 (0.009) | 3.942 (0.013) | 4.019 (0.015) | 3.861 (0.004)
$S^{-}_{2}|H_{S}\sim\text{Weibull}$ | Lin $k$, Exp $\lambda$ | 3.743 | 3.743 (0.002) | 3.743 (0.001) | 3.982 (0.009) | 3.987 (0.004) | 3.993 (0.008) | 4.081 (0.006) | 3.896 (0.002)
$S^{-}_{2}|H_{S}\sim\text{Weibull}$ | Lin $k$, Lin $\lambda$ | 3.578 | 3.584 (0.017) | 3.576 (0.02) | 3.828 (0.017) | 3.855 (0.008) | 3.829 (0.014) | 3.852 (0.015) | 3.729 (0.004)
$S^{-}_{2}|H_{S}\sim\text{Weibull}$ | Lin $k$, Qua $\lambda$ | 3.697 | 3.695 (0.002) | 3.695 (0.002) | 3.885 (0.006) | 3.891 (0.007) | 3.807 (0.014) | 3.897 (0.015) | 3.795 (0.003)
$S^{-}_{2}|H_{S}\sim\text{Weibull}$ | Qua $k$, Exp $\lambda$ | 3.732 | 3.734 (0.007) | 3.735 (0.005) | 3.944 (0.016) | 3.975 (0.009) | 3.95 (0.021) | 4.061 (0.012) | 3.876 (0.004)
$S^{-}_{2}|H_{S}\sim\text{Weibull}$ | Qua $k$, Lin $\lambda$ | 3.718 | 3.715 (0.003) | 3.716 (0.002) | 3.905 (0.011) | 3.909 (0.009) | 3.859 (0.021) | 3.948 (0.021) | 3.824 (0.005)
$S^{-}_{2}|H_{S}\sim\text{Weibull}$ | Qua $k$, Qua $\lambda$ | 3.712 | 3.697 (0.014) | 3.691 (0.014) | 3.87 (0.034) | 3.852 (0.028) | 3.832 (0.067) | 3.944 (0.027) | 3.8 (0.013)
$S^{-}_{2}|H_{S}\sim\text{Lognormal}$ | Exp $\mu_{L}$, Exp $\sigma_{L}$ | 3.887 | 3.855 (0.04) | 3.864 (0.026) | 3.765 (0.118) | 3.924 (0.054) | 3.664 (0.177) | 3.987 (0.108) | 3.849 (0.037)
$S^{-}_{2}|H_{S}\sim\text{Lognormal}$ | Exp $\mu_{L}$, Lin $\sigma_{L}$ | 3.889 | 3.842 (0.04) | 3.802 (0.049) | 3.94 (0.039) | 3.928 (0.042) | 3.93 (0.076) | 3.909 (0.066) | 3.891 (0.018)
$S^{-}_{2}|H_{S}\sim\text{Lognormal}$ | Exp $\mu_{L}$, Qua $\sigma_{L}$ | 3.891 | 3.886 (0.021) | 3.877 (0.029) | 3.976 (0.006) | 3.979 (0.011) | 3.997 (0.003) | 4.082 (0.004) | 3.955 (0.007)
$S^{-}_{2}|H_{S}\sim\text{Lognormal}$ | Lin $\mu_{L}$, Exp $\sigma_{L}$ | 3.894 | 3.893 (0.002) | 3.893 (0.002) | 3.958 (0.008) | 3.964 (0.005) | 3.967 (0.01) | 4.052 (0.006) | 3.946 (0.002)
$S^{-}_{2}|H_{S}\sim\text{Lognormal}$ | Lin $\mu_{L}$, Lin $\sigma_{L}$ | 3.896 | 3.897 (0.001) | 3.897 (0.001) | 3.981 (0.004) | 3.983 (0.003) | 3.999 (0.003) | 4.085 (0.002) | 3.963 (0.001)
$S^{-}_{2}|H_{S}\sim\text{Lognormal}$ | Lin $\mu_{L}$, Qua $\sigma_{L}$ | 3.857 | 3.873 (0.01) | 3.868 (0.008) | 3.932 (0.026) | 3.911 (0.015) | 3.874 (0.038) | 3.935 (0.032) | 3.893 (0.010)
$S^{-}_{2}|H_{S}\sim\text{Lognormal}$ | Qua $\mu_{L}$, Exp $\sigma_{L}$ | 3.880 | 3.862 (0.037) | 3.874 (0.017) | 3.884 (0.063) | 3.843 (0.054) | 3.836 (0.057) | 3.875 (0.075) | 3.865 (0.019)
$S^{-}_{2}|H_{S}\sim\text{Lognormal}$ | Qua $\mu_{L}$, Lin $\sigma_{L}$ | 3.891 | 3.892 (0.003) | 3.892 (0.001) | 3.97 (0.004) | 3.972 (0.003) | 3.976 (0.004) | 4.063 (0.002) | 3.951 (0.001)
$S^{-}_{2}|H_{S}\sim\text{Lognormal}$ | Qua $\mu_{L}$, Qua $\sigma_{L}$ | 3.881 | 3.879 (0.008) | 3.877 (0.006) | 3.925 (0.032) | 3.919 (0.035) | 3.95 (0.037) | 4.03 (0.019) | 3.923 (0.010)
$S^{-}_{2}|H_{S}\sim\text{Gamma}$ | Exp $\alpha$, Exp $\beta$ | 3.801 | 3.858 (0.009) | 3.854 (0.009) | 4.016 (0.023) | 3.998 (0.014) | 3.917 (0.035) | 4.079 (0.011) | 3.932 (0.007)
$S^{-}_{2}|H_{S}\sim\text{Gamma}$ | Exp $\alpha$, Lin $\beta$ | 3.775 | 3.815 (0.021) | 3.821 (0.016) | 3.928 (0.021) | 3.923 (0.015) | 3.832 (0.027) | 3.921 (0.037) | 3.859 (0.007)
$S^{-}_{2}|H_{S}\sim\text{Gamma}$ | Exp $\alpha$, Qua $\beta$ | 3.876 | 3.864 (0.006) | 3.863 (0.004) | 3.971 (0.02) | 3.965 (0.015) | 3.915 (0.023) | 3.982 (0.02) | 3.919 (0.005)
$S^{-}_{2}|H_{S}\sim\text{Gamma}$ | Lin $\alpha$, Exp $\beta$ | 3.872 | 3.872 (0.002) | 3.873 (0.001) | 4.02 (0.005) | 4.021 (0.005) | 3.991 (0.007) | 4.08 (0.006) | 3.961 (0.002)
$S^{-}_{2}|H_{S}\sim\text{Gamma}$ | Lin $\alpha$, Lin $\beta$ | 3.871 | 3.872 (0.002) | 3.873 (0.001) | 4.008 (0.005) | 4.01 (0.008) | 3.974 (0.006) | 4.06 (0.005) | 3.953 (0.002)
$S^{-}_{2}|H_{S}\sim\text{Gamma}$ | Lin $\alpha$, Qua $\beta$ | 3.866 | 3.868 (0.001) | 3.868 (0) | 4.001 (0.005) | 4.004 (0.003) | 3.959 (0.004) | 4.046 (0.004) | 3.945 (0.001)
$S^{-}_{2}|H_{S}\sim\text{Gamma}$ | Qua $\alpha$, Exp $\beta$ | 3.823 | 3.826 (0.004) | 3.825 (0.003) | 3.94 (0.052) | 3.947 (0.01) | 3.881 (0.022) | 3.955 (0.016) | 3.885 (0.009)
$S^{-}_{2}|H_{S}\sim\text{Gamma}$ | Qua $\alpha$, Lin $\beta$ | 3.875 | 3.866 (0.007) | 3.866 (0.006) | 4.017 (0.019) | 4.019 (0.016) | 4.015 (0.017) | 4.081 (0.025) | 3.963 (0.006)
$S^{-}_{2}|H_{S}\sim\text{Gamma}$ | Qua $\alpha$, Qua $\beta$ | 3.837 | 3.844 (0.007) | 3.845 (0.006) | 3.876 (0.024) | 3.918 (0.026) | 3.823 (0.051) | 3.9 (0.046) | 3.863 (0.012)
$S_{2}|H_{S}\sim\text{GEV}$ | Exp $\mu_{G}$, Exp $\sigma_{G}$ | 3.567 | 3.446 (0.035) | 3.479 (0.021) | 3.822 (0.044) | 3.862 (0.023) | 3.919 (0.01) | 3.973 (0.009) | 3.724 (0.010)
$S_{2}|H_{S}\sim\text{GEV}$ | Exp $\mu_{G}$, Lin $\sigma_{G}$ | 3.495 | 3.399 (0.023) | 3.42 (0.013) | 3.586 (0.026) | 3.561 (0.031) | 3.677 (0.013) | 3.734 (0.009) | 3.553 (0.007)
$S_{2}|H_{S}\sim\text{GEV}$ | Exp $\mu_{G}$, Qua $\sigma_{G}$ | 3.555 | 3.532 (0) | 3.531 (0) | 3.82 (0.001) | 3.816 (0.001) | 3.875 (0.001) | 3.954 (0) | 3.726 (0.000)
$S_{2}|H_{S}\sim\text{GEV}$ | Lin $\mu_{G}$, Exp $\sigma_{G}$ | 3.370 | 3.305 (0.016) | 3.316 (0.012) | 3.537 (0.105) | 3.25 (0.075) | 3.025 (0.178) | 2.674 (0.088) | 3.211 (0.032)
$S_{2}|H_{S}\sim\text{GEV}$ | Lin $\mu_{G}$, Lin $\sigma_{G}$ | 3.655 | 3.551 (0.009) | 3.555 (0.013) | 3.645 (0.052) | 3.619 (0.041) | 3.29 (0.172) | 3.482 (0.032) | 3.542 (0.028)
$S_{2}|H_{S}\sim\text{GEV}$ | Lin $\mu_{G}$, Qua $\sigma_{G}$ | 3.688 | 3.697 (0.002) | 3.694 (0.001) | 3.558 (0.003) | 3.554 (0.002) | 3.237 (0.004) | 3.305 (0.002) | 3.533 (0.001)
$S_{2}|H_{S}\sim\text{GEV}$ | Qua $\mu_{G}$, Exp $\sigma_{G}$ | 3.638 | 3.436 (0.067) | 3.567 (0.008) | 3.67 (0.019) | 3.673 (0.017) | 3.575 (0.032) | 3.586 (0.024) | 3.592 (0.012)
$S_{2}|H_{S}\sim\text{GEV}$ | Qua $\mu_{G}$, Lin $\sigma_{G}$ | 3.287 | 3.104 (0.213) | 3.256 (0.089) | 1.966 (0.087) | 1.686 (0.037) | 1.344 (0.017) | 1.323 (0) | 2.281 (0.038)
$S_{2}|H_{S}\sim\text{GEV}$ | Qua $\mu_{G}$, Qua $\sigma_{G}$ | 3.356 | 2.89 (0.129) | 3.094 (0.06) | 2.365 (0.205) | 2.535 (0.122) | 1.979 (0.149) | 2.567 (0.06) | 2.684 (0.047)
$S_{2}|H_{S}\sim\text{Weibull}$ | Exp $k$, Exp $\lambda$ | 3.864 | 3.867 (0.002) | 3.868 (0.002) | 3.966 (0.007) | 3.979 (0.005) | 3.87 (0.011) | 3.978 (0.013) | 3.913 (0.003)
$S_{2}|H_{S}\sim\text{Weibull}$ | Exp $k$, Lin $\lambda$ | 3.880 | 3.885 (0.002) | 3.885 (0.002) | 4.022 (0.006) | 4.024 (0.005) | 3.961 (0.009) | 4.05 (0.007) | 3.958 (0.002)
$S_{2}|H_{S}\sim\text{Weibull}$ | Exp $k$, Qua $\lambda$ | 3.897 | 3.9 (0.002) | 3.9 (0.001) | 4.04 (0.006) | 4.041 (0.004) | 3.995 (0.02) | 4.088 (0.005) | 3.98 (0.003)
$S_{2}|H_{S}\sim\text{Weibull}$ | Lin $k$, Exp $\lambda$ | 3.872 | 3.873 (0.001) | 3.874 (0.001) | 3.984 (0.005) | 3.989 (0.004) | 3.902 (0.004) | 3.994 (0.006) | 3.927 (0.002)
$S_{2}|H_{S}\sim\text{Weibull}$ | Lin $k$, Lin $\lambda$ | 3.888 | 3.891 (0.003) | 3.891 (0.002) | 4.008 (0.012) | 4.013 (0.009) | 3.96 (0.007) | 4.022 (0.015) | 3.953 (0.003)
$S_{2}|H_{S}\sim\text{Weibull}$ | Lin $k$, Qua $\lambda$ | 3.901 | 3.902 (0.006) | 3.903 (0.003) | 4.032 (0.006) | 4.038 (0.006) | 4.006 (0.007) | 4.096 (0.008) | 3.983 (0.003)
$S_{2}|H_{S}\sim\text{Weibull}$ | Qua $k$, Exp $\lambda$ | 3.871 | 3.872 (0.001) | 3.871 (0.001) | 3.978 (0.009) | 3.986 (0.006) | 3.892 (0.009) | 3.989 (0.006) | 3.923 (0.003)
$S_{2}|H_{S}\sim\text{Weibull}$ | Qua $k$, Lin $\lambda$ | 3.889 | 3.89 (0.002) | 3.89 (0.001) | 4.018 (0.005) | 4.021 (0.002) | 3.957 (0.005) | 4.047 (0.004) | 3.959 (0.001)
$S_{2}|H_{S}\sim\text{Weibull}$ | Qua $k$, Qua $\lambda$ | 3.903 | 3.88 (0.044) | 3.9 (0.01) | 4.046 (0.003) | 4.043 (0.002) | 4.006 (0.003) | 4.096 (0.005) | 3.982 (0.007)
$S_{2}|H_{S}\sim\text{Lognormal}$ | Exp $\mu_{L}$, Exp $\sigma_{L}$ | 3.629 | 3.632 (0.007) | 3.636 (0.005) | 3.872 (0.021) | 3.869 (0.01) | 3.8 (0.022) | 3.891 (0.025) | 3.761 (0.006)
$S_{2}|H_{S}\sim\text{Lognormal}$ | Exp $\mu_{L}$, Lin $\sigma_{L}$ | 3.485 | 3.485 (0.006) | 3.497 (0.006) | 3.803 (0.009) | 3.807 (0) | 3.861 (0) | 3.941 (0) | 3.697 (0.002)
$S_{2}|H_{S}\sim\text{Lognormal}$ | Exp $\mu_{L}$, Qua $\sigma_{L}$ | 3.563 | 3.562 (0.002) | 3.562 (0.001) | 3.896 (0.004) | 3.897 (0.002) | 3.926 (0.006) | 4.008 (0.004) | 3.773 (0.001)
$S_{2}|H_{S}\sim\text{Lognormal}$ | Lin $\mu_{L}$, Exp $\sigma_{L}$ | 3.651 | 3.649 (0.003) | 3.65 (0.002) | 3.914 (0.008) | 3.92 (0.007) | 3.878 (0.011) | 3.961 (0.016) | 3.803 (0.003)
$S_{2}|H_{S}\sim\text{Lognormal}$ | Lin $\mu_{L}$, Lin $\sigma_{L}$ | 3.591 | 3.59 (0.003) | 3.588 (0.003) | 3.69 (0.016) | 3.705 (0.012) | 3.584 (0.027) | 3.7 (0.015) | 3.635 (0.005)
$S_{2}|H_{S}\sim\text{Lognormal}$ | Lin $\mu_{L}$, Qua $\sigma_{L}$ | 3.612 | 3.613 (0.001) | 3.613 (0.001) | 3.747 (0.004) | 3.754 (0.003) | 3.563 (0.004) | 3.646 (0.006) | 3.65 (0.001)
$S_{2}|H_{S}\sim\text{Lognormal}$ | Qua $\mu_{L}$, Exp $\sigma_{L}$ | 3.669 | 3.657 (0.007) | 3.661 (0.004) | 3.942 (0.01) | 3.925 (0.016) | 3.947 (0.013) | 3.97 (0.028) | 3.824 (0.005)
$S_{2}|H_{S}\sim\text{Lognormal}$ | Qua $\mu_{L}$, Lin $\sigma_{L}$ | 3.659 | 3.653 (0.011) | 3.653 (0.01) | 3.895 (0.062) | 3.928 (0.01) | 3.895 (0.055) | 3.901 (0.326) | 3.798 (0.046)
$S_{2}|H_{S}\sim\text{Lognormal}$ | Qua $\mu_{L}$, Qua $\sigma_{L}$ | 3.185 | 3.227 (0.021) | 3.212 (0.016) | 2.642 (0.074) | 2.626 (0.035) | 2.191 (0.099) | 2.37 (0.061) | 2.779 (0.765)
$S_{2}|H_{S}\sim\text{Gamma}$ | Exp $\alpha$, Exp $\beta$ | 3.608 | 3.606 (0.008) | 3.605 (0.005) | 3.694 (0.011) | 3.685 (0.01) | 3.66 (0.017) | 3.773 (0.011) | 3.662 (0.003)
$S_{2}|H_{S}\sim\text{Gamma}$ | Exp $\alpha$, Lin $\beta$ | 3.669 | 3.643 (0.006) | 3.645 (0.005) | 3.795 (0.01) | 3.802 (0.006) | 3.752 (0.013) | 3.843 (0.016) | 3.736 (0.004)
$S_{2}|H_{S}\sim\text{Gamma}$ | Exp $\alpha$, Qua $\beta$ | 3.604 | 3.602 (0.015) | 3.601 (0.013) | 3.974 (0.014) | 3.974 (0.015) | 3.96 (0.016) | 4.038 (0.028) | 3.822 (0.004)
$S_{2}|H_{S}\sim\text{Gamma}$ | Lin $\alpha$, Lin $\beta$ | 3.710 | 3.709 (0.004) | 3.71 (0.002) | 3.975 (0.002) | 3.975 (0.003) | 3.972 (0.007) | 4.058 (0.006) | 3.873 (0.002)
$S_{2}|H_{S}\sim\text{Gamma}$ | Lin $\alpha$, Exp $\beta$ | 3.668 | 3.668 (0.003) | 3.668 (0.003) | 3.85 (0.008) | 3.857 (0.005) | 3.816 (0.007) | 3.902 (0.007) | 3.776 (0.002)
$S_{2}|H_{S}\sim\text{Gamma}$ | Lin $\alpha$, Qua $\beta$ | 3.692 | 3.665 (0.025) | 3.662 (0.024) | 3.91 (0.024) | 3.911 (0.019) | 3.898 (0.045) | 3.995 (0.018) | 3.819 (0.009)
$S_{2}|H_{S}\sim\text{Gamma}$ | Qua $\alpha$, Exp $\beta$ | 3.599 | 3.487 (0.138) | 3.433 (0.112) | 3.904 (0.071) | 3.938 (0.011) | 3.922 (0.01) | 3.987 (0.026) | 3.753 (0.027)
$S_{2}|H_{S}\sim\text{Gamma}$ | Qua $\alpha$, Lin $\beta$ | 3.671 | 3.67 (0.007) | 3.668 (0.005) | 3.908 (0.008) | 3.898 (0.005) | 3.836 (0.015) | 3.993 (0.011) | 3.806 (0.003)
$S_{2}|H_{S}\sim\text{Gamma}$ | Qua $\alpha$, Qua $\beta$ | 3.622 | 3.618 (0.005) | 3.618 (0.004) | 3.985 (0.004) | 3.985 (0.003) | 3.974 (0.005) | 4.045 (0.014) | 3.835 (0.003)

Table S1: AIC and cross validation scores of all model combinations for
$S_{2}|H_{S}$ and $S^{-}_{2}|H_{S}$, for all considered distribution parameter
functional forms. The highest scoring model per distribution per column is
marked in bold. Sample standard errors are given to three decimal places.