Statistical Modeling and Applications of Joint Distributions for Significant Wave Height, Spectral Peak Period, and Peak Direction of Propagation: A Case Study in the Norwegian Sea

Abstract

The estimation of long-term extreme response is a crucial task in the design of marine structures. The target extreme responses are typically defined by annual exceedance probabilities of 10⁻² and 10⁻⁴. Various approaches can be employed for this purpose, with preference given to statistical long-term analysis, which involves aggregating the exceedance probabilities of all potential sea states contributing to the exceedance of the target extremes. A joint model encompassing important metocean parameters such as wind, waves, and current is often necessary. This study specifically focuses on waves and wave-induced responses. In characterizing short-term sea state conditions, significant wave height ($H_s$), spectral peak period ($T_p$) and peak direction of propagation (${\mathsf{\Phi }}_p$) are identified as the most important sea state characteristics. The objective of this work is to present the results of the joint model for the three sea state parameters, i.e., $H_s$, $T_p$ and ${\mathsf{\Phi }}_p$, at an offshore site in the Norwegian Sea. The conditional modeling approach is applied using long-term hindcast data, and different statistical models are discussed for fitting the marginal and conditional distributions. The fitted parameters for all directional sectors are provided, offering a comprehensive representation of the joint model for direct use in long-term response analysis. Two case studies are included to illustrate the application of the fitted joint model in long-term response analyses. The case studies identify the governing wave directions and the most important combinations of short-term sea state characteristics regarding the estimation of long-term extreme responses.

1. Introduction

The design and operation of marine structures necessitate a comprehensive understanding of the long-term variations in essential metocean processes to mitigate overloading and fatigue effects. This research study specifically focuses on structures that endure continuous operation at specific locations for the operational life, typically ranging from 25 to 50 years. The predominant metocean processes of major interest include waves, wind, and current. Examples of such structures encompass fixed and floating installations utilized in the oil and gas industry, fixed and floating wind turbines, and floating structures employed in the fish-farming industry; the latter are anticipated to be increasingly deployed in more exposed offshore environments in the future. It should be noted that the desired safety level may vary across different applications. To verify that a structural design fulfils the required target safety level against overload failures, Ultimate Limit State (ULS) and Accidental Limit State (ALS) action effects must satisfy the limit state equation below [1]:

$${\gamma }_p{x}_p+{\gamma }_v{x}_v+{\gamma }_e{x}_e \le \frac{y_c}{{\gamma }_m}$$

(1)

where various action effects are denoted by x and the indices p, v and e refer to permanent action effects (caused by weight of structure and permanent equipment), variable operational actions effects (caused by the loading capacity specified by the vessel owner) and the external action effects (caused by wind, waves and current, etc.), respectively. $y_c$ is the limiting capacity for the action effect under consideration. Involved partial safety factors are denoted by $\gamma$ with indices as given above and $m$ representing the partial safety factor for the limiting capacity. National rules and regulations define how the various quantities of Equation (1) shall be determined.

In this study, the focus is on estimating the target characteristic values of the externally generated action effect $x_e$. For ULS control, $x_e$ is frequently defined as the value corresponding to an annual exceedance probability of 10⁻². This value is multiplied by a safety factor, ${\gamma }_e$, typically 1.3–1.4. For ALS control, $x_e$ is often defined as the value exceeded by a probability of 10⁻⁴ per year. In most cases, the safety factors are set to 1 in connection with the ALS control.

To estimate the target response extremes for ULS and ALS, it is necessary to add up the exceedance probabilities of all possible sea states that can be faced at the site or, rather, all metocean conditions that contribute significantly to the exceedance of target extremes. At the Norwegian Continental Shelf (NCS), the most common approach has been to use the all-sea-state approach [2,3]. In this approach, the long-term probability of exceeding a response level is calculated as a weighted sum of the short-term probabilities of exceeding the response level, where the weights are the probability of occurrence of the various sea states. This has been the default approach at NCS since the late seventies/early eighties. Here, the short-term sea state refers to a stationary or weakly stationary random wave process for a given length of time, where the characteristics of the random wave such as total variance, distribution of variance versus wave number and frequency, and direction of propagation are close to constant. For Norwegian waters, it is commonly assumed that stationarity holds for a duration of 3 h. A short-term sea state typically exhibits a combined nature, where a wind-generated sea coexists with an incoming swell system generated by a distant storm event or a dying wind sea due to an abrupt change in wind conditions [4]. Accounting for the combined systems can be important in the planning and designing of marine operations, as wind sea and swell sea can have comparable severity in terms of significant wave height for mild sea states. When considering extreme responses in relation to ULS and ALS design, the influence of swell is generally limited, as highlighted in studies such as [5,6]. Therefore, this study will treat the sea surface as a single system characterized by an analytical wave spectrum model.

In the statistical long-term response analysis, a joint model of important sea state parameters is required. For short-term sea state conditions, significant wave height ($H_s$), spectral peak period ($T_p$) and peak direction of propagation (${\mathsf{\Phi }}_p$) are identified as the most important sea state characteristics. Several methods, such as maximum likelihood model [7], the conditional modeling approach [8] and the Nataf model [9] can be used to derive a joint model for these characteristics. Discussions of various methods can be found in Bitner-Gregersen et al. [10]. Among these methods, the conditional modeling approach is straightforward and has been widely used for various design purposes. However, when the available wave data are limited and the characteristics included in the joint model are too many, the joint modeling will lead to large uncertainties [11]. In such cases, a simplified method involves utilizing an omni-directional wave climate for all direction sectors. A joint distribution of $H_s$ and $T_p$ from all directions is firstly established and then applied in the response analysis for all directions. This method has often been used in situations when only a limited amount of high-quality wave data are available. While this approach is rather conservative for the mildest sectors, the adequacy of the approach may well be acceptable for the most severe directional sectors. In the literature, one can find previously published joint models of $H_s$ and $T_p$ at various locations, e.g., Li et al. [12] and Johannessen et al. [13], which can be used in the simplified method for various design purposes.

To predict extreme responses for structures sensitive to wave directions, it is necessary to incorporate the wave direction of propagation (${\mathsf{\Phi }}_p$) into the joint model. Given the fact that presently good quality hindcast data spanning over 60 years are available at NCS, it is feasible to create a more realistic conditional model that considers the direction of propagation (${\mathsf{\Phi }}_p$). However, the number of published joint models that incorporate the direction characteristic is limited. Hence, the motivation behind this work is to present the comprehensive results of a joint model encompassing three sea state characteristics, $H_s$, $T_p$ and ${\mathsf{\Phi }}_p$, at a site in the Norwegian Sea. This study provides fitted joint models of $H_s$ and $T_p$ for all directional sectors, offering a comprehensive representation of the joint model that can be directly utilized in long-term response analysis for various design purposes. To illustrate the practical application of the joint probability model, two case studies of the long-term response analyses are conducted. The case studies focus on evaluating the critical directions and dominating sea states contributing to the extreme responses.

