Estimating Metocean Environments Associated with Extreme Structural Response
Matthew Speers David Randell Jonathan Angus Tawn Philip Jonathan







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


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.




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

0000 1111


0000 1111








§ 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., ). 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., ) where a deterministic relationship is assumed. The forward approach thus seeks a multivariate distribution F_𝐗 for environmental variables 𝐗, as well as a distribution F_R_L|𝐗 which characterises the maximum stochastic response R_L induced on the structure by the environment 𝐗 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|𝐗 over the environment space. In addition, F_R_L|𝐗 may be used to find the probability of the event of structural failure ℱ_L within a sea state (of duration L) with environmental conditions 𝐗. 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|𝐗 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., ) as well as the conditional multivariate extremes model developed by <cit.>, conditional simulation of extreme time series (of waves and wave kinematics) proposed by <cit.>, and efficient importance sampling from distributions (e.g., ). We demonstrate the direct estimation of the distribution of base shear, approximated using the loading equation of <cit.> on simple structures, following the procedure of <cit.>. 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 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., ); 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 <cit.>, 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 𝐗 (e.g., ). A recent review of statistical methodologies for metocean design is provided by <cit.>.

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_𝐗. 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., ). 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_𝐗 used to obtain it. In this work, we consider the IFORM environmental contour method (e.g., ), recommended by both <cit.> and <cit.> standards. We will assess the relative predictive performance of different hierarchical models for F_𝐗 and their respective IFORM environmental contours. The second objective of this work is to compare estimates for the  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 <ref>, 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 <ref> 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_𝐗 are given in Section <ref>. In Section <ref>, we present estimates of  for three variants of a simple stick structure, and use these to quantify the performance of different IFORM environmental contours. In Section <ref>, 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.




§ 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 <cit.>, 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 <ref> 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.


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π H_S / (g T_2^2), where g is the acceleration due to gravity. Figure <ref> shows the resulting storm peak data 𝐗=(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., ). 







§ ENVIRONMENT AND RESPONSE MODELLING
 



 §.§ Outline of the forward model and CDE_P
 



  §.§.§ 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 CDE_P. For generality, and alignment with the work of others (e.g., ), 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 <ref>, so as to emphasise the key methodological steps only. 

We assume access to metocean data for storm peak variables 𝐗^sp (e.g., H_S and S_2 from Section <ref>). Since the characteristics of storm peak events generally vary with respect to covariates (e.g., ) we also assume access to storm peak covariates Θ^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 ∈𝒮_T = {1, 2,…, 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 {𝐗_s}_s ∈𝒮_T over a given storm (themselves not identically distributed given the sea state covariates {Θ_s}_s ∈𝒮_T) and (c) the resulting maximum response R_S|{(𝐗_s, Θ_s)}_s ∈ S_T induced over sea states within the storm.

As sea states are dependent over time, therefore R_L,s and R_L,s', for any s,s' ∈𝒮_T where s≠ s', 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 {(𝐗_s, Θ_s)}_s ∈𝒮_T evolving over time.  This is because a load R_L,s from the model of <cit.> 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 (𝐗_s, Θ_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', given (𝐗_s, Θ_s) and (𝐗_s', Θ_s'), is negligible. We exploit this reasoning to facilitate step (c), where we assume that the random variables {R_L,s}_s∈𝒮_T are conditionally independent given sea states {(𝐗_s, Θ_s)}_s ∈𝒮_T, which leads to the following simplification, for response r>0


    F_R_S|{(𝐗_s, Θ_s)}_s ∈ S_T(r|{(𝐱_s, θ_s)}_s ∈ S_τ) = ∏_s ∈𝒮_T F_R_L|(𝐗_s, Θ_s) (r|𝐱_s, θ_s),


where we omit the subscript s in writing F_R_L|(𝐗_s, Θ_s) since the dependence of R_L on s is contained through the vales of (𝐗_s, Θ_s).



  §.§.§ 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 𝐗^sp given storm peak covariates Θ^sp, (b) characterisation of the conditional time-varying within-storm evolution of sea state characteristics {(𝐗_s,Θ_s)}_s ∈𝒮_T given storm peak characteristics 𝐗^sp, Θ^sp, and (c) the estimation of the maximum response R_L, s given sea state characteristics (𝐗_s, Θ_s) for an L-hour sea state. We must also consider the variability in the duration T of a storm event. See <cit.> for previous discussions of similar models. Given knowledge of the above, and exploiting assumption (<ref>) made in Section <ref>, the distribution of R_S can be written as 

    F_R_S(r)    =  ∫_θ^sp∫_𝐱^sp∫_({(𝐱_s, θ_s)}_s ∈𝒮_τ, τ)∏_s ∈𝒮_τ F_R_L|(𝐗_s, Θ_s) (r|𝐱_s, θ_s)
       × f_({(𝐗_s, Θ_s)}_s ∈𝒮_T, T) | 𝐗^sp, Θ^sp({(𝐱_s, θ_s)}_s ∈𝒮_τ, τ|𝐱^sp, θ^sp) 
       × f_𝐗^sp|Θ^sp(𝐱^sp|θ^sp) × f_Θ^sp (θ^sp)   ({(𝐱_s, θ_s)}_s ∈𝒮_τ, τ)  𝐱^sp  θ^sp,

for r >0, where f_Θ^sp is the joint probability density of storm peak covariates, f_𝐗^sp|Θ^sp is the joint probability density of storm peak variables given storm peak covariates, and f_({(𝐗_s, Θ_s)}_s∈𝒮_T, T )| 𝐗^sp, Θ^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 λ 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)= ∑_m=0^∞ [F_R_S(r)]^m λ^m e^-λ/m! = exp[-λ(1-F_R_S(r))].


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).



  §.§.§ Reduced forward model
 

