The DTU 10-MW reference wind turbine (RWT) is used in this paper, designed from the NREL 5-MW RWT. The wind turbine was designed per the International Electrotechnical Commission (IEC) Class 1A wind regime and is a traditional three-bladed, clockwise rotation-upwind turbine equipping with a variable speed and collective pitch control system. The DTU 10-MW RWT numerical model has been successfully developed and studied in many academic works, e.g., (Hu et al., 2021; Muggiasca et al., 2021; Wang et al., 2022; Yu et al., 2022). The summary of the DTU 10-MW RWT is shown in Table 1.

This work uses a semi-submersible floating structure to support the 10-MW RWT. It was introduced by Dr. techn. Olav Olsen AS in the LIFES 50 + project (Yu et al., 2017). The floater comprises post-tensioned concrete, hosting a central column with three outer columns. The four columns are mounted on a star-shaped pontoon, where a slab is attached at the bottom. Three catenary mooring lines are used to maintain the floater in position, and in each line, a clumped mass is attached, separating the line into two segments. Greater details of the OO-Star Wind Floater and the mooring system are shown in Table 2 and Table 3, respectively.

FAST (Fatigue, Aerodynamics, Structures and Turbulence) (version8, v8.16.00a-bjj), an open-source WT simulation tool developed by the National Renewable Energy Laboratory (NREL), is utilised in this work for the fully coupled aero-hydro-elastic-servo dynamic analysis for the 10-MW FWT. The FAST code couples together five computer codes: AeroDyn (Moriarty and Hansen, 2005), HydroDyn (Jonkman et al., 2014), ServoDyn, and MoorDyn (Hall, 2015), to account for the aerodynamic loads on rotor blades, hydrodynamic loads on floaters, control dynamics, structural dynamics and mooring system dynamics. In addition, FAST provides the interface for reading the time-varying stochastic wind for time-domain simulations. The FAST simulation tool has been successfully used in other well-known projects such as OC3: Offshore Code Comparison Collaboration (Jonkman and Musial, 2010) and OC4: IEA Task Wind 30 (Robertson et al., 2014), and its modelling capability has been authenticated using multiple floating structures in the Netherlands (Coulling et al., 2013).

Consider a long term global response process (such as anchor tension force or FWT surge motion) X (t) of the FWT, measured over a time interval   ( 0 , T ) . Let X 1 , … , X N be measurements of the process X ( t ) at discrete time moments t 1 , … , t N in   ( 0 , T ) . The target is to estimate the distribution function of the extreme value   M N = max  { X j   ; j = 1 , … , N } accurately namely, to estimate CDF (cumulative density function) P ( η ) = Prob ( M N ≤ η ) for large values of the response   η . The following random functions are introduced: A k j ( η ) = 1 { X j > η , X j − 1 ≤ η , … , X j − k + 1 ≤ η }   ,     j = k , … , N ,   k = 2 ,   3 , ... (1) and B k j ( η ) = 1 { X j − 1 ≤ η , … , X j − k + 1 ≤ η }   ,     j = k , … , N ,   k = 2 ,   3 , …   , (2) Where 1 { A } = 1 if A is true, while it is 0 if not. As shown in (Wang, 2001; Bak et al., 2013; Naess and Karpa, 2015a; Naess and Karpa, 2015b; Gaidai et al., 2016; Jian et al., 2018; Gaidai et al., 2019a; Gaidai et al., 2019b; Xu et al., 2019b; Gaidai et al., 2020; Xu et al., 2021) P k ( η ) ≈ exp   ( − ∑ j = k N E [ A k j ( η ) ] E [ B k j ( η ) ] ) ≈ exp   ( − ∑ j = k N E [ A k j ( η ) ] ) ,     η → ∞ (3)

The measured time series can be subdivided into K subsequent (short term) blocks such that E [ A k j ( η ) ] remains approximately constant within each block and such that ∑ j ∈ C i E [ A k j ( η ) ] ≈ ∑ j ∈ C i a k j ( η ) for a sufficient range of η -values, where C i stands for the set of indices for the block with number i ; with   i = 1 , … , K , and where a k j ( η ) are the realised values of A k j ( η ) for the measured time series, then ∑ j = k N E [ A k j ( η ) ] ≈ ∑ j = k N a k j ( η ) . Thus, for a given stationary process (short term sea current state), one has P k ( η ) ≈ exp  ( − ( N − k + 1 ) ε ^ k ( η ) ) (4) Where ε ^ k ( η ) = 1 N − k + 1 ∑ j = k N a k j ( η ) (5)