2. Statistical Assessment of Action and Action Effects

2.1. Description of Short-Term Sea States for Response Estimation

Assuming that the weakly stationary random field, $\mathsf{\Xi }\left(x,y,t\right)$ is Gaussian, the wave spectrum is, in a statistical sense, a complete description of the random process. The 3-dimensional wave spectrum, $S_{\mathsf{\Xi }\mathsf{\Xi }}$, is often written in the form [14]

$${\mathrm{s}}_{\mathsf{\Xi }\mathsf{\Xi }}\left(\omega ,\phi +\psi ;h,t,\phi \right)=s_{\mathsf{\Xi }\mathsf{\Xi }}\left(\omega ;h,t,\phi \right)d_{\mathsf{\Xi }}\left(\psi \right)$$

(2)

It is tacitly assumed that an analytical model parameterized in terms of significant wave height, $H_s$, and spectral peak period, $T_p$, can be used for the frequency spectrum, e.g., JONSWAP spectrum [15], where the frequency is denoted by $\omega$. $d_{\mathsf{\Xi }}\left(\psi \right)$ is the wave spreading function representing the variability of the propagation direction, $\psi$, of the various frequency components around the mean direction of propagation, $\phi$, of the wave system. A spreading function often used is given below [16]:

$$d_{\mathsf{\Xi }}\left(\psi \right)=\frac{\mathsf{\Gamma }\left[\frac{n}{2}+1\right]}{\sqrt{\pi } \mathsf{\Gamma }\left[\frac{n}{2}+\frac{1}{2}\right]}co{s}^n\left\{\psi \right\}; \mathrm{where} \psi \in \left[-\frac{\pi }{2}, \frac{\pi }{2}\right]$$

(3)

$\mathsf{\Gamma }\left( \right)$ is the Gamma function and the factor in front of the cosine term ensures that the integral over the spreading function is 1. In general, the spreading function is also a function of $\omega$. This is neglected here. It is also assumed that the spreading function is the same for all sea states. The latter assumption is rather good for the severe sea states, but it may represent a rather crude approximation for lower sea states where the underlying combined nature of the sea is often more clearly reflected. For the response problem discussed here (long-term extremes) this is not expected to be important. The parameter, n, is defining the width of the spreading function. Low values suggest large spreading, whereas very large values would suggest long-crested sea. For severe sea states, n around 10 may be a good approximation.

For a linear response problem, the response can be performed in the frequency domain. The response–wave relation is given by the transfer function, $h_{\mathsf{\Xi }R}\left(\omega \right)$, where $\mathsf{\Xi }$ denotes the stationary wave process, $\mathsf{\Xi }\left(t\right)$, at a reference position (for example, the projection of the center of gravity of the structure under consideration) and R denotes a selected response process, $R\left(t\right)$, of the structure. The transfer function is a complex-valued function including both amplitude scaling and phase shift between the response component and the underlying wave component. For a linear response quantity, the response process will also be Gaussian and completely characterized by the response spectrum given by

$$s_{RR}\left(\omega , \phi +\psi ;h, t, \phi \right)={\left|h_{\mathsf{\Xi }R}\left(\omega , \phi +\psi \right)\right|}^2{s}_{\mathsf{\Xi }\mathsf{\Xi }}\left(\omega ;h, t,\phi \right)d_{\mathsf{\Xi }}\left(\psi \right)$$

(4)

In order to determine the frequency spectrum, we can integrate over $\psi$:

$$s_{RR}\left(\omega ;h, t, \phi \right)=s_{\mathsf{\Xi }\mathsf{\Xi }}\left(\omega ;h, t,\phi \right){{\displaystyle \int }}_{-\frac{\pi }{2}}^{\frac{\pi }{2}}{\left|h_{\mathsf{\Xi }R}\left(\omega ,\phi +\psi \right)\right|}^2{d}_{\mathsf{\Xi }}\left(\psi \right)d\psi$$

The $k$th spectral moment for the response spectrum is given by

$$m_R^{\left(k\right)}\left(h, t, \phi \right)={{\displaystyle \int }}_0^{\infty }{\omega }^k{s}_{RR}\left(\omega ;h, t, \phi \right)d\omega$$

(5)

From the spectral moments, we can determine the variance of the response process, ${\sigma }_R^2=m_R^{\left(0\right)},$ expected zero-up-crossing frequency, ${\omega }_{z,R}=\sqrt{m_R^{\left(0\right)}/m_R^{\left(2\right)}}$, expected number of-zero-up-crossings in a sea state of duration $\Delta T$ (s), and $n_{z,R}=\frac{{\omega }_{z,R}}{2\pi }\Delta T$. These quantities are needed for estimating short-term response extremes. In the following, short-erm duration is considered as 3 h (10,800 s) and short-term extremes of a response quantity C are denoted as $C_{3h}$.

2.2. Long-Term Response Analysis

At the Norwegian Continental Shelf (NCS), the most common approach has been to use the all-sea-state approach. In this approach, the long-term probability of exceeding a response level is calculated as a weighted sum of the short-term probabilities of exceeding the response level, where the weights are the probability of occurrence of the various sea states. Given the sea-state characteristics, significant wave height $H_s$, peak period $T_p$ and peak direction of propagation ${\mathsf{\Phi }}_p$, the long-term distribution for the largest response in a short-term sea state (3 h), denoted as. $C_{3h}$, is given by [2]

$$F_{C_{3h}}\left(c\right)={\displaystyle {\int }_h}{\displaystyle {\int }_t}{\displaystyle {\int }_{\phi }}{F}_{C_{3h}|H_s{T}_p{\mathsf{\Phi }}_p}(c|h, t,\phi ) f_{H_s{T}_p{\mathsf{\Phi }}_p}\left(h,t,\phi \right)dhdtd\phi$$

(6)

where $F_{C_{3h}|H_s{T}_p{\mathsf{\Phi }}_p}(c|h, t,\phi )$ is the short-term distribution function for the 3 h maximum response for a given sea state and $f_{H_s{T}_p{\mathsf{\Phi }}_p}\left(h,t,\phi \right)$ is the joint probability density function for the selected sea-state characteristics.