The objective of the current work is to compare estimates for  made using the forward approach, summarised by (<ref>), 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 Θ^sp and Θ_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 (<ref>) for F_R_S reduces to 


    F_R_S(r) = ∫_𝐱^sp F_R_L|𝐗^sp(r| 𝐱^sp) f_𝐗^sp(𝐱^sp)  𝐱^sp,
 

for r >0, where F_R_L|𝐗^sp is the cumulative distribution function of the maximum response over an L-hour sea state given storm peak variables 𝐗^sp. For brevity, we henceforth omit the storm peak superscript and write


    F_R_S(r) = ∫_𝐱 F_R_L|𝐗(r| 𝐱) f_𝐗(𝐱)  𝐱,


for r  >0, where f_𝐗 is now the joint density of storm peak variables 𝐗=(X_1,…,X_p), and F_R_L|𝐗 is the distribution of maximum response over an L-hour storm peak sea state with variables 𝐗. This brings us back to the setting of Section <ref>, where we introduced observations of storm peak sea state data 𝐗 with p=2. The key steps in evaluating the reduced forward model in (<ref>) therefore become the estimation of f_𝐗 and F_R_L|𝐗. The first of these inferences is achieved using the conditional extremes model of <cit.>, as described in Section <ref>. The second inference involves conditional simulation of environmental time-series following <cit.>, and importance sampling from the tail of F_R_L|𝐗 following <cit.>, as described in Section <ref>. 



  §.§.§ Probability of structural failure
 

When designing offshore structures, it is often desirable to determine the probability of the event of structural failure ℱ_L for a given L-hour sea state with variables 𝐗=𝐱. 
Given an estimate for the distribution of R_L|(𝐗=𝐱), we can evaluate this probability using the expression

    ℱ_L|𝐗=𝐱 = ∫_r_Lℱ_L|R_L=r_L f_R_L|𝐗 (r_L |𝐱)  r_L,

where the failure probability ℱ_L|R_L=r_L depends on the nature of the structure, as illustrated in Figure <ref> in Section <ref>.



  §.§.§ Conditional density of the environment
 

The conditional density of the environment f_𝐗|R_L(·|r_P) describes the joint density of the environmental variables 𝐗, conditional on the appearance of a P-year maximum response r_P within a sea state of length L hours. Given estimates for f_𝐗, F_R_L|𝐗, F_R_L and r_P,  CDE_P can therefore be evaluated using Bayes' rule 


    f_𝐗|R_L(𝐱|r_P) = f_R_L|𝐗(r_P|𝐱) f_𝐗(𝐱)/f_R_L(r_P),


for 𝐱∈^p, where f_R_L|𝐗 and f_R_L are the densities corresponding to the distributions F_R_L|𝐗 and F_R_L respectively. Examples of  for the central North Sea application are given in Section <ref>.



 §.§ Joint modelling of storm peak conditions
 



  §.§.§ Outline of the conditional extremes model

The upper extremes of the marginal and joint distributions of the environmental variables 𝐗=(X_1,…,X_p), corresponding to the storm peak sea state, are described using the conditional extremes model of <cit.>. For our illustrative example, in Section <ref>, 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 𝐗, 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., ).

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., ). Estimation of the conditional extremes model is thus performed in two stages: (a) marginal extreme value modelling of each variable X_j (j=1,…, p) in turn, followed by the marginal transformation X_j ↦ 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 𝐘=(Y_1, …, Y_p). Subsequently, we estimate f_𝐗 within an extreme joint tail region using the fitted conditional extremes model. The above steps are discussed below. 



  §.§.§ Marginal modelling and marginal transformation to Laplace scale
  

We adopt the approach of <cit.> for marginal modelling of storm peak variable X_j for j=1,…,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) = 
    F_X_j(x)     x ≤ u_j 
    F_X_j(u_j) + {1-F_X_j(u_j)} F_GPD,j(x; u_j, σ_j, ξ_j)     x > u_j,


where F_X_j is the empirical distribution of X_j and

    F_GPD,j(x;u_j, σ_j, ξ_j) = 1 - ( 1+ ξ_j (x-u_j)/σ_j)_+^-1/ξ_j,

for x>u_j, with scale and shape parameters σ_j>0 and ξ_j ∈ℝ and where y_+ = max(y, 0) for y ∈. The values of u_j (j=1,…,p) are selected using the univariate extreme threshold selection methods summarised by <cit.>, see SM. Parameters σ_j and ξ_j are jointly estimated using standard maximum likelihood techniques. The probability integral transform 

    Y_j= log{2 F_X_j(X_j)}    for  X_j<F_X_j^-1(0.5) 
     -log{2[1-F_X_j(X_j)]}    for  X_j > F_X_j^-1(0.5),