In the above equations, an assumption of ergodicity has been used for each short term segment of the recorded time series in order to estimate the short term expected values by using observed values of the a k j ( η ) functions. An alternative way of expressing the long term extreme value distribution in Eqn. 4, is obtained by considering the empirical probability distribution of m = 1 , .. , M sea current states having probabilities   p m , so that   ∑ m = 1 M p m = 1 . Next, introduce the long term ACER function of order k: ACE R k ( η ) ≡ ∑ m = 1 M ε ^ k ( η , m ) p m (6) Where ε ^ k ( η , m ) is the same function as in Eqn. 5 but restricted to a specific sea state with number m . As shown in (Gaidai et al., 2020)- (Xu et al., 2019b), the long-term extreme value distribution of   M ( T ) , can then be expressed as follows based on the ACER function of order k: P ( η ) ≈ exp ( − N . ACE R k ( η ) ) (7) Where A C E R k ( η )   is the long term empirical ACER function of order k, with k " N ; N is the total number of data points from the recorded time series used to estimate the ACER functions. Typically, these could be local peaks from the measured time series.

As the order k increases, the accuracy of Eqn. 7 improves and A C E R k ( η ) functions converge conveniently fast with growing k, see (Gaidai et al., 2020)- (Xu et al., 2019b). Note that by increasing the conditioning level k, possible data clustering effects, e.g. with narrow-band response components in the recorded time series, can be accounted for. This is an essential advantage of the ACER method as it increases the accuracy of its extreme predictions and avoids resulting in over-conservative design values.

The A C E R k as functions of the level, η are generally quite regular in the tail, i.e., for high response values of   η . More specifically, for η ≥ η 0 , the tail behaves very closely like e x p { − ( a η + b ) c + d } with a ,   b ,   c ,   d being suitable constants.

It is suggested to do the optimisation on the log-level by minimising the following mean square error function F with respect to the four arguments   a k ,   b k ,   c k , p k , q k . F ( a k ,   b k ,   c k , p k , q k ) = ∫ η 0 η 1 ω ( η ) { ln ( ACE R k ( η ) ) − d k + ( a k η + b k ) c k } 2 d η ,     η ≥ η 0 (8) Where   η 1 is a suitable data cut-off value, i.e., the largest η response value with the confidence with acceptable interval width. Weight function ω is defined as ω ( η ) = { l n C + ( η ) − l n C − ( η ) } − 2 with ( C − ( η ) ,   C + ( η ) ) being a 95% CI, empirically estimated from the measured data. Detailed procedure for further optimising parameters   a k ,   b k ,   c k ,   d k   has been outlined in (Gaidai et al., 2020)- (Xu et al., 2019b).

As already mentioned, the method proposed in this paper will be based on the ACER methodology. This involves both the univariate ACER functions and the bivariate ACER functions. The unique feature of the ACER functions is that they provide the possibility to portray the exact extreme value distribution inherent in the data time series, both the univariate and the bivariate (Wang, 2001)- (Xu et al., 2019b). Hence, the ACER method is fundamentally different from the traditional approach relying on the fitting of hardly asymptotic data to asymptotic extreme value distributions, which are based on the assumption of stationary time series instead of the ACER method. The empirical ACER functions are represented in nonparametric numerical functions with uncertainty bounds. The accuracy obtained depends, of course, on the amount of data available to estimate these functions. It is also an essential feature of the ACER method that it is not limited to stationary time series. It is entirely valid for nonstationary time series as long as the measured data reflects this nonstationarity. For the sake of easy reference, the univariate ACER methodology has been briefly outlined below. The bivariate case follows quite closely the univariate case.

This section presents a statistical bivariate integral correction based on the bivariate ACER method coupled with the Gumbel logistic model (Gaidai et al., 2019a)- (Xu et al., 2019b). Note that this correction is not limited to only extreme value estimates. However, it can be applied with appropriate bivariate models for any statistical values of interest to improve their accuracy based on synchronously measured longer, highly correlated data records.

Let X = max { R 2 ( t ) ; t ∈ [ 0 ,   T return ] } ,   Y = max { R 3 ( t ) ; t ∈ [ 0 ,   T ] } , where X = Surge and Y = Tension are responses simulated synchronously, see Section 3.1, with T being return period of interest, and let F X Y ( ξ , η ) = P r o b ( X ≤ ξ ,     Y ≤ η )   be the joint bivariate cumulative distribution function (CDF) of ( X ,   Y ) .