If the sea surface is a stationary Gaussian process and the target response can be considered as linear function regarding waves. Then, the global maxima (largest maximum between adjacent zero-down crossings) of the responses follow a Rayleigh distribution and the distribution of $C_{3\mathrm{h}}$ is obtained by raising the Rayleigh distribution to the expected number of global maxima within 3 h. For nonlinear response problems, the common approach is to perform 3 h time-domain simulations for each sea state using different realizations of the wave process. The number of wave realizations for each sea state must be large enough to reflect the probabilistic structure of the 3 h extreme value. The maximum response from each simulation is identified and a proper distribution function is fitted to the sample of 3 h extremes. The simulations must be performed for a large number of different sea states in order to cover the sample space of $H_s$ and $T_p$. Very often, the 3 h extreme value is assumed to follow a Gumbel distribution:

$$F_{C_{3\mathrm{h}}|H_s{T}_p{\mathsf{\Phi }}_p}(c|h, t,\phi )=\exp \left\{-exp\left(-\frac{c-\lambda \left(h,t\right)}{\kappa \left(h,t\right)}\right)\right\}$$

(7)

where $\lambda \left(h,t\right)$ and $\kappa \left(h,t\right)$ are estimated Gumbel parameters. A continuous function of these parameters as functions of $H_s$ and $T_p$, often referred to as response surfaces, can be established based on the estimated parameters for many different combinations of sea-state characteristics. An example of the obtaining response surfaces for Gumbel parameters is shown in Li et al. [17].

Regarding the joint probability density function of the sea-state characteristics $f_{H_s{T}_p{\mathsf{\Phi }}_p}\left(h,t,\phi \right)$, the estimation of this probability function will be detailed in Section 3. Here, it is mentioned that the peak direction of propagation, ${\mathsf{\Phi }}_p$, is modeled as a discrete variable in this study. It is described by the probability mass function for 12 sectors of width 30°, ${\phi }_1,{\phi }_2, \dots \dots {\phi }_{12}$ with corresponding probabilities $p_{{\phi }_1},p_{{\phi }_2}, \dots , p_{{\phi }_{12}}$. Accordingly, Equation (6) is rewritten to the form

$$\begin{array}{c}\hfill F_{C_{3h}}\left(c\right)={\displaystyle \sum _{i=1}^{12}}{F}_{C_{3h}|{\mathsf{\Phi }}_p}(c|{\phi }_i) p_{{\phi }_i}\\ \hfill ={\displaystyle {\sum }_{i=1}^{12}} \left\{{\displaystyle {\int }_h}{\displaystyle {\int }_t}{F}_{C_{3h}|H_s{T}_p{\mathsf{\Phi }}_p}(c|h,t,{\phi }_i)f_{H_s{T}_p|{\mathsf{\Phi }}_p}\left(h,t|{\phi }_i\right)dhdt\right\} p_{{\phi }_i}\end{array}$$

(8)

From the conditional long-term distribution given the sector ${\phi }_i,$ one can estimate the conditional response for this sector corresponding to a return period of $M$ (years) by solving

$$1-F_{C_{3\mathrm{h}}|{\mathsf{\Phi }}_p}\left(c_{\frac{1}{K_i}}|{\phi }_i\right)=\frac{1}{M·K_{1Y}·p_{{\phi }_i}}=\frac{1}{K_i}$$

(9)

$K_{1Y}=2920.25$ is the number of short-term events of 3 h duration per year, and $K_i$ is the expected number of occurrences in sector i during M years. The usefulness of conditional response extremes can be discussed if they are much lower than the marginal extremes determined below. It is usually other sectors that are important regarding extremes if the target sector is much milder than the governing sectors.

The full long-term analysis is conducted by first performing a conditional long-term analysis for all sectors. Thereafter, the marginal long-term distribution is found as a weighted sum of the conditional long-term distributions with sector probabilities as weights. The response value corresponding to a return period of M years or, equivalently, an annual exceedance probability of $1/M$, $c_{1/M}$, is estimated by solving

$$1-F_{C_{3\mathrm{h}}}\left(c_{\frac{1}{M}}\right)=\frac{1}{M·K_{1Y}}$$

(10)

In an ideal case with an infinite amount of data, $c_{\frac{1}{K_i}}\le c_{\frac{1}{M}}, i=1, 2, \dots ., 12$. A condition for establishing the long-term response distribution using the all-sea-states approach is the availability of a good quality joint density function for the wave characteristics $H_s, T_p$, and ${\mathsf{\Phi }}_p$.

3. Joint Probabilistic Modeling of Sea State Characteristics

3.1. Database and Reference Site

The wave data used in the present study are obtained from the database Norwegian Reanalysis 10 km (NORA10) developed by the Norwegian Meteorological Institute, covering over 60 years of wind and wave characteristics for every 3 h [18,19]. The database was validated by comparing hindcast data with available measurements and it was concluded that the hindcast data were adequate for design purposes, see, e.g., Bruserud and Haver [20]. The dataset includes wind speed and direction, $H_s$, $T_p$ and ${\mathsf{\Phi }}_{\mathrm{p}}$ for the total sea, wind sea and swell sea. In this study, we consider the data of the total sea. The NORA10 data of one reference site located in the Norwegian Sea with coordinates 67.05° N, 7° E (see Figure 1) were chosen for the present study. The hindcast wave data cover a period of 61 years, from September 1957 to October 2018, and the total number of 3 h sea states in the dataset is 178,725. The omni-directional scatter diagram is given in Figure A1 in Appendix A.

To model the effect of the direction of propagation, the data are divided in 12 directional sectors with an equal width of 30°. Information about the 12 sectors is given in

Table 1. Main characteristics of the 12 directional sectors for the reference site and the probability mass function for the sectors, $p_{{\phi }_i}$.

$\mathbf{Sector} {\mathit{\phi }}_{\mathit{i}}$	Sector’s Range	Sector’s Midpoint	${\mathit{p}}_{{\mathit{\phi }}_{\mathit{i}}}$
1	[345° 15°]	0°	0.1094
2	[15° 45°]	30°	0.1790
3	[45° 75°]	60°	0.0263
4	[75° 105°]	90°	0.0046
5	[105° 135°]	120°	0.0069
6	[135° 165°]	150°	0.0267
7	[165° 195°]	180°	0.0162
8	[195° 225°]	210°	0.1094
9	[225° 255°]	240°	0.2904
10	[255° 285°]	270°	0.1369
11	[285° 315°]	300°	0.0453
12	[315° 345°]	330°	0.0490