is then applied to each variable in turn, obtaining the standard Laplace-scaled equivalents 𝐘 of 𝐗, such that, for j=1,…, p, the distribution of Y_j is

    F_Y_j(y) = 
    1/2exp(y)      y ≤ 0 
     1-1/2exp(-y)      y > 0.




  §.§.§ Joint dependence modelling
 

Having transformed environmental variables 𝐗 to standard Laplace-scale equivalents 𝐘, we now apply the model of <cit.> to estimate the joint distribution F_𝐘 in the upper tail region. It is shown by <cit.> and <cit.> that, for j=1,…, p, there exist unique values of parameters α_|j∈ [-1, 1]^p-1, β_|j∈ (-∞, 1]^p-1, satisfying the constraints of <cit.>, and 𝐳_|j∈^p-1, y>0, such that


    lim_v_j→∞𝐘_-j -α_|j Y_j/Y_j^β_|j < 𝐳_|j, Y_j - v_j >y | Y_j >v_j  = e^-y G_|j(𝐳_|j),


where the (p-1)-vector 𝐘_-j denotes the p-vector 𝐘 with jth element Y_j removed, G_|j is a (p-1)-dimensional distribution function with non-degenerate marginals, and componentwise operations are assumed. Property (<ref>) can be leveraged by assuming a non-linear regression of 𝐘_-j onto Y_j holds for all values of 𝐘 within the region {𝐘∈^p : Y_j > v_j}, for some suitably large finite threshold v_j>0. For conditioning variable Y_j (j=1,…, p), the form of this regression is


    (𝐘_-j|{Y_j=y}) = α_|j y + y^β_|j𝐙_|j ,    y > v_j,


where 𝐙_|j∼ 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 α_|j and β_|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

    𝐙_|j = 𝐘_-j - α_|j Y_j/Y_j^β_|j,   for   Y_j > v_j,

as in <cit.>.



  §.§.§ Simulation under the conditional model and estimation of the environment joint density


Inferences in ^p using the fitted conditional extremes model are typically made by careful combination of Laplace-scale simulations in each of the upper tail regions {𝐘∈^p : Y_j > v_j}, for j=1,…,p, together with empirical estimation in the remaining region {𝐘∈^p : Y_j ≤ v_j  ∀ j}, as described in <cit.>, to give a set of size N_sim realisations from the estimate of the joint density f_𝐘. The p fitted marginal models (<ref>) can then be used componentwise to transform this sample of 𝐘 with Laplace-scale marginals to a sample of 𝐗 on the original physical scale.

We can use the simulated sample to estimate the probability of 𝐗 being in sub-regions of ^p. Specifically, if D is the set of feasible 𝐗 values such that 𝐗∈^p∖ D=0, then we consider a partition (D_1, …, D_M) of D. Then, if N_sim,i is the number of realisations in set D_i, we can estimate 𝐗∈ D_i, for any i=1, … ,M, as (𝐗∈ D_i)=N_sim,i/N_sim. To obtain an estimate f_𝐗(𝐱) for any 𝐱∈ D we exploit the property that, if 𝐱∈ D_i and |D_i| is sufficiently small that f_𝐗 is reasonably constant for all 𝐱∈ D_i, then


    (𝐗∈ D_i) =∫_𝐱'∈ D_if_𝐗(𝐱') 𝐱'≈ |D_i|f_𝐗(𝐱),


yielding the estimate f̂_𝐗(𝐱) = N_sim, i/(N_sim|D_i|) for all 𝐱∈ 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|∝ M^-1. 



 §.§ Estimation of maximum response in a storm peak sea state given storm peak variables
 



  §.§.§ Outline of estimation of F_R_L|𝐗


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 𝐗. However, to estimate the distribution F_R_I|𝐗 of the maximum response due to the action of an individual wave in sea state 𝐗, we first evaluate the distribution F_R_I|(𝐗, C) of the maximum response R_I to an individual wave in the sea state 𝐗 with crest elevation C. This is achieved by simulation of wave fields under sea state conditions 𝐗 with known crest elevation C, followed by propagation of the resulting stochastic wave fields through to the structural response model; see Section <ref> for details. Then, integrating out C, we have

    F_R_I|𝐗(r|𝐱) = ∫_^+ F_R_I|𝐗, C(r|𝐱, c) f_C|𝐗(c|𝐱) c,


for r  >0, where f_C|𝐗 is the density of crest elevation in the sea state 𝐗, where we assume that crests are Rayleigh-distributed, with density


    f_C|𝐗(c|𝐱) = 16c/h(𝐱)^2exp(-8c^2/h(𝐱)^2).


for c>0, with sea state  𝐗=𝐱 with significant wave height h(𝐱)=x_1. Computationally efficient estimation of F_R_I|𝐗 following (<ref>) is achieved using importance sampling; see Section <ref>.  Finally, we obtain F_R_L|𝐗 from F_R_I|𝐗 by assuming that (a) there are a fixed number of waves per L-hour sea state 𝐱, given by Q_L(𝐱) =60^2L/ t_2(𝐱) where t_2(𝐱) 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 <ref> 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|𝐗(r|𝐱) = {F_R_I|𝐗(r|𝐱)}^Q_L(𝐱),

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


  §.§.§ 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|(𝐗, C) in two stages: (a) simulation of realisations of wave fields under sea state conditions 𝐗 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 <cit.> 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 𝐗=𝐱 with wave spectrum S(·; 𝐱), using linear wave theory. The JONSWAP spectrum of <cit.> is chosen as the form of S(·; 𝐱) due to its applicability to the North Sea wave conditions <cit.>, see SM for further details. <cit.> provides expressions for linear crest elevation E(t;𝐱, c), horizontal velocity U(t; z, 𝐱, c), and horizontal acceleration U̇(t; z, 𝐱, c), at time t∈ and vertical position z∈, 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; 𝐱, c)=c, both at the location of the structure. The forms of E(t; 𝐱, c) and U(t;z, 𝐱, c) are


    E(t; 𝐱, c)   =∑_n=1^N {(A_n + 𝒬σ^2_n) cos(ω_n t )+(B_n + ℛσ^2_n ω_n)sin(ω_n t )},