Raw response time series with d t = 0.0246 sec were blocked into k = 10 consecutive discrete time points maxima to reduce the neighbouring data points correlation effect. The conditioning level k was chosen according to the response power spectral density as it was observed that ACER functions have converged at that level in the distribution tail, see Section 3.2.

F X ( ξ ) and   F Y ( η )   denote the corresponding univariate marginal CDFs for X   and Y , respectively. In this paper, it is assumed that the bivariate couple ( Surge ( t ) ,   Tension ( t ) )   has been observed over a period of time   t ∈ [ 0 ,   T ˜ ] , where the observation duration   T ˜ is not long enough for accurately predicting the univariate extreme response levels with a target low probability of interest. Now, consider the case when a 'long' record of Tension ( t )   is available over a time   t ∈ [ 0 ,   T ] , with T " T ˜ , with a corresponding estimated CDF   F X l o n g ( ξ ) of the CDF F X ( ξ ) , which has a probability density function (PDF) p X = F X ′ . In this paper T return ≈ 5 years is the return period that corresponds to the extreme probability p of interest, p = 0.5 Δ 10 − 7 , T = 20 h, T ˜ =   T / 600 see Section 4. Then for any Y - response level of interest η ∗ , w i t h   Δ → 0 , F Y ( η ∗ ) =   F X Y ( ¥ , η ) = ∫ 0 + ∞ Prob ( Y ≤ η ∗ | X = ξ ) p X ( ξ ) d ξ = ∫ 0 + ∞ Prob ( Y ≤ η ∗ , X ∈ [ ξ ,   ξ + Δ ] ) Prob ( X ∈ [ ξ ,   ξ + Δ ] ) p X ( ξ ) d ξ =   ∫ 0 + ∞ F X Y , X ′ ( ξ , η ∗ ) d ξ ,     (9) with F X Y , X ′ ( ξ , η ) = ∂ ∂ ξ F X , Y ( ξ , η ) . See (Gaidai et al., 2019a)- (Xu et al., 2019b) for the details related to these equations. The following copula model for the bivariate extreme value distribution is referred to as the Gumbel logistic model (Wang, 2001)- (Xu et al., 2019b):     F X Y ( ξ , η ) = exp { − [ ( − ln   F X ( ξ ) ) m + ( − ln   F Y ( η ) ) m ] 1 / m } (10)

In this model, it is seen that m = 1 corresponds to the case when X and   Y are independent. When 0 < m < 1, X and   Y become dependent, (Wang, 2001)- (Xu et al., 2019b). However, this dependence structure is of a special kind since it only involves the marginal distributions. Still, it appears to be useful in some practical cases. Another popular extreme value copula is the Asymmetric Gumbel logistic model, see (Wang, 2001)- (Xu et al., 2019b). This bivariate extreme value model will not be discussed further in this paper since it gives identical results as the Gumbel logistic model. The Gumbel logistic model has been verified to be useful for various offshore engineering practical applications, provided the marginal extreme value distributions are estimated using the univariate ACER method instead of standard asymptotic extreme value distributions (Wang, 2001)- (Xu et al., 2019b). This is because the asymptotic distributions typically used are often not accurate enough in the tails, being the results of fitting real, subasymptotic data to asymptotic distributions. If Eqn. 35 is differentiated with respect to x, it is obtained that F X Y , X ′ ( ξ , η ∗ ) = F X Y ( ξ , η ∗ ) [ 1 + ( ln   F Y ( η ∗ ) / ln   F X ( ξ ) ) m ] 1 m − 1 d d ξ ln ⁡ F X ( ξ ) (11)

The numerical estimates F ^ Y ( η ∗ ) and F ^ X Y ( ξ , η ∗ ) of F Y ( η ∗ ) and F X Y ( ξ , η ∗ ) , respectively, based on the available time series of recorded data, are now used in the following expression to obtain the corrected estimate F Y c o r r ( η ∗ ) : F Y c o r r ( η ∗ ) = ∫ 0 + ∞ F ^ X Y ( ξ , η ∗ ) [ 1 + ( ln   F ^ Y   ( η ∗ ) / ln   F ^ X l o n g ( ξ ) ) m ] 1 m − 1 d d ξ ln ⁡ F ^ X l o n g ( ξ ) d ξ . (12)

Note that all quantities on the right side of Eqn. 12 are known from the available time series of recorded data. The Gumbel copula parameter m has been calibrated to fit joint empirical distribution F ^ X Y . For the latter optimisation task, the Trust-region-reflective non-linear least-squares optimisation algorithm can be used, and the interior-point algorithm to find the minimum of constrained non-linear multivariable function; for details, see (Wang, 2001)- (Xu et al., 2019b).