. Most of the waves comes from south-southwest to west of the reference site, i.e., sectors [195° 225°], [225° 255°] and [255° 285°]. The west-northwest sector, [285° 315°], has just 5% of occurrences, but as will be seen later, sea conditions from this sector are severe. It must therefore be included among the worst sectors. North and north-northeast sectors [345° 15°] and [15° 45°] are also associated with high percentages of wave data, but as will be shown later, of less severity than the southwest-northwest sectors. More than 85% of the data are coming from the 6 sectors listed above. To illustrate the methodology, a detailed description for determining the joint distribution of $H_s$ and $T_p$ for the sector [225° 255°] will be given, but the fitted parameters of the joint distribution for all sectors will be presented. The adequacy of the fitted models will be shown by comparing model characteristics with the empirical distribution based on the data.

3.2. Conditional Distribution of $H_s$ Given ${\Phi }_p$

In this study, the approach to obtain the joint description of the metocean parameters is based on the Conditional Modeling Approach (CMA). For a given sector ${\mathsf{\Phi }}_p=\phi$, the expression of the conditional probability density function (PDF) of the joint model of $H_s$ and $T_p$ for a given sector, $f_{H_s{T}_p|{\mathsf{\Phi }}_p}\left(h,t|\phi \right)$, is given in Equation (11).

$$f_{H_s{T}_p|{\mathsf{\Phi }}_p}\left(h,t|\phi \right)=f_{H_s|{\mathsf{\Phi }}_p}\left(h|\phi \right)·f_{T_p|H_s,{\mathsf{\Phi }}_p}\left(t|h,\phi \right)$$

(11)

where $f_{H_s|{\mathsf{\Phi }}_p}\left(h|\phi \right)$ is the conditional distribution of $H_s$ given ${\mathsf{\Phi }}_p$, and $f_{T_p|H_s,{\mathsf{\Phi }}_p}\left(t|h,\phi \right)$ is the conditional distribution of $T_p$ given $H_s$ and ${\mathsf{\Phi }}_p$. The current method is commonly used in the literature and design codes, examples can refer to Li et al. [12] and Johannessen et al. [13].

Several distributions have been considered to model the marginal distribution of $H_s$ for the 12 sectors. The most suitable model for the data from the reference site is found to be the 3-parameter Weibull, which is the distribution function that is most frequently used as the long-term distribution of $H_s$ when using the all-sea-states approach. The cumulative distribution function (CDF) of the 3-parameter Weibull model is shown in Equation (12).

$$F_{H_s|{\mathsf{\Phi }}_p}(h|\phi )=1-\exp \left[-{\left(\frac{h-c_{3w}}{{\alpha }_{3w}}\right)}^{{\beta }_{3w}}\right]; h>c_{3w}$$

(12)

where ${\alpha }_{3w}$ is the scale parameter, ${\beta }_{3w}$ is the shape parameter and $c_{3w}$ is the location parameter. The parameters of the directional distributions are estimated based on the method of moments. The fitted results for sector 9 corresponding to [225°–255°] are shown in Figure 2, where the data and the fitted distribution are plotted in the Gumbel probability paper.

As indicated in Figure 2, the fitted distribution shows good agreement with the empirical distribution for $H_s$ within the range of 1 m to 11 m. When it comes to estimating extreme responses, a possible deviation for less than 1 m is not a significant concern. The primary focus lies on the upper tail of the distribution. There is a slight curvature in the upper tail as observed in the empirical distribution. The exact cause of this curvature is difficult to identify; it can be an underlying trend, a consequence of clustering in correlated data (in the all-sea-states approach, adjacent data are highly correlated) or, simply, a consequence of limited data availability. It is noted that the fitted distribution remains rather accurate for estimating the conditional 10⁻² annual probability value. Regarding the estimation of the 10⁻⁴ annual probability value, it will be more uncertain. It is, however, expected that the values estimated from the fitted 3-parameter Weibull will be on the safe side. This is good for the ALS design since the partial safety factor is usually 1 for the ALS limit state. What is presented here is representative for most sectors, but there are exceptions for the coastal sectors (Sectors 2–4). This will be further discussed in Section 3.4. The fitted parameters for the 3-parameter Weibull distribution for all sectors are given in

Table 2. Estimated parameters of 3-parameter Weibull distribution for marginal distribution of $H_s$ for each sector (see Equation (12)), and corresponding estimated $H_s$ extremes. The last row provides the extremes by combining all sectors (omni-directional results).

$\mathbf{Sector} {\mathit{\phi }}_{\mathit{i}}$	${\mathit{a}}_{3\mathit{w}}$	${\mathit{b}}_{3\mathit{w}}$	${\mathit{c}}_{3\mathit{w}}$	${\mathit{H}}_{\mathit{s}}$,1-year (m)	${\mathit{H}}_{\mathit{s}}$,100-year (m)	${\mathit{H}}_{\mathit{s}}$,10,000-year (m)
1	1.99	1.38	0.78	8.29	11.66	14.98
2	1.39	1.19	0.85	7.80	11.26	14.88
3	1.26	1.22	0.77	5.32	8.33	11.38
4	1.44	1.56	0.63	3.57	5.72	7.62
5	1.98	1.56	0.84	5.23	8.11	10.69
6	2.31	1.65	0.68	6.68	9.43	11.94
7	2.47	1.47	0.91	7.56	11.40	14.99
8	2.75	1.46	0.68	10.27	14.27	18.14
9	2.36	1.31	0.73	11.37	15.78	20.25
10	2.32	1.25	0.71	11.04	16.07	21.21
11	2.35	1.25	0.61	9.58	14.80	20.08
12	2.37	1.37	0.75	8.91	13.11	17.22
All sectors	-	-	-	12.36	16.75	21.29

.

As the conditional distribution functions, $F_{H_s|{\mathsf{\Phi }}_p}\left(h|{\phi }_i\right)$, are estimated for $i=1, 2, \dots , 12$, the long-term marginal distribution for $H_s$ is obtained using the two first terms Equation (8) when replacing $C_{3h}$ with $H_s$. Directional extremes of $H_s$ corresponding to a return period of $M$ (years), $h_{s,M}^{({\phi }_i)}$, can be estimated by solving Equation (9) after replacing $C_{3h}$ with $H_s$. For the omni-directional case, target extremes can be found by Equation (10) after replacing $C_{3h}$ with $H_s$.