and 

    U(t; z, 𝐱, c) =
                ∑_n=1^N  ω_n cosh(k_n(d+z))/sinh(k_nd){(A_n + 𝒬σ^2_n) cos(ω_n t ) + (B_n + ℛσ^2_n ω_n) sin(ω_n t)},


for z <  E(t; 𝐱, c) and zero otherwise, for a regular grid of angular frequencies ω_1,…, ω_N>0 with spacing δω >0 and N∈ specified below, where A_n,B_n∈ (n=1,…, N) and 𝒬, ℛ∈ are random coefficients, d is water depth and k_n (n=1,…, N) are wave numbers given implicitly by ω_n^2 = g k_n tanh (k_n d). The equation for U̇(t;z, 𝐱, c) can be found by differentiation of U(t;z, 𝐱, c) with respect to t. Coefficients A_n,B_n (n=1,…, N) are a series of independently and identically distributed N(0, σ^2_n) random variables with variance σ^2_n= S(ω_n;𝐱)δω, the integrated spectral density in the frequency band (ω_n - δω/2, ω_n + δω/2) of the discretised wave spectrum. The random coefficients 𝒬 and ℛ are defined as


    𝒬 =1/∑_n σ^2_n( c - ∑_n=1^N A_n )  and  ℛ = 1/∑_n ω_n^2  σ^2_n( - ∑^N_n=1ω_n B_n).


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 <cit.>, 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 𝐗=(X_1,X_2)=(H_S,S_2) from Section <ref> 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, 𝐱, c) induced on a cylinder by the wave at time t and vertical position z, and is given by

    M(t;z, 𝐱, c) =ρ c_m(z) V U̇(t;z, 𝐱, c)+1/2ρ c_d(z) A U(t;z, 𝐱, c)|U(t;z, 𝐱, c)|,

where c_m(z), c_d(z)>0 are inertia and drag coefficients, ρ=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 <ref>. To evaluate the total base shear B_S(t; 𝐱, c) on the structure at time t, we integrate M(t;z, 𝐱, c) to give


    B_S(t; 𝐱, c) = ∫_-d^E(t; 𝐱, c) M(t;z, 𝐱, c) z,


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

The response R_I|(𝐗=𝐱, C=c) may be obtained by considering the portion of the time series {B_S(t;𝐱, c)}_t ∈ that corresponds to the central wave conditioned to attain E(t=0; 𝐱,c)=c; that is, the period of time t ∈𝒯_0⊂, with 0 ∈𝒯_0, for which the wave surrounding the conditioning crest of elevation c at time t=0 acts on the structure. We define

    R_I|(𝐗=𝐱,C=c) = max_t ∈𝒯_0B_S(t;𝐱, c).

We obtain realisations of R_I|(𝐗=𝐱, C=c) from a time series of the base shear response (<ref>) evaluated using Morison loads (<ref>), in turn calculated from wave fields simulated according to expressions (<ref>) and (<ref>). 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 <cit.>, meaning here 𝒯_0 ⊂ [-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 ∈ [-60, 60] and z ∈ [-100, 150]. We set N=n_t in expressions (<ref>) and (<ref>), 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, 𝐱, c)}_i=1,j=1^n_t,n_z. Numerical evaluation of integral (<ref>) with respect to z yields a realisation of the time series {B_S(t_i; 𝐱, c)}_i=1^n_t. A realisation of the maximum individual wave response R_I|(𝐗=𝐱, C=c) is then obtained by applying (<ref>), using the set t ∈{t_i}_i=1^n_t∩𝒯_0 as an approximation to 𝒯_0. Given conditioning crests {c_i}_i=1^k, the above procedure can be used to map c_i ↦ r_i, for i=1,…, k, obtaining a set of maximum responses {r_i}_i=1^k, for a given sea state 𝐱.



  §.§.§ Importance sampling of simulated maximum responses
 

The procedure discussed in Section <ref> is used to obtain realisations {r_i}_i=1^k of R_I|(𝐗, C), for a set of k conditioning crests {c_i}_i=1^k and specified values of storm peak variables 𝐗. These are then used in integral (<ref>) to estimate the distribution F_R_I|𝐗 of the maximum response to an individual wave in sea state 𝐗, assuming a Rayleigh distribution (<ref>) for C|𝐗. However, evaluation of integral (<ref>) via Monte Carlo methods sampling from the Rayleigh density is inefficient in targeting the tail of the response distribution F_R_I|𝐗. Given our interest in the extreme structural response on the structure, we therefore employ the importance sampling approach described by <cit.>, writing integral (<ref>) as


    F_R_I|𝐗 (r|𝐱) = ∫_r F_R_I|𝐗, C(r|𝐱, c)   f_C|𝐗(c|𝐱)/g^(ϵ)_C|𝐗(c|𝐱) g^(ϵ)_C|𝐗(c|𝐱) c,


for r  >0, where g^(ϵ)_C|(𝐗= 𝐱) is the density of the Uniform[0,ϵ h(𝐱)] distribution, for significant wave height h(𝐱) and some ϵ>0. For a fixed number of conditional wave simulations, sampling of the conditioning crest c from g^(ϵ)_C|𝐗 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|𝐗. Therefore, the upper tail of the distribution F_R_I|𝐗 is more efficiently estimated using the sampling distribution g^(ϵ)_C|𝐗. The value of ϵ in (<ref>) is selected so that g^(ϵ)_C|𝐗 provides adequate coverage of the domain of f_C|𝐗, i.e., the exceedance probability


    {C>ϵ h(𝐱)|𝐗=𝐱} = ∫_c>ϵ h(𝐱)f_C|𝐗(c|𝐱)c = exp( -8ϵ^2),


is sufficiently close to zero. We set ϵ=2, which gives a value of probability (<ref>) in the order of 10^-14. 

Integral (<ref>) is then estimated as follows. The k conditional crests {c_i}_i=1^k are sampled from the uniform proposal density g^(ϵ)_C|𝐗. Corresponding realisations  {r_i}_i=1^k of single-wave maximum responses are obtained using the procedure described in Section <ref>.  The distribution F_R_I|𝐗 is then estimated as


    F̂_R_I|𝐗 (r|𝐱) = ∑_i=1^k1_{r_i≤ r}ι^(ϵ) (c_i|𝐱)/∑^k_i=1ι^(ϵ) (c_i|𝐱),


for r >0, where ι^(ϵ)(c|𝐱) = f_C|𝐗(c|𝐱)/g^(ϵ)_C|𝐗(c|𝐱) is the importance sampling ratio and 1_{r_i≤ r}=1 if r_i≤ r, zero otherwise, for i=1,…, k. Estimate (<ref>) is an empirical cumulative distribution function of the simulated responses, weighted to remove bias introduced from sampling crests from g^(ϵ)_C|𝐗 rather than f_C|𝐗. We use estimate (<ref>) to evaluate the distribution of maximum response per L-hour sea state using relation (<ref>). Given an estimate for f_𝐗 obtained as in Section <ref>, we may then calculate the marginal maximum response distribution using integral (<ref>).




§ ENVIRONMENTAL CONTOURS
 



 §.§ 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., ) contours assume a convex form for this boundary, whereas ISORM (e.g., ) 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 <ref>) and the fitting approach we employ to estimate it for our example dataset (Section <ref>). 



 §.§ IFORM design contours
 

Section <ref> details how, given an estimate for the distribution of R_L|𝐗 from Section <ref>, in principle we can evaluate the probability of structural failure ℱ_L|𝐗=𝐱 for a sea state (of duration L) with variables 𝐗=𝐱. IFORM offers an approach to structural design which avoids direct calculation of ℱ_L|𝐗=𝐱, by attempting to make conservative assumptions. For an ocean environment represented by a set of random variables 𝐗 transformed to independent standard Gaussian random variables 𝐔, IFORM assumes that ℱ_L|𝐔=𝐮 is deterministic for all 𝐮, taking values {0,1}, contrary to the failure probability discussed in Section <ref> which takes any value in [0,1]. Writing the region of environmental space corresponding to failure ℱ_L as ℱ, IFORM assumes that the boundary ∂ℱ of ℱ is linear, and lies tangential to a contour of constant transformed environmental density f_𝐔, making the assumed location of ∂ℱ dependent on the joint Gaussian distribution f_𝐔. The assumption of a failure boundary of this type is typically conservative, in that estimates for (𝐔∈ℱ) using it have positive bias. 

The transformation of 𝐗↦𝐔 is achieved via the method of <cit.>, which proceeds as follows. For storm peak variables 𝐗=(X_1, …, 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), …, F_X_p|(X_1, …, 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 <ref> for an example approach. The Rosenblatt transformation maps a realisation 𝐱= (x_1, …, x_p) of 𝐗 to the realisation 𝐮=(u_1, …, u_p) of 𝐔 via u_1 = Φ^-1{F_X_1(x_1)} and u_j = Φ^-1{F_X_j|(X_1, …, X_j-1)(x_j|x_1, …, x_j-1)}, for j=2, …, p, where Φ is the standard Gaussian cumulative distribution function.

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

    β_I = Φ^-1(1-1/N_anP),

where N_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_an P). 
Consider failure region ℱ_β_I with boundary δℱ_β_I tangential to the hypersphere centred at the origin with radius β_I. For any angle ψ, in spherical polar coordinates for a point on the hypersphere, we have

    𝐔∈ℱ_β_I = ∑_i=1^p w_i (ψ) U_i > β_I,

where ∑^p_i=1w_i^2(ψ)=1, due to the linearity of δℱ_β_I and it being tangential to the hypersphere at radius β_I. So, as W=∑_i=1^pw_i(ψ)U_i ∼ N(0,1) given the independence of (U_1, …, U_p), it follows that 𝐔∈ℱ_β_I = W > β_I, which directly gives expression (<ref>) for β_I.

 Once the environmental contour has been estimated in 𝐔-space, it can be represented in the original 𝐗-space via the inverse Rosenblatt transformation x_1 = F^-1_X_1{Φ(u_1)} and  x_j = F^-1_X_j|(X_1, …, X_j-1){Φ (u_j)|x_1, …, x_j-1}, for j=2, …, p. Unlike the hypersphere-shaped contours in the 𝐔-space, these contours are not guaranteed to be convex (see Section <ref>). 