Directional and omni-directional extremes for $M=1, 100, 10000$ are estimated and given in Table 2. The highest extreme significant wave heights are associated with the western sectors [225° 255°], [255° 285°] and [285° 315°]. The sector that accounts for the highest percentage of observations (sector 9 [225° 255°]) results in the second highest extreme significant wave height among all sectors. It is seen that the extreme significant wave height decreases eastwards mainly because of reduced fetch. One can also note that for the studied site, the omni-directional value is slightly larger than the corresponding worst directional value. The effect reduces from about 9% at the 1-year level to about 4% for the 100-year level and to less than 1% for the 10,000-year level.

3.3. Conditional Distribution of $T_p$ Given $H_s$ and ${\Phi }_p$

To estimate the conditional distribution of $T_p$ given $H_s$ and ${\mathsf{\Phi }}_p$, the data for $T_p$ are firstly binned for different $H_s$ classes for each sector. Since the sectors have different amounts of data, the sample size of the conditional distribution differs for the various sectors; see the percentage for all sectors in Table 1. Two probabilistic models were considered for this study: the Lognormal distribution, that is mostly used for the present problem, and the 3-parameter Weibull distribution. Selecting the appropriate model for the classes of $H_s$ with sufficient data is challenging. However, it is essential for the conditional distribution of $T_p$ given $H_s$ and ${\mathsf{\Phi }}_p$ to be applicable for the range of $H_s$ beyond available observations. In fact, this is the most important part of the joint distribution. This means that within the framework of conditional modeling, the distribution parameters must be extrapolated based on findings from $H_s$-classes with sufficient data. The Lognormal distribution proves to be more amenable considering these aspects, and thus in this context. The PDF of the Lognormal model is given as

$$f_{T_p|H_s{\mathsf{\Phi }}_p}\left(t|h,\phi \right)=\frac{1}{\sqrt{2\pi }·{\sigma }_{ln\left(T_p|H_s{\mathsf{\Phi }}_p\right)}·t}\exp \left[-\frac{1}{2}{\left(\frac{\ln \left(t\right)-{\mu }_{ln\left(T_p|H_s{\mathsf{\Phi }}_p\right)}}{{\sigma }_{ln\left(T_p|H_s{\mathsf{\Phi }}_p\right)}}\right)}^2\right]$$

(13)

where the distribution parameters ${\sigma }_{ln\left(T_p|H_s{\mathsf{\Phi }}_p\right)}$ and ${\mu }_{ln\left(T_p|H_s{\mathsf{\Phi }}_p\right)}$ are the standard deviation and mean of $ln\left(T_p\right)$ for given $H_s$ class and sector ${\mathsf{\Phi }}_p$, respectively. These parameters must be estimated for all classes of $H_s$ and sectors of ${\mathsf{\Phi }}_p$ with data. For Sector 9 [225° 255°], the mean and variance of $ln{T}_p$ are given in

Table 3. The mean, $m_{ln{T}_p}$, and variance, ${\sigma }_{ln{T}_p}^2$, of $ln\left(T_p\right)$ for Sector 9 [225° 255°] (see Equation (13)).

${\mathit{H}}_{\mathit{s}}$ (m)	<1	1–2	2–3	3–4	4–5	5–6	6–7	7–8	8–9	9–10	10–11	>11.0
$m_{ln{T}_p}$	2.12	2.25	2.37	2.47	2.53	2.57	2.61	2.64	2.68	2.71	2.73	2.79
${\sigma }_{ln{T}_p}^2$	0.06	0.06	0.06	0.04	0.04	0.03	0.02	0.02	0.01	0.01	0.01	0.01

for the $H_s$-classes with data. There are so few data above 11 m that the data are pooled into one class. The $H_s$ for the merged class, 12.9 m, is the mean $H_s$ value of the data in the class.

When establishing the joint distribution of $H_s$ and $T_p$ for an omni-directional case or for each sector, the conditional distribution of $T_p$ given $H_s$ and ${\mathsf{\Phi }}_p$ must be available for all values of $H_s$ up to the design return period, say 100–10,000 years. This is a challenge regarding the 10,000-year level, which is well above the period that can be observed. The Lognormal parameters for $T_p$ for the upper possible range of $H_s$ must therefore be obtained by extrapolation. The functions given in Equations (14) and (15) are frequently adopted for the conditional mean and conditional variance of $ln{T}_p$, respectively.

$${\mu }_{ln\left(T_p\right)|H_s{\mathsf{\Phi }}_p}=a_1+a_2{h}^{a_3}$$

(14)

$${\sigma }_{ln\left(T_p\right)|H_s{\mathsf{\Phi }}_p}^2=b_1+b_2\exp \left(-b_{3}h\right)$$

(15)

$a_1$, $a_2$, $a_3$, $b_1$, $b_2$ and $b_3$ are the fitted parameters. For Sector 9, these functions are fitted to the characteristics given in Table 3 by a least-squares approach. The fitted curves are compared to the class estimates in Figure 3. The fitted curves follow the class estimates rather well, but there are uncertainties associated with extrapolated ranges. A similar fitting procedure is carried out for all sectors. The obtained coefficients of Equations (14) and (15) are given in

Table 4. Estimated parameters of Lognormal distribution for conditional distribution of $T_p$ given $H_s$ and ${\mathsf{\Phi }}_p$ (see Equations (14) and (15)).

$\mathbf{Sector} {\mathit{\phi }}_{\mathit{i}}$	${\mathit{a}}_{1,\mathit{i}}$	${\mathit{a}}_{2,\mathit{i}}$	${\mathit{a}}_{3,\mathit{i}}$	${\mathit{b}}_{1,\mathit{i}}$	${\mathit{b}}_{2,\mathit{i}}$	${\mathit{b}}_{3,\mathit{i}}$
1	0.95	1.05	0.21	0.001	0.05	−0.25
2	0.57	1.44	0.15	0.001	0.05	−0.22
3	1.41	0.42	0.42	0.001	0.05	−0.33
4	1.52	0.15	0.94	0.001	0.18	−1.43
5	−0.02	1.75	0.11	0.001	0.03	−0.11
6	1.34	0.48	0.39	0.001	0.06	−0.50
7	1.29	0.57	0.37	0.001	0.04	−0.16
8	1.32	0.77	0.24	0.001	0.06	−0.23
9	1.51	0.70	0.24	0.001	0.07	−0.18
10	1.49	0.65	0.26	0.001	0.08	−0.23
11	1.53	0.46	0.40	0.001	0.04	−0.21
12	0.07	1.90	0.14	0.001	0.04	−0.21