The procedure for construction of a P-year IFORM contour in terms of environmental variable 𝐗=(X_1, …, X_p) is summarised in Algorithm <ref>.






 §.§ Joint parametric models for storm peak variables
 
Construction of the IFORM contour for environmental variables 𝐗= (X_1, X_2) = (H_S, S_2) via Algorithm <ref> 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 <cit.>, we select the GPD tail model (<ref>) 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(μ_L, σ_L), as in <cit.>, Gamma(α, β), Weibull(λ, k) and the Generalised Extreme Value distribution GEV(μ_G, σ_G, ξ); 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 {𝐱_i}_i=1^N_all (see Section <ref>) is partitioned into a `body' where H_S ≤ v and a `tail' where H_S > v for v>0, denoted {𝐱_i^B}_i=1^N_B and {𝐱_i^T}_i=1^N_T respectively, with N_B + N_T = N_all. The tail portion {𝐱_i^T}_i=1^N_T is itself partitioned into K subsets {𝒮_j}_j=1^K each of sizes ⌊ N_T/K⌋ or ⌈ N_T/K⌉. A K-fold cross validation is performed using training set {𝐱_i}_i=1^N ∖𝒮_j and test set 𝒮_j at fold j∈{1, …, 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 (<ref>) (see Section <ref>) 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. 




§ RESULTS
 



 §.§ Estimating the joint density of storm peaks
 

We employ the forward methodology of Sections <ref> and <ref> to estimate the environmental density f_𝐗 and response distribution F_R_L|𝐗 for our motivating dataset of 𝐗=(H_S, S_2) introduced in Section <ref>. In turn, these are used to evaluate the distributions F_R_S and F_R_A (as in Sections <ref> and <ref>) and subsequently  (as in Section <ref>). The  is then compared with various IFORM contour estimates using the methods of Section <ref>.

We model the joint environment 𝐗 using the method of <cit.> discussed in Section <ref>, fitted to data in the region A_ν = {𝐗∈^2 : H_S > ν} for conditioning threshold ν. We take ν = 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_𝐗 is estimated in A_ν as described in Section <ref>. The density in the complement ^2 ∖ A_ν 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.


In the 2-dimensional case 𝐗=(H_S, S_2), the parameters of model (<ref>) reduce to α∈ [-1, 1] and β∈ (-∞, 1] when conditioning on H_S. We find corresponding estimates α̂=0.378 and β̂=0.533, as well as an estimate ξ̂=-0.02 for the H_S marginal shape parameter ξ in (<ref>), the latter indicating a near exponential upper tail for H_S. The effect of α̂ is seen in Figure <ref>, 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 𝒜_ν. 



 §.§ 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 <ref>) each with two distribution parameters modelled as functions of H_S. We do not model the GEV shape ξ 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 <ref>, yielding results given in full in the SM. Table <ref> summarises these results, showing the optimal model for each of S_2|H_s and S_2^-|H_S, with unique a, b, c ∈ 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.




The models in Table <ref> are used with Algorithm <ref> to construct the IFORM contours in Figure <ref>. Each contour corresponds to a return period of P=1000 years.  Contours are labelled 𝒞_P^1 to 𝒞_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  in Section <ref>. 



 §.§ Estimating the conditional density of associated environmental variables
 



  §.§.§ Estimation of  for example structure models