for all directional sectors. The adequacy of the fitted curves for selected sectors is shown in Figure A3 in Appendix A. For most sectors, a reasonable agreement between the data and the fitted model is seen. For some coastal sectors, Sector 2, Sector 3, and Sector 4, the fitted models of the conditional variance of $ln{T}_p$ are rather uncertain due to fewer data. However, the sea conditions for these sectors are much milder than for the other sectors. Thus, in most cases, these uncertainties have little effect regarding estimating design loads.

The conditional distribution of $T_p$ given $H_s$ and ${\mathsf{\Phi }}_p$ for two $H_s$ classes is given in Figure 4. In this study, the Lognormal model was selected due to the adequacy of the simple smoothed functions for the distribution parameters. It is seen from Figure 4 that the fitted distribution is reasonably accurate for a central 90% band. As conditional cumulative probability approaches around 0.98 (equivalent to a value of 4 on the Gumbel scale), there are very limited data available, leading to increased uncertainty in the upper tail of the distribution. Nevertheless, it is not expected that the extreme upper tail of the conditional distribution will be significant for most response cases. This will be illustrated later in the case study with response quantity. The conditional distribution of $T_p$ is obtained by using the smoothed functions for distribution parameters, as presented in Table 4.

3.4. Discussion of Joint Model for Response Analyses

The complete joint model for the reference site has been established in this paper, and all the fitted parameters are presented in Table 1, Table 2 and Table 4. Since there exists strong correlation among closely spaced (in time) data, a classical hypothesis test is not applicable in this context. Instead, the goodness of fit for the fitted models is assessed by comparing them with the corresponding empirical distributions. Results for Sector 9 [225° 255°] are presented in Figure 2, Figure 3 and Figure 4 in the main text and more fitted models can be found in Figure A2, Figure A3 and Figure A4 in the Appendix A.

Regarding the fitted models for $H_s$ for different sectors, a reasonable fit is obtained for most sectors as seen from Figure A2 in the Appendix A. It is observed for some few sectors the fitted model seems to be non-conservative, i.e., Sector 3 and Sector 5. This is not a concern here because none of these sectors are expected to be important regarding design responses. For Sector 11 and Sector 12, the fitted models seem slightly conservative in the upper tail.

To investigate further the omni-directional model for $H_s$, Figure 5 compares the fitted models using two methods. Method A is a weighted sum of the distributions from all sectors, while Method B is a 3-parameter model fitted to the omnidirectional data directly. The empirical distribution for omni-directional data is also shown in Figure 5. Good agreement for both fitted models for almost the full range of data can be observed. For $H_s>14 \mathrm{m}$, the fitted models seem to be rather conservative. However, the number of observations above this level is merely 9 based on around 60 years of data. Thus, the shape of the upper tail is associated with large uncertainties. The shape of the upper tail is crucial for the estimation of ULS- and ALS-extremes, and the model behavior should at least not be unconservative. This study shows that the tail behavior of Method A is a valid model for practical applications of the all-sea-states approach. There are few data points in the upper tail, and it is uncertain if the observed upper tail data will be representative for the long-term behavior. It is worthwhile to note that the model fitted to all data, Method B, agrees very well with the model obtained as a weighted sum of the directional distributions, Method A. Exposing the structure for the Method B model for the worst sectors covering an accumulated width of 90°–120° would most likely yield reasonable, but conservative design action effects for both ULS and ALS. This can be a much faster approach rather than using Method A.

Compared with fitting of models for $H_s$ where all data in each sector are used, the modeling of the conditional distribution of $T_p$ given $H_s$ and ${\mathsf{\Phi }}_p$ is more challenging for long-term analyses of the wave-induced response. First, there are much less data for estimating the parameters of the selected conditional model as sea-state severity increases. Second, there is a limit for how low $T_p$ can be for a given $H_s$. It is, however, not an easy task to establish a robust lower limit for $T_p$, in particular since the steepest sea states often are of a combined nature, i.e., a dominating growing wind sea riding on a more or less long period swell sea [21]. To model a mixed population with a single probabilistic model is challenging. For this study, it is a challenge for low and moderate values of $H_s$.

Moreover, the parameters of the conditional distribution of $T_p$ are estimated for lower values of $H_s$ compared to the extreme values for ULS or ALS analysis. Extrapolation of these parameters for long-term analysis is necessary but challenging. This is one reason for selecting the lognormal distribution as the conditional probabilistic model for $T_p$ in the present study. With this model, the smoothed functions, Equations (14) and (15), for the two parameters were reasonably well defined. Other distributions, e.g., a 3-parameter Weibull, may be equally good, but extrapolation becomes very difficult. The comparison of fitting the conditional distribution of $T_p$ for various $H_s$ classes using Lognormal distribution and Weibull distribution for sector 9 are presented in Figure A4 in Appendix A.

The fitting of parameters associated with the conditional distribution of $T_p$ for selected sectors is shown in Figure A3 in Appendix A. The fitted curve for the conditional expected value of $ln\left(T_p|H_s{\mathsf{\Phi }}_p\right)$ is generally good in the range of data deduced from observations, and the non-linearity of the relation with $H_s$ is rather limited; therefore, the relation is expected to give reasonable values for both ULS and ALS levels of $H_s$. More scatter is observed for the function of the conditional variance of $ln\left(T_p|H_s{\mathsf{\Phi }}_p\right)$. However, the data clearly indicate that the variance is expected to approach zero for very large values of $H_s$. In this study, a lower limit of 0.001 is introduced for the variance, as given in Table 4. It is expected that this will not have any effect regarding long-term response extremes.

4. Case Studies

The purpose of establishing the joint distribution function of $H_s$, $T_p$ and ${\mathsf{\Phi }}_p$ is that it shall be used for obtaining characteristics for loads and responses for ULS-design and ALS-design. Target probabilities of exceeding these characteristics are 0.01 and 0.0001 per year for ULS and ALS, respectively. From the discussions of the joint model, it is expected that the model is accurate or slightly conservative for severe sea states of the important directions of propagation, i.e., sectors openly exposed to weather from the southwest to northwest. For other domains, e.g., sea states propagating from the coast and rather low sea states irrespective of direction of propagation, the model might be less accurate. An intriguing question is whether the model is adequate for ULS and ALS response predictions with sufficient accuracy. To provide further discussions on this issue, two response problems are presented below to illustrate the application of the established joint model. From the analysis, the most important directional sectors, and the dominating sea states for estimating ULS and ALS response extremes, are identified.

4.1. Case Study on the Wave Crest Height

The response quantity for the first case study is the wave crest height, which is a key quantity for the air gap assessment for fixed platforms [17]. The short-term crest height distribution is assumed to follow the Forristall second-order crest height model [22]:

$$F_{C|H_s{T}_1}\left(c|h, t_1\right)=1-exp\left\{-{\left(\frac{c}{{\alpha }_{F}h}\right)}^{{\beta }_F}\right\}$$

(16)

where h is significant wave height, $t_1$ is sea-state average wave period, and ${\alpha }_F$ and ${\beta }_F$ are the scale and shape parameters of the distribution, respectively. For long-crested seas, the following expressions are recommended [16]:

$$\begin{gathered} {\alpha }_F=0.3536+0.2892 s_1+0.1060 Ur \\ {\beta }_F=2-2.1597 s_1+0.0968 U{r}^2 \end{gathered}$$

(17)

where $s_1$ is the mean sea-state steepness, $s_1=\frac{2\pi \ast h}{g\ast T_1^2}$, $Ur$ is the Ursell number, $Ur=\frac{h}{k_1^2{d}^3}$, $k_1$ is the wave number corresponding to $T_1$ and d is water depth. In this study, the wave spectrum for all sea states is a JONSWAP spectrum with peakedness factor ${\gamma }_p=2.8$. This is a crude approximation when $H_s$ and $T_p$ characterize total sea; however, it is acceptable for the most extreme seas. A short-term sea state in Norwegian waters is typically assumed to last for 3 h. The distribution function of the largest maximum in 3 h can be approximated by assuming all global crest heights (largest crest between zero-down-crossings) are identically distributed and statistically independent:

$$F_{C_{3h}|H_s{T}_p}\left(c|h,t\right)={\left[1-exp\left\{-{\left(\frac{c}{{\alpha }_{F}h}\right)}^{{\beta }_F}\right\}\right]}^{m\left(t\right)}$$

(18)

where $m\left(t\right)$ is the expected number of global crest heights in 3 h. Introducing Equation (18) for $F_{C_{3h}|H_s{T}_p{\mathsf{\Phi }}_p}\left(c|h, t,{\phi }_i\right)$, the long-term distribution of $C_{3h}$ is obtained from Equation (8). From this distribution, the crest height can be estimated with the return period, M, from Equation (10). Similarly, crest heights with a return period, M, for the various sectors can be derived by Equation (9). The results for extreme crest heights corresponding to ULS and ALS analyses are shown in

Table 5. The 100-year (ULS) and 10,000-year (ALS) extreme crest heights for various sectors (see Equation (9)). The last row provides the extremes by combining all sectors (omni-directional results) (see Equation (10)).

$\mathbf{Sector} {\mathit{\phi }}_{\mathit{i}}$	${\mathit{K}}_{1\mathit{Y},}\left(\mathit{i}\right)$	100-Year (m)	10,000-Year (m)
1	319.4	13.1	17.4
2	522.7	12.6	17.2
3	76.8	9.57	13.5
4	13.4	6.61	8.89
5	20.1	8.32	11.2
6	78.0	11.0	14.4
7	47.3	12.8	17.4
8	319.4	16.2	21.4
9	848.0	17.6	23.5
10	399.7	17.9	24.5
11	132.3	16.1	22.3
12	143.1	14.4	19.6
All sectors	2920.25	19	25.5

, where the extremes for each directional sector as well as for omni-directional results are compared. It is seen that the extreme crest heights using an omni-directional distribution are around 5% higher than those crest heights for the most severe sector, Sector 10.

To conduct a more in-depth analysis of the sea states that have the highest impact on the extremes, the exceedance probability of the extreme values from each individual sea state can be quantified. Denote the ULS extreme value as $c_{ULS}$; the exceedance probability of $c_{ULS}$ for an arbitrary sea state $q\left(h_k,t_l,{\phi }_i\right)$ reads

$$q\left(h_k,t_l,{\phi }_i\right)=\left(1-F_{C_{3h}|H_s{T}_p}\left(c_{uls}|h_k{t}_l\right)\right)·p\left(h_k,t_l|{\phi }_i\right)·p\left({\phi }_i\right)$$

(19)

Summing the exceedance probabilities over $k$ = 1, 2, …, $k_{max}$, $l$ = 1, 2, …, $l_{max}$ and $i$ = 1, 2, …, 12 results in the target 100-year probability per 3 h, which is known a priori to be $q_{0.01}=\frac{0.01}{2920.25}$. Thus, the relative contribution (%) for an arbitrary 3 h sea state is defined as

$$r_{cont}=\frac{q\left(h_k,t_l,{\phi }_i\right)}{q_{0.01}}·100\%$$

(20)

By summing up the contributions from all sea states for each sector, the relative contributions for all sectors are obtained. Figure 6 presents the relative contribution from different sectors for both ULS and ALS scenarios. It is seen that Sector 10 [255° 285°] gives the highest contribution to the extreme responses, and the primary sectors contributing to the responses are concentrated within a 90° sector ranging from 210° to 300°. This suggests that it may not be necessary to utilize the distribution model for all sectors. Instead, focusing on a detailed modeling of this main 90° sector may be sufficient to obtain a realistic design response for a marine structure. Furthermore, the results in Table 5 show that the extremes based on the omni-directional distribution are slightly larger than the characteristics for the worst sector based on the distribution for the given sector. If the proposed 90$°$ sector is used instead, the extreme values will be very close to those based on the omni-directional distribution.

For a more comprehensive investigation of the most critical sector, Equations (19) and (20) can be applied to Sector 10, i.e., removing $p\left({\phi }_i\right)$ and using $q_{0.01}$ (ULS) and $q_{0.0001}$ (ALS) for Sector 10, respectively. The sea states which are of most importance regarding exceeding the ULS characteristic and ALS characteristic are identified, respectively. The results are presented in Figure 7, where only sea states giving a relative contribution larger than 0.1% are included.

It is seen that a relatively narrow range of spectral peak period of 1–2 s, associated with high and steep sea states, plays a significant role for the extreme values. These important domains are close to the steep side of the ULS and ALS metocean contours, respectively. For sufficiently steep sea states, effects of order higher than the second order may affect the extremes. However, as the steepness increases further, wave breaking tends to decrease the crest extremes to a lower level than the second-order crest height.

4.2. Case Study on the Airgap of a Semi-Submersible