We evaluate  for three examples of the stick structure model (Section <ref>), 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 <ref>. 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∈ [-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.


Figure <ref> shows the corresponding estimates for  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 <ref>.





  §.§.§ Comparison between  and IFORM contours


Provided we have accurate models for the series of p nested conditional distributions F_X_1, …, F_X_p|(X_1, …, 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 <cit.>. 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 , 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 <ref>. 

Assuming that  provides a valid estimate of the environment density conditional on the P-year response, we reason that a well-estimated IFORM contour 𝒞_P should intersect with , and that the discrepancy between 𝒞_P and  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 <ref>, which shows the estimated  for structures A, B and C, for P=1000. Contour 𝒞_P^2 in the left panel overlaps with darker points of  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 <ref> shows 𝒞_P^1 taking (H_S, S_2) values more extreme than the darker points of , 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 𝒞_P^2 taking (H_S, S_2) values less extreme than the darker points of , 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 , 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 <ref> that the overlap between the region of non-zero  and the region 𝒜_P bounded below by the IFORM contour 𝒞_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 , indicating non-conservatism. Conversely, if the overlap includes all of the non-zero region of , the contour appears to be conservative. To formalise this observation, we define the following metric, which quantifies the level of this overlap. Consider


    ζ(P, 𝒜_P) = 2∫_𝐱1_𝒜_P(𝐱) f_𝐗|R_L(𝐱|r_P) 𝐱 - 1,


for f_𝐗|R_L as in Section <ref>, with 1_𝒜_P(𝐱)=1 if 𝐱∈𝒜_P and 1_𝒜_P(𝐱)=0 otherwise, and where r_P is the P-year response of the structure. The metric ζ(P, 𝒞_P) takes values in [-1, 1]. Here, ζ(P, 𝒞_P)>0 indicates conservatism of 𝒞_P due to overestimation of the P-year response, ζ(P, 𝒞_P)< 0 indicates underestimation of the P-year response and so non-conservatism, and ζ(P, 𝒞_P)≈0 indicates accurate estimates for the P-year response and so optimal overall structural design given by IFORM. 

To justify our findings regarding ζ(P, 𝒞_P) in a general case, we apply the following heuristic argument. First, we observe from Figure <ref> that  exhibits approximate marginal symmetry with respect to both H_S and S_2 (i.e., little marginal skew). We also assume 𝒜_P_1⊂𝒜_P_2 for any P_2>P_1≥ 1, and that a contour 𝒞_P varies smoothly with H_S and S_2 for a given P. Then, (a) if ζ(P, 𝒞_P) ≈ 0, the integral over  within 𝒜_P is approximately equal to the integral over  within 𝒜^c_P. We interpret this as indicating that points on 𝒞_P coincide with high values of , as in the cases seen in Figure <ref>, hence structural responses corresponding to points on 𝒞_P will be of similar magnitudes to those from the high density regions of . However, (b) if ζ(P, 𝒞_P)>0, 𝒜_P contains the environmental region where  is non-zero, hence points on 𝒞_P produce structural responses beyond the P-year level, resulting in conservative design using IFORM contour 𝒞_P. Further, (c) if ζ(P, 𝒞_P) <0, the intersection between 𝒜_P and the non-zero  region of the environment space is negligible. Given our assumptions, this arises only when 𝒜_P occupies a region of the environmental space less extreme that that with non-zero . Under these circumstances, points on 𝒞_P will produce structural responses corresponding to return periods less that P, resulting in a lack of conservatism in design.

Estimates for ζ(P,𝒞_P) in Table <ref> support this interpretation, relative to our estimates in Figure <ref>. For structure A, 
only contour 𝒞_P^2 appears to pass through the highest density region of . The other contours lie beyond the highest density region of , corresponding to ζ(P,𝒞_P)>0. Observations for structure B are similar, since  does not vary considerably between structures A and B. For structure C, relative to A and B, the highest density regions in  occur at higher H_S but lower S_2. Hence, despite no change in the locations of contour 𝒞_P^j (j=1,2,3), now only contour 𝒞_P^3 is conservative. Contours 𝒞_P^1 and 𝒞_P^2 are clearly non-conservative, as confirmed in Table <ref>.








Figure <ref> shows the estimated log probability of exceeding the P-year response 3-hour sea state, given by log{1-F_R_L|𝐗(r_p|𝐱) } where F_R_L|𝐗 is estimated as in Section <ref>, for P=1000 and for all 𝐗=(H_S, S_2) within a subset [0, 25] × [0, 0.08] of the environment space ^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  in Figures <ref> and <ref>, 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 <ref> and Figure <ref>. For instance, contours 𝒞_P^1 and 𝒞_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 <ref>, and so not react to the increased severity of response seen at lower values of S_2.




§ 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  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 <cit.> and conditional extremes <cit.>. Simulation of short term conditional wave fields then exploits the conditional method of <cit.> with efficient sampling <cit.> to estimate the distribution of Morison base shear (see ). Findings are summarised in terms of the conditional density of the environment (), the conditional joint density of environmental variables given occurrence of the P-year response. We use  to evaluate the usefulness of corresponding P-year IFORM contours (see e.g., ), 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 ; 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 ); 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 <cit.> and <cit.> that can lead to subjective choices for each. Recent work by <cit.> 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 <cit.> exist, such as testing for the independence between exceedance and residual values (discussed by ) or bootstrap sampling to quantify the uncertainty in parameter estimates due to threshold choice (see e.g., ). 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 <cit.>.

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 <cit.> and <cit.>. Use of more complex structure and environment models, however, incurs higher computational cost, and so an approach for efficient estimation of  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, <cit.> identifies that a number of the features of the methodology of the LOADS joint industry project <cit.> 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.


elsarticle-harv 









    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 <ref> 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 <ref>, 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 <ref> 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.



§ 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  for examples). These methods are summarised below. 

First, we consider the choice of threshold u_H_S for H_S. Figure <ref> shows the values of the GPD standardised scale σ^*=σ_u-ξ u (for σ_u estimated with threshold u>0) and shape parameter ξ for choices of threshold non-exceedance probability p_u_H_S = 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= 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 <ref> 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 = 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  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 <ref> 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.









§ WAVE FIELD SIMULATION




 §.§ Airy wave theory
 