For the second case study, the response quantity is the relative surface elevation at a given point on a semi-submersible, which is assumed to operate at the reference site in the Norwegian Sea (see Figure 1). The relative surface elevation is the surface elevation relative to the motions at the target position onboard the platform (wave surface seen from the target position onboard). This is the key quantity to assess the airgap problem for semi-submersibles. The major difference between the previous case study of the wave crest height and this case study is that the relative wave elevation depends not only on the surface waves. It is also significantly affected by the platform motions in waves. The detailed characteristics on the platform and the key derivations of the airgap and relative surface elevation for the target point on the semi-submersible can be found in Haver and Patiño [5].

The chosen point at the semi-submersible has coordinates $x’$ = −47.62 m and $y’$ = −47.62 m in the local body-fixed coordinate system, $\left(x^{\prime }, y^{\prime }\right);$ see Figure 8. The selected point is close to the western corner along the East–West diagonal, i.e., the orientation of the platform is selected to be such that this corner faces the worst wave direction, Sector 10 [255° 285°]. The most critical platform motions include heave of center of gravity (CoG), pitch (rotation about local $y^{\prime }$-axis) and roll (rotation about local $x^{\prime }$-axis). The combined pitch and roll platform motion will lead to a rotation about the South–North diagonal with respect to the earth fixed coordinate system. Figure 9 presents transfer functions for the relative surface elevation for two critical wave directions.

Long-term analyses on the relative surface elevation for this platform have been performed following a similar procedure as the previous case study. First, the short-term distribution of the response for a given sea state is established, following the methods applied in Haver and Patiño [5]. In this reference, a Gaussian assumption is applied for the sea surface elevation process with a linear transfer function of the platform motions from hydrodynamic analysis. Effects of non-linearities in the surface process are accounted for by an asymmetry factor of 1.3 for the points close to the platform projection as recommended by DNV [23]. The wave spectrum is JONSWAP with a peakedness of 2 and a long-crested sea is used. The short-term distribution is thus established given the above assumption and linear motion transfer functions of the platform. Then, using the established joint distribution of the wave characteristics, the long-term distribution of the relative elevation can be obtained following Equation (8), i.e., $C_{3h}$ is now the 3 h maximum relative crest height at the studied target point.

From the long-term distribution, the relative crest height corresponding to an annual probability of exceedance equal to 10⁻² accounting for all directions reads 17.22 m, while the value corresponding to annual probability of 10⁻⁴ reads 22.65 m. The still water airgap for the reference platform is given as 20 m. This means there will be no wave-deck impacts at the ULS level, while a 2–3 m submergence of cellar deck bottom must be expected for the western corner at ALS level.

Following the same procedure as discussed in the case study in Section 4.1, the contributions from different directional sectors and from individual sea states on the extreme responses can be obtained. The relative contribution from the various directions regarding the exceedance of the estimated marginal extremes is shown in Figure 10. It is seen for this case that the results are completely dominated by the western sector, i.e., sector 10 [255° 285°]. This result differs greatly from that observed in the previous case study. This is a combined effect of high crest heights and a large transfer function value of the relative wave elevation for this sector.

The relative contribution from individual sea states to exceed the marginal extremes for the Western sector (Sector 10) is also obtained and presented in Figure 11. It is seen that the most important sea states are along the steep side of contours. The most important sea states are sea states between 10 and 15 s spectral peak periods, which is consistent with the shape of the transfer function for relative surface elevation in the target point, as presented in Figure 9. It can also be found from Figure 9 that the transfer function for relative wave elevation for Sector 10 (central direction of 270°) is much higher than those for Sector 11 (central direction of 300°). This explains why Sector 10 is dominating the extreme responses.

5. Conclusions

The paper presents a joint probabilistic model for the significant wave height, spectral peak period, and spectral peak direction of propagation at an offshore site in the Norwegian Sea. The model is fitted using a dataset spanning 61 years of 3-hourly sea state characteristics provided by the Norwegian hindcast database NORA10. The adequacy of the fitted model is assessed by comparing fitted conditional distributions with empirical distributions estimated directly from corresponding data samples. All the fitted parameters in the joint model encompassing the three sea state parameters are provided, which can be directly utilized in long-term response analysis for various design purposes. Detailed discussions regarding the fitted joint model are made in Section 3.4. The fitted model presented in the paper is valid for an area around the site of observations. The methodology to fit the joint distribution model in this work can be applied to a wide range of locations in the Barents Sea, the Norwegian Sea, and the North Sea. For offshore areas where the extreme response is governed by the occurrence of rare extreme hurricanes, long-term response analyses should be based on a storm-over-threshold method.

To explore the practical implications of the established joint model, it is applied in two response case studies. These case studies aim to identify the most critical directional sectors and sea states contributing to long-term extreme responses during structural design. The findings highlight that high and steep sea states have the most detrimental impact. It is recommended to specifically investigate the sea states that have the highest significance. The major findings from the case studies are summarized below:

• The studies reveal that the characteristic response for design is governed by the three to four critical sectors with a direction width of 30°. Furthermore, the results indicate that sea states with a significant wave height lower than 8 m do not significantly impact the estimated ULS and ALS extreme responses. This suggests that it may not be necessary to develop a distribution model for all sectors. Instead, it is recommended to focus on a joint probabilistic model covering the worst sectors, which can be represented by a width of 90°–120°, and includes all sea states with a significant wave height above 8 m.
• In the second case study, when the structure is oriented with the most unfavorable heading against the worst incoming weather, the directional sensitivity of the extreme responses narrows. However, it is still advisable to focus on a wider design sector spanning 90°–120° around the worst direction. This ensures that the structural design remains appropriate for all directions, accounting for potential variations in weather conditions.

As a final recommendation, based on the observations from the case studies indicating that lower sea states with significant wave height below 8 m have minimal impact on extreme responses, it is advisable to compare the results obtained using the presented all-sea-states approach with an all-storms-over-threshold approach. The all-storms-over-threshold approach estimates the long-term extremes by considering only storms that exceed a specified wave severity threshold; see, e.g., Tromans and Vanderschuren [24] and Stanisic et al. [25]. At present, both methods can be used at NCS, but in recent years, there has been an increasing trend in the application of the all-storms-over-threshold approach due to the availability of very long series of metocean data. However, it is important to note that challenges exist when using the storms-over-threshold approach, including the selection of an appropriate threshold value, and establishing the distribution of extremes over the threshold. Considering these challenges, it is crucial to conduct further case studies and comparisons with the all-sea-states approach in future research. This will contribute to a better understanding of the strengths and limitations of each method and facilitate their improved application and interpretation for the design of marine structures.