Section 3.3.2 introduces the wave model of <cit.> which is itself derived from the work of <cit.>, commonly referred to as linear wave theory, which provides physically based models for the surface elevation of the wave surface and associated kinematics. <cit.> models the stochastic surface elevation E at time t∈ and location x ∈ as 

    E(t; x)=∑_n=1^N{A_n cos(ω_n t- k_n x)+B_nsin(ω_n t - k_n x)},

 with contributing angular frequencies ω_n>0, coefficients A_n,B_n>0 and wave numbers k_n as defined in Section 3.3.2, for n=1,…, N. This is the non-conditioned form of expression (10), for an arbitrary location x ∈. This model ensures conservation of mass, momentum satisfaction of necessary boundary conditions in a simple setting; see <cit.> for details.
 
 <cit.> also introduces the `velocity potential' equation

    ϕ(t;x, z)=∑_n=1^N{ω_n cosh [k_n(d+z)]/k_nsinh (k_n d)[ -A_nsin (ω_n t- k_n x) + B_ncos (ω_n t- k_n x)]},

for time t ∈, horizontal position x ∈ and vertical position z ∈ relative to mean surface level, with water depth d. As in the main text, all physical quantities are given in SI units. Expression (<ref>) possesses the property 

    U(t; x, z)=∂ϕ/∂ x  and   V(t; x, z) = ∂ϕ/∂ z,

where U(t; x, z) and V(t;x, z) are horizontal and vertical velocities respectively. Hence, differentiation of (<ref>) yields 

    U(t;x, z)    = ∑_n=1^N [ {A_n cos(ω_n t  - k_n x) + B_n sin(ω_n t  - k_n x)}ω_n cosh[k_n(d+z)]/sinh(k_nd)], 
    
        V(t; x, z)    = ∑_n=1^N [ {-A_n sin(ω_n t - k_n x) + B_n cos(ω_n t  - k_n x)}ω_n sinh[k_n(d+z)]/sinh(k_nd)], 
    U̇(t;x,  z)    = ∑_n=1^N [ {-A_n sin(ω_n t - k_n x) + B_n cos(ω_n t  - k_n x)}ω_n^2 cosh[k_n(d+z)]/sinh(k_nd)], 
    V̇(t; x,  z)    = ∑_n=1^N [ {-A_n cos(ω_n t - k_n x) - B_n sin(ω_n t - k_n x )}ω_n^2cosh[k_n(d+z)]/sinh(k_nd)],

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



 §.§ JONSWAP wave spectrum
 

In Section 3.3.2, we assume the wave spectrum S(ω; 𝐗) has a JONSWAP parametric form, which is shown by <cit.> to be suitably flexible to capture measure offshore spectral behaviour. The JONSWAP spectrum has spectral density function 

    S(ω ; 𝐗=𝐱)=αω^-rexp{-r/4(|ω|/ω_p(𝐱))^-4}γδ(ω ; 𝐱),

for ω>0, where 𝐗=(H_s, S_2) and ω_p(𝐱) = 2π/t(𝐱), where t(𝐱) is the observed value of the second spectral moment wave period T_2 = √((2π H_S)/(g S_2)) in sea state 𝐗=𝐱, with 

    δ(ω ; 𝐱)=exp{-1/2(0.07+0.02 · I{ω_p(𝐱)>|ω|})^2(|ω|/ω_p(𝐱)-1)^2},

and constants r, α, γ > 0. The normalising constant α is chosen so that 

    4·{∫_-∞^∞ S(ω; 𝐱)ω}^1/2 = h(𝐱),

where h(𝐱) is the observed value of significant wave height H_S in sea state 𝐗=𝐱.



 §.§ Efficient wave simulation using the fast Fourier transform


We apply the fast Fourier transform (FFT) algorithm <cit.> 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'_n = A_n + 𝒬σ^2_n   and   B'_n = B_n + ℛσ^2_n ω_n,

we write

    E(t_j)    = ∑^N_n=1 g^(η)_n exp{-2π i/Njn},    g^(η)_n = A'_n + B'_n i, 
    
        U(t_j, z)    = ∑^N_n=1 g^(u)_n exp{-2π i/Njn},     g^(u)_n = (A'_n + B'_n i) ·ω_n  ·cosh{k_n(z+d)}/sinh(k_nd),
    U̇(t_j, z)    = ∑^N_n=1 g^(u̇)_n exp{-2π i/Njn},    g^(u̇)_n = (B'_n-A'_ni) ·ω_n^2 ·cosh{k_n(z+d)}/sinh(k_nd),

where i^2 = -1 and j=1, …, N. This allows us to utilise the fast Fourier transform <cit.> to efficiently compute the linear wave behaviour given by expressions (10) and (11) for z ∈ 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' = 
        z    for   z≤ 0, 
    
        0    for   z>0,

rather than z.



§ STORM PEAK HIERARCHICAL MODEL SELECTION
 


 §.§ 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_GEV(x)=exp{-[1+ξ(x-μ_G/σ_G)]^-1 / ξ},

defined on the set {x: 1+ξ(x-μ) / σ>0}, with location μ_G ∈, scale σ_G >0 and shape ξ∈. Second, the Weibull distribution, which has distribution function

    F_Weibull(x)=1-e^-(x / λ)^k,

for x ≥ 0, with scale λ>0 and shape k>0. Third, the Lognormal distribution, with density function

    f_Lognormal(x) = 1/x σ_L √(2 π)exp(-(ln x-μ_L)^2/2 σ_L^2),

for x >0, with location μ_L and log-scale σ_L>0. Fourth, the gamma distribution, with density function

    f_Gamma(x) = β^α/Γ(α) x^α-1 e^-β x,

where Γ is the gamma function, for x ≥ 0, with shape α>0 and rate β>0.  



 §.§ Detailed cross validation results

Here we present the full details of the results summarised by Table 1 in Section 5.2, summarised in Table <ref> 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_vCVK' contain K-fold cross validation scores for an extreme threshold v= 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. 






