Assessing some statistical and physical modelling uncertainties of extreme responses for monopile-based offshore wind turbines, using metocean contours

This study examines the influence of probabilistic models for wave parameters in the joint environmental model and hydrodynamic/soil models on extreme mudline bending moments for monopile-based wind turbines at representative wind speeds, using the environmental contour method. For significant wave height, the 3-parameter Weibull model using the method of moments (MoM) provides the best fit to hindcast data across different wind classes, for the statistical models and data considered in the study. The hybrid Log-normal-Weibull (LonoWe) model also provides a reasonable fit but is sensitive to the transition point between distributions. Both models yield the largest extreme responses, with differences of approximately 0.5–3.5%. The 3-parameter Weibull model with maximum likelihood estimation (MLE) and the 2-parameter Weibull model result in less conservative contours, leading to up to 13% lower extreme responses, compared to LonoWe and Weibull (MoM). Regarding peak period, both the Log-normal and 3-parameter Weibull models provide reasonable fits, with the latter being more accurate near the steepness (breaking) limit. The stochastic variation among maxima due to seed variability and the uncertainty in quantile estimates as a function of number of samples was found to be crucial, particularly for severe sea states at the cut-out speed. Soil modelling is particularly important when the turbine is parked and encounters peak wave periods close to the turbine’s natural periods, while the effect of soil modelling on the extremes during turbine operation is negligible. Additionally, the impact of diffraction becomes relatively important for short wave periods. However, it is worth noting that the choice of load models has less impact on extreme responses compared to variations in the contours caused by different statistical models or seed variability.


Introduction
Offshore wind turbine (OWT) foundation design shall ensure structural integrity against extreme environmental conditions (IEC, 2019).Ultimate limit state (ULS) design checks require estimating long-term extreme responses, represented by characteristic values with specified target return period (IEC, 2019; DNV-GL (Det Norske Veritas -Germanischer Lloyd), 2019).One approach is the full long-term analysis (FLTA) (Naess and Moan, 2012), which considers the variability in the environment and the short-term response, accounting for the contribution of all short-term conditions to the long-term response.The method directly integrates the product of the joint probability density function of a given environmental condition and the corresponding short-term response probability distribution.FLTA is computationally prohibitive, as it requires all combinations of environmental variables to be considered, i.e. in time-domain simulations or model tests, to capture the stochastic nature of the relevant responses (Baarholm et al., 2010).Furthermore, it does not account for serial correlation in environmental conditions, resulting in positive bias in estimated extreme responses (de Hauteclocque et al., 2022;Mackay et al., 2021).For OWTs, the wind shall be considered in addition to waves, increasing the number of short-term conditions to be analysed (Agarwal and Manuel, 2009;Li et al., 2016Li et al., , 2019)).Therefore, alternative approaches are commonly applied to estimate extreme responses, such as the environmental contour method (ECM) (Winterstein et al., 1993;Haver and Winterstein, 2009).
ECM is a simplification of general inverse reliability methods that include the response conditional on the environment as an additional variable.The method is widely used in offshore structural design, and it is computationally efficient as it decouples the probabilistic https://doi.org/10.1016/j.apor.2024.103880Received 16 June 2023; Received in revised form 6 December 2023; Accepted 2 January 2024 description of the environment from structural design, approximating long-term extreme responses by considering only a few critical shortterm conditions along the N-year contour.ECM typically consists of constructing a joint distribution of metocean parameters, and then establishing environmental contours.Different approaches exist for both steps (described in Section 2.1), introducing uncertainties on the resultant contours, and therefore, on the sea states to be analysed for response estimation.The present study uses the global hierarchical modelling approach (Haver, 1985;Haver and Nyhus, 1986;Haver, 1987;Mathisen and Bitner-Gregersen, 1990a;Bitner-Gregersen and Haver, 1991a;Bitner-Gregersen, 2015a;Horn et al., 2018), and Inverse first-order reliability method (IFORM), as those approaches are recommended in current engineering design guidelines (IEC, 2019;DNV-GL (Det Norske Veritas -Germanischer Lloyd), 2019;NORSOK, 2017).
Furthermore, using ECM on OWTs poses additional challenges, discussed in Section 2.2.Firstly, OWTs are dynamically sensitive structures, and their response strongly depends on both wind and wave conditions.Therefore, neglecting the stochastic nature of metocean variables, e.g. of the peak period, introduces further uncertainty on resultant responses (Haselsteiner et al., 2022;Ross et al., 2020).Secondly, OWTs have non-monotonic behaviour to wind loads, and extreme responses can occur when maximum wind loads, i.e. around the rated speed, are combined with moderate but steep waves, close to the OWT natural period (Valamanesh et al., 2015;Velarde et al., 2019).This effect becomes increasingly important for larger OWTs, due to the combined effect of higher wind loads, and their larger natural periods that are close to primary wave periods.Consequently, high dynamic excitation can lead to extreme responses for conditions with low peak period/wave height, and high steepness.Therefore, representative joint distribution models and contours are needed for a wide range of metocean combinations, and not only for the tail of the distributions.
Having identified the critical sea states along the contours, the magnitude of extreme responses is also largely dependent on the physical load models that are used in the simulation tools.The present study focuses on the impact of hydrodynamic load and soil-structure interaction modelling on the extreme responses in the foundation.
Morison's equation (Morison et al., 1950) is widely used for hydrodynamic loads, assuming slender structures.For structures with large dimensions relative to the incident wave wavelength, diffraction effects become important, and design standards (IEC, 2019) suggest using MacCamy and Fuchs' model (MacCamy and Fuchs, 1954), which can more precisely estimate hydrodynamic loads for larger monopile foundations.Nevertheless, MacCamy and Fuchs is strictly applicable only to linear waves, while extreme responses are typically governed by nonlinear steep waves.Approaches such as higher order potential flow simulations or computational fluid dynamics, which are able to capture both near-field diffraction and nonlinear wave kinematics, are too computationally expensive for use in aero-hydro-servo-elastic simulation tools for assessing thousands of design load cases.To tackle this limitation while maintaining computational efficiency, a new load model that incorporates a frequency-dependent mass coefficient and can be applied with higher-order wave kinematics has been developed and validated against experimental results for irregular severe sea states (Dadmarzi et al., 2024).In the present study, the relative significance of diffraction on extreme responses is assessed by comparing the conventional Morison's equation with the new load model.
Soil-structure interaction between the foundation -the part of the monopile below mudline -and the surrounding soil (hereafter referred to as foundation modelling) for the design of monopile-based OWTs is often modelled using  −  curves, as outlined in design standards (IEC, 2019;DNV-GL, 2018).In this approach the pile is represented as a beam, and the relation between the soil resistance  to lateral displacement  is modelled using a series of discrete, uncoupled, nonlinear elastic springs along the pile.This method was primarily developed for slender, small-diameter piles with high lengthto-diameter ratios, subject to dominant axial loads (Reese and Matlock, 1956).Despite being widely employed (Leblanc et al., 2010;O'Kelly and Arshad, 2016), the method's adequacy to fully capture the magnitude and character of the loads of monopile foundations has been questioned (Leblanc et al., 2010;Doherty and Gavin, 2012;Bhattacharya, 2014).It also lacks the ability to model important aspects such as hysteretic soil damping.Consequently, alternative numerical foundation models have been proposed for fully-integrated analyses of monopile-based OWTs (Page et al., 2016).In the present study, the relative importance of foundation modelling on extreme responses is assessed by comparing the commonly used  −  approach with a nonlinear elasto-plastic macro-element model that incorporates hysteretic effects, developed as part of the REDWIN project (Løkke, 2018), specifically for predicting the response of monopile-based OWTs in integrated time-domain analyses (Page et al., 2017(Page et al., , 2018b(Page et al., ,a, 2019)).
Using ECM involves various statistical uncertainties that affect the long-term extreme responses, with particular challenges when used for OWTs.Concurrently, due to the continuous growth in the size of monopiles, commonly used physical load models introduce uncertainties in the extreme response estimates.The main aim of this study is to investigate the impact of some uncertainties related to statistical and physical modelling on extreme responses of large diameter monopile-based OWTs, and to provide a more thorough understanding of their relative importance, when characteristic design responses are estimated using ECM.For that purpose, two monopile-based OWT models were considered, the DTU 10 MW (Bak et al., 2013) and the IEA 15 MW (Gaertner et al., 2020) turbines.The focus is on understanding the variability in extreme responses when employing metocean contours derived from different probabilistic models for joint distributions.Additionally, state-of-the-art foundation and hydrodynamic load models, better suited for large-diameter monopiles, are compared with models commonly used for small-diameter piles.
The paper is organized as follows: Section 2 gives a more thorough literature review on the environmental contour method, and the challenges related to the hydrodynamic load and foundation models used for OWTs.Section 3 summarizes the hindcast data, the joint environmental model and the procedure to establish the contours used in the study.Section 4 presents the OWT simulation models, and the load models compared in the study.Finally, Section 5 summarizes the results and Section 6 concludes the paper with recommendations for future work.

Environmental contours
Establishing environmental contours typically firstly involves estimating a statistical model to describe the joint distribution of metocean parameters.This model is then used to establish the -year contours, which represent combinations of parameters with -year return period.It is worth noting that there are also methods that do not require a joint distribution model to estimate contours (Derbanne and de Hauteclocque, 2019;Mackay and de Hauteclocque, 2023).Once the contours are established, the relevant response is estimated for various conditions along the -year contour.Typically, multiple 1-hour or 3hour stochastic time-domain simulations are carried out for different sea states along the contours, and the characteristic response in each sea state is then defined as the 85 th -95 th percentile of the distribution of the (1-hour or 3-hour) maximum response conditional on sea state to indirectly account for the short-term variability (Baarholm et al., 2010;Haver et al., 1998;Kleiven and Haver, 2004;Muliawan et al., 2012).Alternative approaches (NORSOK, 2017;Derbanne et al., 2017) accounting for short-term variability are not discussed here.Finally, the method assumes that the -year maximum response can be approximated by the highest characteristic response amongst the conditions along the contour.
Various methods can be employed to determine the joint distribution of metocean variables, including global hierarchical models (Haver, 1985;Haver and Nyhus, 1986;Haver, 1987;Mathisen and Bitner-Gregersen, 1990a;Bitner-Gregersen and Haver, 1991a;Bitner-Gregersen, 2015a;Horn et al., 2018), copula models (Vanem, 2016;Zhang et al., 2018;Fazeres-Ferradosa et al., 2018;Lin et al., 2020), kernel density estimates (Ferreira and Guedes Soares, 2002;Eckert-Gallup and Martin, 2016), and conditional extreme value models (Jonathan et al., 2010(Jonathan et al., , 2014)).The global hierarchical model that is applied in the present work is widely used and also recommended by the design standards, while it poses several challenges.Firstly, the selection of specific statistical models for the metocean parameters and their dependence structure is mostly based on subjective user-judgment, rather than solid theoretical principles.Furthermore, fitting a statistical model to all observations does not guarantee a representative fit to the tail of the distribution, which is the range of most interest for estimating extremes, while serial correlation in the data is ignored, resulting in positive bias in extreme estimates (de Hauteclocque et al., 2022;Mackay et al., 2021).
Similarly, different methods exist for establishing metocean contours, based on different concepts of exceedance probability.For IFORM (Giske et al., 2017;Winterstein et al., 1993;Haver and Winterstein, 2009), and the direct sampling method (Bang Huseby et al., 2013;Huseby et al., 2015;Vanem, 2019), the contour exceedance probability is defined as a marginal exceedance probability under various rotations of the coordinate axes.In contrast, for recent methods such as inverse second-order reliability method (ISORM) (Chai and Leira, 2018;Giske et al., 2018), and highest density contour method (HDC) (Haselsteiner et al., 2017), contours are defined in terms of the total exceedance probability, being significantly more conservative (Mackay and Haselsteiner, 2021).Another approach, referred to as Direct-IFORM, replaces the multi-dimensional joint distribution by a series of univariate fits, which can be applied in higher dimensions for arbitrary numbers of variables, without loss in performance.The method was firstly introduced by Derbanne and de Hauteclocque (2019) and further presented by Mackay and de Hauteclocque (2023), where its applicability to the design of offshore wind turbines is discussed.It is based on the same notion of marginal exceedance probabilities under rotations of the axes, that is used in standard IFORM and direct sampling contours, and its main advantage is that there is no need to establish a joint distribution to estimate contours, as the marginal exceedance probabilities can be estimated directly from the data.
Several studies have compared methods for establishing metocean contours (Manuel et al., 2018;Vanem, 2017;Huseby et al., 2017;Vanem et al., 2020;Raed et al., 2020).Mackay and Haselsteiner (Mackay and Haselsteiner, 2021) studied the impact of the choice of contour on simple structural reliability problems, highlighting that some understanding about the shape of a structure's failure surface is required to choose an appropriate contour method.Haselsteiner et al. (2021) presented an extensive study comparing contours established from various methods.Large variations were observed with no clear outcomes regarding how the different steps followed to establish the contours contribute to these variations.de Hauteclocque et al. (2022) extended the study, and compared extreme responses obtained from various contour methods to reference estimates for different offshore structures.The study highlighted the importance of statistical model accuracy for accurate response estimates, while it was found that serial correlation introduces a significant positive bias into long-term response estimates, especially for lower return periods.Ross et al. (2020) presented an extensive survey of the metocean user community regarding metocean contours, with recommendations about when and how they should be used.

Challenges of ECM related to OWT responses
ECM has been used for OWTs; however, its direct applicability is challenging.IEC 61400-3-1 (IEC, 2019) suggests establishing a 2dimensional (2D) wind speed (  ) -significant wave height (  ) contour, and the peak period (  ) is chosen deterministically, as the value leading to the largest response for a particular   −   combination.For dynamically sensitive structures, neglecting the stochastic nature of the wave period in establishing the contours can lead to non-conservative extreme response estimates.As stated by Ross et al. (2020), including all dominant metocean variables is crucial to ensure accurate -year extreme response estimation.Valamanesh et al. (2015) calculated larger mudline bending moments for moderate sea states close to the natural period, compared to load cases where   was chosen for the most severe   −   combination.Moreover, Velarde et al. (2019) demonstrated that resonant loads during parked or idling situations under operational wind speeds could dominate the design loads for monopile-based OWTs.The observation that critical conditions for offshore structures may be in the vicinity of the natural period has also been mentioned or discussed in other studies (Winterstein et al., 1993;Vanem, 2017).Different approaches have been used to account for the three dominant variables (  −   −   ) in a probabilistic manner.Li et al. (2015) proposed a 3-dimensional (3D) contour surface, Horn and Winterstein (2018) used multiple 2D contours, (  −   and   −   ), while other studies used   −   contours representative for different wind speeds (Li et al., 2016).It is still unclear which approach is best suited to deal with the environmental variables relevant to the design of OWTs (Haselsteiner et al., 2022).
Another challenge using the ECM for OWTs is their non-monotonic behaviour to wind loads.In particular, OWTs are designed with controllers that maximize power extraction below the rated wind speed, while between rated and cut-out speed the blades are pitched to maintain constant power output and minimize the mean thrust load.The mean wind loads gradually decrease until the cut-out speed.Beyond the cut-out speed, the blades are fully pitched, and the turbine enters a parked state where wind loads and fore-aft aerodynamic damping abruptly drop.This behaviour results in maximum wind loads near the rated speed, and consequently extreme responses, such as overturning moments, can occur in operational conditions.The importance of evaluating operational and parked states has been highlighted in the literature, e.g., by Saranyasoontorn and Manuel (Saranyasoontorn and Manuel, 2004).Li et al. (2016) proposed a modified ECM by establishing multiple   −   metocean contours for different wind speeds between rated and cut-out wind speed to find the largest response.The method has been applied to bottom-fixed and semi-submersible OWTs (Li et al., 2016(Li et al., , 2017) ) and integrated offshore renewable energy devices (Li et al., 2019), with reasonable accuracy, compared to FLTA.The term ''modified ECM'' refers to the approach of Li et al. (2016) where multiple   −  metocean contours are employed over the operational range of the wind turbine, but does not imply any modification to the process of obtaining these contours.

Foundation modelling of OWTs
The  −  method is widely employed for OWTs and can be easily integrated into simulation tools.Contrary to the axially loaded piles for which  −  curves were developed initially, the loading of monopiles in OWTs is characterized by horizontal thrust loads, which, combined with the height of the structure, lead to dominant bending moments at the seabed, whereas the vertical loads are relatively unimportant (Byrne and Houlsby, 2004;O'Kelly and Arshad, 2016).Furthermore, large-diameter monopiles are characterized by slenderness ratios (ratio between embedded length and diameter of monopile) below 5, which leads to relatively rigid behaviour compared to the slenderness ratios of 15-20 and flexible behaviour assumed in the development of the − method (Lombardi et al., 2013).Another essential aspect of foundation modelling is accurately representing soil damping, which arises from various sources (Page et al., 2016;Bhattacharya, 2019).Among these sources, hysteretic damping becomes the primary contributor when analysing monopile foundations (Page et al., 2017).The relative importance of soil damping notably increases in situations where aerodynamic damping is negligible, i.e. when the OWT is parked, therefore its accurate representation in the simulation tools becomes of high importance.
An alternative approach is macro-element modelling, which was initially introduced by Roscoe and Schofield (1957) and has been widely applied to offshore applications (Bienen et al., 2006;Correia, 2011;Houlsby and Cassidy, 2002;Skau et al., 2017;Tistel et al., 2017).Unlike  −  curves, macro-element models condense the response of the foundation and surrounding soil into a force-displacement relationship at a single point (Correia, 2011).A macro-element was developed (Løkke, 2018) specifically for predicting the response of monopile-based OWTs (Page et al., 2017(Page et al., , 2018b(Page et al., ,a, 2019).In the model, the relationship between horizontal forces, moments, and displacements, as well as rotations, is described using the framework of multisurface plasticity, which was originally proposed by Iwan (1967).Details about how hysteretic behaviour is incorporated in the model are given in Page et al. (2018b).The complete nonlinear response of the soil-foundation system can be captured through load-displacement curves obtained from nonlinear pushover analyses conducted in FEA.This macro-element model incorporates soil resistance components, including base and side shear.Furthermore, unlike the traditional  −  formulation, it is capable of reproducing the change in overall stiffness due to nonlinear hysteretic behaviour and the hysteretic damping observed in monopile-based OWTs in integrated time-domain analyses.The model also accounts for the influence of multi-directional loading, which has been observed to affect the foundation stiffness and hysteretic damping (Page et al., 2018a).To validate its performance, the model's behaviour was compared against results from large-scale pile tests, full-scale field measurements, and Finite Element Analysis (FEA) simulations conducted on monopile foundations of an OWT installed in the North Sea (Page et al., 2018b(Page et al., ,a, 2019)).The comparison demonstrated that the macro-element model can accurately predict the nonlinear hysteretic behaviour observed in field tests and simulated in FEA.

Hydrodynamic modelling
Different methods exist to estimate extreme loads on vertical cylinders.Those include the constrained wave approach, where a regular stream function wave is embedded into a background irregular linear wave realization, or fully nonlinear computations (either considering potential flow or including viscous effects in computational fluid dynamics).However, these methods have limitations in terms of accuracy or computational time (Pierella et al., 2021;Suja-Thauvin, 2019).Several higher-order theories have been developed, such as Faltinsen et al. (1995), Malenica and Molin (1995), and Kristiansen and Faltinsen (2017).Comparisons of these higher-order models and Rainey's model (Rainey, 1989) for extreme responses of monopile foundations can be found in the literature (Suja-Thauvin, 2019;Bredmose et al., 2013) but are not considered for the premise of the study.Among the models for hydrodynamic loads, Morison's equation (Morison et al., 1950) is recommended by design standards (IEC, 2019; DNV-GL (Det Norske Veritas -Germanischer Lloyd), 2019), and it is commonly used in practice for OWT global analysis in aero-servo-hydro-elastic simulations tools.Morison's equation is a semi-empirical formulation based on the assumption that the cylinder is slender, meaning its diameter is small compared to the wavelength of the incident waves.This assumption implies that the cylinder does not generate a significant diffracted wave field.
For structures with dimensions large relative to the incident wavelength (typically   > 0.2), significant diffraction effects occur, and the application of Morison's equation alone may not provide accurate results.In such cases, design standards (IEC, 2019) recommend using the analytical solution developed by MacCamy and Fuchs (1954), which addresses the diffraction of long-crested waves on vertical piles.The computation of the total horizontal force follows a similar approach to Morison's equation, with the drag component unchanged.However, the inertia force component considers the size of the pile and the frequency of the incoming waves, accounting for the diffraction effects.Notably, as monopile foundations increase in size, diffraction effects become more prominent.
Morison's equation allows for the use of wave kinematics from various theories, whereas the MacCamy and Fuchs model is strictly applicable only to linear waves.Although this is not problematic for fatigue analysis, which is mainly influenced by moderate waves well represented by linear wave theory, the maximum loads relevant to ULS design are typically governed by steep waves.In these cases, nonlinearities become significant, especially in intermediate and shallow water depths, and cannot be ignored.A new model that combines traditional Morison's formulation with a frequency-dependent mass coefficient (based on the formulation from MacCamy and Fuchs) has been developed (Dadmarzi et al., 2024).This computationally efficient model approximates diffraction effects for larger diameters when higher wave kinematic theories are employed.

Hindcast data
A numerical hindcast model from the National Kapodistrian University of Athens (NKUA) was used to generate 10-year hindcast data for several locations for the Marina Platform project (Li et al., 2015).The hindcast data have 1-hour resolution for the period 2001 to 2010 for a site located at the Norwegian Continental Shelf (55.11 • N, 3.47 • E) and 30 m water depth.The dataset provides information about metocean parameters such as mean wind speed 10 m above sea level ( 10 ), significant wave height (  ), wave peak period (  ) and wind-wave directionality.The mean wind speed at the hub height ( ℎ ) of each turbine used in the study was estimated using the power law for wind shear with exponent  = 0.14 (IEC, 2019).The combined sea data, which include wind-driven and swell components, have been used.Fig. 1 shows scatter density plots for the main metocean parameters.Current and directional variability were not considered.
The metocean data were assumed to be independent and identically distributed.While this assumption is commonly used, it is more applicable when the data are derived from simulations such as Monte Carlo methods using the established joint environmental model, noting that even then, they are still not representative of the environment.When fitting a parametric model to measured or hindcast data, as done in this study, this assumption is not appropriate, as winds and waves are properties of large-scale physical processes, and always exhibit serial correlation.The autocorrelation function for   during different winter seasons (November-February) showed significant temporal serial correlation, dropping below 0.8 only after 10 h of data.Neglecting this serial correlation can potentially overestimate long-term extreme responses (de Hauteclocque et al., 2022;Mackay et al., 2021).However, since the goal of this study is to compare and evaluate statistical and modelling uncertainties in extreme responses using a specific dataset, rather than providing design values for OWTs, the strict validity of this assumption is not expected to significantly affect the study's findings.

Joint environmental model and dependence structure
To establish environmental contours, a joint probabilistic model was fitted to the metocean data, by adopting the conditional modelling approach (Bitner-Gregersen and Haver, 1989).Wave data were sampled in wind speed classes with bin size of 2 m/s, associated with the mean wind speed at hub height,  ℎ .
A marginal distribution was fitted to the hub height wind speed data, denoted as   ℎ (), and conditional models were fitted to   for given  ℎ , and to   for given  ℎ and   , denoted as    | ℎ (ℎ|) and    |  , ℎ (|ℎ, ), respectively.The joint probability is expressed as the product of marginal and conditional probability density functions, as given in Eq. ( 1), (𝑡|𝑢, ℎ). (1)

Probabilistic models for metocean parameters
The 1-hour mean wind speed  ℎ is modelled using a 3-parameter Weibull distribution defined by Eq. ( 2).
Parameters   ,   and   denote the shape, scale and location parameters, respectively.A 2-parameter model (  = 0), which is more common in the literature, was also examined (Bitner-Gregersen and Haver, 1991b;Bitner-Gregersen, 2005, 2015b).Good agreement between the data and the fitting models was obtained, with negligible difference for wind speeds lower than 3.0 m/s, which are not considered important for the purposes of the study.For significant wave height,   , previous studies have shown that it can be modelled reasonably well by a 3-parameter Weibull distribution (Bitner-Gregersen and Haver, 1991b;Mathisen and Bitner-Gregersen, 1990b).In this study, three statistical models have been adopted for the conditional distribution of   given  ℎ : 1.A 2-parameter Weibull distribution, given in Eq. ( 3), where  2 ,  2 are the shape, and scale Weibull parameters, respectively.2. A 3-parameter Weibull model, given in Eq. ( 4), where  3 ,  3 , and  3 denote the shape, and scale and location Weibull parameters, respectively.3. A hybrid model consisting of a Log-normal and a Weibull distribution, known as the LonoWe distribution (Haver, 1985), The LonoWe distribution was first introduced and is commonly used with a 2-parameter Weibull distribution.In the present study, a 3-parameter model was used.For each wind class,    | ℎ (ℎ|) was modelled by a Log-normal distribution for ℎ ≤  and a 3-parameter Weibull distribution for ℎ >  (Eq.( 5)).In Eq. ( 5),  ℎ is the mean and  2 ℎ is the variance of ln(ℎ).Parameters  ℎ ,  ℎ and  ℎ denote the shape, scale and location parameters of the Weibull distribution, respectively.
The dependence function for the scale, shape, and location parameters for Weibull models are described by the power function given the wind speed data within a wind class, shown in Eq. ( 6), where parameters  =1,2,3 ,  =1,2,3 ,  =1,2,3 are estimated using nonlinear curve-fitting based on least-squares, and  represents the wind speed data within each wind class.Fig. 2 shows an example of the fitted parameters of the dependence function for the Weibull models.The dependence functions for the Log-normal parameters of the LonoWe distribution are given in Eq. ( 7), where  =1,2,3 ,  =1,2,3 are the parameters of the power and exponential functions, respectively.It should be noted that the dependence functions used in this study are typically used in literature and recommended in the design standards (IEC, 2019; DNV-GL (Det Norske Veritas -Germanischer Lloyd), 2019).However, they do not provide any physical interpretation.Haselsteiner et al. (2020) present recent work on physical interpretations of the dependence functions.
Based on previous studies, when only wave data are considered, the conditional distribution for   given   is usually assumed to follow a Log-normal distribution (Haver, 1985;Mathisen and Bitner-Gregersen, 1990a;Bitner-Gregersen andHaver, 1989, 1991b).A method to determine the conditional distribution for   given   and   was developed by Johannessen et al. (2002), and it has been used in other studies (Horn et al., 2018;Li et al., 2015).As mentioned by Li et al. (2015), the process of obtaining the distribution of   conditionally on both   and   is challenging and not straightforward following the methods described by Johannessen et al. (2002), therefore Johannessen's parametrization method was not used.In the present study, two statistical models were evaluated for the   conditional distribution given  ℎ and   using the available data: a Log-normal (Eq.( 8)) and a 3-parameter Weibull distribution (Eq.( 9)).
In Eq. ( 8),   and  2  are the mean and variance of ln() of the Lognormal model and   ,   ,   denote the shape, scale and location Weibull parameters, respectively.The dependence functions used for the peak period are the same as shown in Eqs. ( 6) and ( 7), where  is substituted with the wave height, ℎ, within the corresponding wind and wave class.

Parameter estimation
The parameters of the marginal distribution of  ℎ were estimated using the maximum likelihood estimate (MLE).
For the   2-parameter Weibull model, MLE was used to estimate the distribution parameters.To estimate the 3-parameter Weibull model parameters for   , two methods were used, namely MLE and Method of Moments (MoM).To estimate the distribution parameters for the MoM method, the expected value (   ), the variance ( 2

𝐻 𝑠
) and the skewness coefficient (   ) of   data for each wind class are estimated and then replaced in Eqs.(A.1), (A.2), and (A.3) respectively, given in Appendix.
For the LonoWe model, the following procedure was followed.For different values of transition point , first a Log-normal distribution was fitted for ℎ ≤  within the wind class.Then, the Weibull parameters were estimated using non-linear least square method, by solving the continuity condition for    | ℎ (ℎ|) and    | ℎ (ℎ|) between the adopted Log-normal and Weibull distributions at ℎ = .The LonoWe models established for different , were then compared and the distribution that minimized the root mean square error (RMSE) between the model and the data was adopted for further comparisons.
For the   Weibull and Log-Normal models, the MLE method was used to estimate the distribution parameters.

Metocean contours
Based on the established joint distribution, a contour corresponding to a given return period was estimated for each wind class.In the present study, the estimation of contours is based on the IFORM (Winterstein et al., 1993) using the Rosenblatt transformation (Rosenblatt, 1952), where a joint probability model in the physical space (X-space) is transformed to the standard Gaussian normalized U-space.The environmental variables in the physical space ( ℎ ,   ,   ) are mapped to uncorrelated normally distributed variables ( 1 ,  2 ,  3 ) based on Eq. ( 10), where  is the standard normal cumulative distribution function.Functions   ℎ ,    | ℎ ,    | ℎ ,  denote the cumulative distribution functions (CDFs) of  ℎ ,   ,   , respectively.Eq. ( 10) establishes a unique connection between variables (, ℎ, ) and ( 1 ,  2 ,  3 ), and the U-space variables can be back-transformed to the original physical space based on Eq. ( 11).
In the Gaussian space, points of constant probability density define a sphere with radius equal to  (the reliability index), shown in Eq. ( 12), and as a result, the   −   contours for each wind class  are based on a circle (in  2 −  3 plane) with radius   given in Eq. ( 13), The variable   is chosen as the mean value of the wind speed data in class .Assuming data to be statistically independent, highlighting that this is not true in reality, parameter  is related to the target exceedance probability   required to estimate the contours by Eq. ( 14), In Eq. ( 14),   is the target return period (expressed in total number of hours of stationary sea states) and  is the duration of a stationary sea state (Kleiven and Haver, 2004).Considering a 50-year return period   and the duration of a sea state event  = 1 hr based on the hindcast data, Fig. 3 shows the contour sphere of radius  = 4.58 and contour circles of radius   for various wind classes .

Assessment of contours compared to wave steepness limit
The physical limit of the wave steepness is generally not considered for the establishment of the probability distributions and the environmental contours.Therefore, non-physical   −   combinations can be present in the contours, resulting in non realistic loading conditions.This may result from multiple reasons, such as that the marginal and conditional models may be misspecified, i.e. they are not an accurate fit for the data or the assumed dependence structure for conditional model The average wave steepness for irregular sea states, denoted as   , is defined as the ratio of   to the peak wavelength (  ).The limiting values of wave steepness (DNV-GL (Det Norske Veritas -Germanischer Lloyd), 2019) are shown in Eq. ( 15), and linear interpolation can be used between the boundaries.
Wave length can be calculated from the deep-water dispersion relation, shown in Eq. ( 16).However, due to the shallow depth of the location ( = 30 m), the deep-water wave length cannot be used directly, and the intermediate depth dispersion relation should be used to estimate the wave length   , given in Eq. ( 17).
In the present study, the resultant contours are compared to the steepness curves as obtained using the deep-water and intermediate wavelength relationship, to evaluate the adequacy of the   conditional models mainly at the steep side of contours, and ensure that subsequent design sea states remain within realistic ranges.

Simulation models
Two monopile-based OWT models were used in the study: the DTU 10 MW (Bak et al., 2013) and the IEA 15 MW (Gaertner et al., 2020), with main characteristics shown in Table 1.The rotor-nacelle assembly, tower and monopile above seabed were modelled in SIMA (SIMO, 2017;RIFLEX, 2017), an aero-hydro-servo-elastic software developed by SIN-TEF Ocean.All wind inflow fields for the time-domain simulations were performed using the stochastic, full-field, turbulence simulator TurbSim (Jonkman, 2009), using the Kaimal spectrum and exponential coherence model (IEC, 2019).
The wind turbine blades were modelled using flexible beam elements, and the structural and aerodynamic coefficients were based on the reference document by Bak et al. (2013).The corresponding generator-torque and blade-pitch controllers were implemented based on the Basic DTU Wind Energy Controller (Hansen and Henriksen, 2013) and NREL's Reference OpenSource Controller (NREL, 2020), respectively.Aerodynamic loads were calculated using the blade element momentum theory, incorporating the Glauert and Prandtl corrections (Burton et al., 2011), and dynamic stall and wake effects.
The tower and monopile above the seabed were modelled using axisymmetric beam elements.The monopile had a constant diameter (  ) and wall thickness (  ) along its length for each wind turbine model.The dimensions of the monopile for the 10 MW turbine were based on the work of Velarde and Bachynski (2017), considering a water depth of 30 m and a diameter-to-thickness ratio of approximately 80.The monopile for the 15 MW wind turbine was designed to achieve a target natural frequency around 0.19 Hz, which is approximately in the middle of the soft-stiff frequency range for that turbine.Table 1 shows the monopile dimensions for the two OWTs.
The embedded length of the monopile was chosen to have a lengthto-diameter ratio of 4, similar to the 10 MW monopile, and a diameterto-thickness ratio of 100.For each turbine model, the tower was represented by sections of specified length with a constant diameter and wall thickness, increasing linearly from the top to the base of the tower.The sectional properties of the tower for the 10 MW OWT followed the preliminary design presented by Velarde (2015), considering a monopile with a 9 m diameter and 0.11 m thickness at a water depth of 30 m.The tower design for the 15 MW OWT was based on the original design (Gaertner et al., 2020).

Foundation modelling
Two modelling approaches were employed in the simulation tool to represent the foundation, i.e., the monopile below the mudline.The first model is the − curves, and the second is a macro-element model, both described in Section 2.3.Both models have been calibrated using the results of full 3D continuum modelling of the soil volume and the foundation.Details of the soil layer properties used for the calibration were given by Katsikogiannis et al. (2021) In the  −  curves method, the pile was modelled as a beam, and the resistance of the soil to lateral displacement of the pile was represented by a series of discrete, uncoupled, axisymmetric nonlinear elastic springs along the pile.The nonlinear elasto-plastic macro-element model condenses the soil-foundation system into a single node located at the mudline, representing the interface between the foundation and the rest of the structure.The model is integrated into SIMO-Riflex through a dynamic link library (DLL), which serves as an interface between the macro-element model and SIMO-Riflex.At each calculation time step, the DLL retrieves the nodal displacement (  ) and rotation fields (  ), at the mudline from SIMO-Riflex to estimate the external forces (  ,   ) acting on the macro-element model.These external forces are then computed by the macro-element model and returned to SIMO-Riflex.This iterative process continues until specific convergence criteria are satisfied.Fig. 4 shows the 15 MW OWT SIMO-Riflex model using  −  springs and macro-element to model foundation.
Fig. 5 shows the response at the mudline obtained from a decay test for the two foundation models, for the 15 MW wind turbine.The comparison highlights the non-linear elastic behaviour captured by both foundation models.While the  −  curves follow the same loading curve during load reversals, the macro-element model exhibits non-linear hysteretic behaviour, resulting in a loop.The area enclosed by the loop represents the energy dissipated in the soil, which can be interpreted as hysteretic damping at the foundation level, as also shown in Page et al. (2017).

Hydrodynamic models
Two hydrodynamic load models are employed to assess the significance of diffraction effects in extreme loads.The first model employs the conventional Morison's equation, and the second is new load model that incorporates a frequency-dependent mass coefficient and can be applied with higher-order wave kinematics (Dadmarzi et al., 2024).Morison's equation, combines viscous drag and inertia load components, and it can be expressed as follows,  where  is the force per unit length,   ,   are the drag and inertia coefficients respectively,  is the water density, and   is the pile diameter.  is the relative horizontal fluid particle velocity normal to the member (  =   −   ), and u , u are the flow and structure accelerations normal to the member, respectively.The inertia coefficient can be split in two contributions where   is the added mass coefficient.The first term represents the Froude-Krylov force from the pressure of the undisturbed wave, and the second term represents the diffraction effects.The drag and inertia coefficients can vary based on different factors, and experimental results have shown considerable variation in their values even under similar conditions (IEC, 2019).For design calculations, it is common practice to assume constant values of   = 0.9 and   = 2.0, respectively (DNV-GL (Det Norske Veritas -Germanischer Lloyd), 2019).Notably, as monopile foundations increase in size, diffraction effects become more prominent.This is illustrated in Fig. 6, where the inertia load coefficient   exhibits a significant decrease for shorter waves and larger diameters, resulting in a reduction of the total horizontal force on the structure.To compute the inertia force per unit length, the undisturbed incoming fluid particle acceleration time-series obtained using second-order waves is transformed to the frequency domain using Fourier transform.Assuming that the inertia force can be represented as a linear combination of  Fourier components derived from the nonlinear acceleration signal, the force is given by Eq. ( 20) Each Fourier component of the acceleration signal, u , is associated with a corresponding frequency-dependent mass coefficient,    , calculated using MacCamy-Fuchs diffraction theory as shown in Eq. ( 21), where are the first order Bessel functions of the first and second kind, respectively, with the prime denoting the derivative with respect to the argument,  =   ∕2 is the monopile radius, and   is the wave number.For a Fourier component at frequency   ,   is the positive real solution calculated from the linear dispersion relation (Eq.( 22)) for intermediate water depth ℎ:

Results and discussion
The results section of the paper is divided into three main sections.Sections 5.1 and 5.2 examine the impact of the conditional statistical models for   and   , described in Section 3.2.1, on the resulting contours.The objective is to understand how different statistical models affect the contours used to estimate extreme responses.Then, in Section 5.3, the variation of extreme responses across representative wind classes using different statistical models is explored.The influence of stochastic variability among realizations is also evaluated in this context.Finally, Sections 5.4 and 5.5 investigate the impact of foundation and hydrodynamic load modelling on extreme responses, for critical sea states.

Effect of 𝐻 𝑠 conditional distributions on resultant contours
As mentioned in Section 3.2, wave data were sampled in wind speed classes, associated with the mean wind speed at hub height,  ℎ .Due to the different hub height between the turbines (see Table 1), the data contained within the same wind class may slightly differ between the 10 MW and 15 MW turbine models.Consequently, the fitted conditional distributions and resulting contours also exhibit slight variations between the two turbines.However, similar conclusions apply to both cases.
For the   conditional distribution given  ℎ , different statistical models were compared.Fig. 7 shows the empirical and the fitted distributions for different wind classes.The logarithm of the exceedance probability, (1 −    ), is plotted in the  axis, to better capture the behaviour of the models at the tail probability of the distributions.The 2-parameter Weibull distribution resulted in large discrepancies between the analytical model and the hindcast data in the upper tail.This can be observed by the steep slope in the tail regions, indicating that the probability of observing values beyond a certain threshold decreases quickly, and the model is significantly non-conservative.The discrepancies were notable particularly in lower and intermediate wind classes, as the model was inadequate to capture the mixed nature of the data (wind-sea and swell).For higher winds, it provided a more reasonable fit, as the sea states are mostly wind-driven for the site considered.
A more accurate fit of the upper tail was obtained using the 3parameter Weibull distribution; however, its behaviour depends on the method used to estimate the fitted parameters.Using the maximum likelihood estimate (MLE), less weight is given to the tail regions of the data to estimate the distribution parameters.This resulted in nonconservative contours disregarding some severe sea states, mainly for lower and intermediate wind classes, where many small and moderate sea states are observed, with only a few severe sea states.This is one of the main limitations of global hierarchical models used to fit all observations, when the upper tail is of main interest.In contrast, this limitation seemed to have less impact when using the method of moments (MoM), as it gives more weight to the upper-tail data to estimate the distribution parameters.Consequently, this model provided a better fit for the upper tail and more conservative contours for all wind classes.Generally, the MoM resulted in relatively high estimates for the location parameter ( ℎ ), providing a worse fit at the lower tail than the MLE method; however, this is not of great concern when evaluating the extremes.
The LonoWe model accurately fitted the full range of   data within most wind classes.However, its performance was sensitive to the transition point () between the Log-normal and Weibull distributions, which had a notable impact on the fit and the resultant contours, as shown in Fig. 8.This sensitivity was particularly important for the upper tail of the data.As mentioned in Section 3.3, the value of  was chosen based on minimizing the RMSE between the data and the analytical model.However, the variations in RMSE became negligible for increasing values of , making it challenging to estimate an optimal value.If an accurate probabilistic model for the whole range of   within a wind class is of interest, then it is essential to select a value for  that ensures an adequate amount of data is available for fitting both the Log-normal and Weibull distributions, avoiding extremely low or high values of .When extreme values of   are of interest, focus should be given on fitting only the high quantiles of   , see e.g.work from Haselsteiner and Thoben (2020), omitting the bulk of the data, while  can be chosen based on a quantity that better indicates the quality of the fit at the upper tail, and not RMSE, which is a measure of fit to all observations.Fig. 9 shows the differences between the fitted models for three wind classes and the resultant 50-year contours with the   −  scatter density plots.The Weibull model for   is used for the plots in this section.Overall, the data seem to lie inside the 50-year contours for the 3-parameter Weibull and LonoWe distributions without excessive out-crossings.In some cases, the models resulted in slightly higher (or lower) extreme wave heights than the hindcast data suggest.This is a consequence of model misspecification, the statistical uncertainty due to the limited number of observations of higher   classes for a given wind speed and the method used to estimate the distribution parameters.Furthermore, it should be noted that although the differences in the upper tail seem negligible in the Weibull probability plot, they are significant in the resultant contours.

Effect of 𝑇 𝑝 conditional distributions on resultant contours
The peak period conditional distribution is relatively challenging to establish, mainly for two reasons.Firstly, fitting an analytical model to combined sea data, which includes wind-driven and swell components, enhances the modelling difficulty, due to their different nature.Fig. 10(a) illustrates the scatter density plots of the site-specific hindcast data for wind-driven, and swell components and the 50-year contours for various wind classes.As depicted, for the site considered in the study, for low and moderate wind speeds, sea states are often of a combined nature with the most severe sea states originating from swell, while in high wind speeds severe sea states are purely wind-driven.The wave components that dominate in different wind speeds are also reflected in the shape of the contours (Fig. 10(b)), where for low wind speeds, contours are characterized by long periods representing mainly the swell and moderate wind-driven seas, while gradually their shape transforms representing the wind-driven severe sea states that dominate in high wind speeds.
Log-normal and Weibull models fitted the wave peak period data reasonably well.An example is shown in Fig. 11 for different wave classes in wind class 10-12 m/s, where both models depict similar behaviour with minor differences.Nevertheless, the variations in contours were distinct.Fig. 12 shows the 50-year contours for three wind classes obtained using Weibull and Log-normal distributions for   .The steepness limits using deep-and intermediate-water depth dispersion relation are also indicated.Using a 3-parameter Weibull model introduces a lower limit, the location parameter, capturing more adequately the lower tail of the   data within   class, being closer to the steepness curve, which can be observed directly from the scatter density plots.The resultant contours using the Log-normal distribution did not follow the data trend close to the wave steepness limit, resulting in relatively high   values above the steepness for a given peak period.This is seen for wind class 10-12 m/s, for   between 3-8 s.Considering that the natural period of large diameter monopile-based OWTs is in this range, a high dynamic response is expected even for moderate values of   , which can result in overestimating responses.
The difference between the steepness criteria is negligible for wave periods lower than 8 s and water depth of 30 m.However, the steepness limits differ considerably for peak periods larger than 8 s.The effect is pronounced for the contour in wind class 26-38 m/s.Here, it should be noted that the contour for wind class 26-38 m/s was used only for illustration purposes, and not for any response calculation.The use of such a large bin was due to insufficient data to obtain a finer resolution, which is a key limitation of the conditional modelling approach.Wave conditions over this range of wind speeds will vary considerably, and such large bins should therefore be avoided for response calculations.As it can be seen, using the Log-normal distribution resulted in a contour with   −  combinations, which, although reasonable for deep water depths, are notably above the steepness limit when the effect of the water depth is considered.This depicts the importance of evaluating the resultant contours based on the relevant steepness criteria, to avoid unrealistic   −   combinations.

Extreme responses along the contours
To investigate the variation of extreme responses along wind speed contours, two wind speeds were selected: the rated speed (10-12 m/s) and the cut-out speed (24-26 m/s).The response under consideration is the fore-aft bending moment at the mudline.The same soil and wave load models were used, specifically the NGI macro-element model and the 2 nd -order waves with the Morison equation using a constant value for   .TurbSim was employed to generate turbulent wind time series for the aeroelastic simulations.The mean wind speeds were set to the mid-values of the rated (11 m/s) and cut-out (25 m/s) wind classes.The turbulence intensity was defined based on the normal turbulence model (NTM) for IEC Class C.
Multiple sea states were analysed along the contours obtained from different statistical models.For each sea state, 20 one-hour time-domain simulations were conducted, considering random wind and wave seeds.The number of simulations was chosen to minimize  statistical uncertainty arising from seed variability, and was considered a reasonable compromise between accuracy and computational effort for assessing response variations along the contours.To mitigate the impact of wind and wave seed variability between different sea states, the same set of 20 seeds was employed for each load case analysed along the contours.It should be noted that using the same set of seeds keeps the same surface elevation profile for fixed   and different   , however it will give a different set of waves at different wave periods.Therefore, even using the same set along the contours the effect of seed variability is not eliminated.The influence of seed variability on the response is discussed in Section 5.3.3.
The short-term global maxima of the fore-aft bending moment at the mudline for each sea state were fitted using the Gumbel distribution (Bury, 1975).Eq. ( 23) shows the cumulative distribution function    () of the short-term global maxima, characterized by parameters   and   .These parameters are estimated using the method of moments, based on the expected value and standard deviation of the Gumbel distributed variable as in Eq. ( 24), where  ≈ 0.5772 is the Euler-Mascheroni constant.The representative extreme response for a sea state (ℎ  ,   ) along a 50-year contour for a wind class associated with a mean wind speed   at hub height was then found by solving Eq. ( 25), where   represents the percentile of the 1-hr extreme value distribution divided by 100.For most practical offshore problems dominated by wave loading, reasonable results are obtained 0.85 ≤   ≤ 0.95 (Baarholm et al., 2010;Haver et al., 1998;Kleiven and Haver, 2004) for 3-hr simulations, for 100-year return period.For the purposes of the study,    was obtained from the extreme value distribution for   = 0.9 using Eq. ( 26),

Rated speed (11 m/s)
The extreme mudline fore-aft bending moment responses at the rated speed were found to be equally important as those at the cut-out speed for both the 10 MW and 15 MW OWTs, due to the combination of wind loads induced by the maximum thrust force and moderate but relatively steep waves.The extreme responses at the rated speed were influenced by factors such as wave height, wave steepness, and the natural period of the structure.Figs.13(a) and 14(a) present the metocean contours at the rated speed for different conditional models of significant wave height (  ) for the 10 MW and 15 MW OWTs, respectively.The colour scale for the dots showing each state represents the 90 th percentile of the bending moment at the mudline.Additionally, Figs.13(b) and 14(b) display individual wave events extracted from the 20 one-hour simulations of the worst sea state identified for each OWT.The red line indicates the natural period of the OWTs.The wave height () is defined as the crest-to-trough height, and the period ( ) represents the time, between two consecutive zero up-crossings.The  (c-d) Wave events from 20 1-hr analyses for two sea states that resulted in the same  90% (10 MW).Black crosses indicate the events that resulted in the highest mudline bending moment.Red vertical line indicates the natural period of the 10 MW OWT.
wavelength of each wave event was calculated using Eq. ( 17), where the corresponding wave period was used instead of   .Subsequently, the wave steepness was calculated by the ratio of wave height to wavelength using  = ∕.
For the 10 MW OWT, the highest response was observed at the top of the contour.Along the contours,   gradually increases while the waves became less steep, and the wave energy content around the natural period (3.7 s) decreases.As depicted in Fig. 13(a) and in Table 2, the choice of conditional   model influenced the estimated extreme responses.The largest value of  90% was obtained using the 3parameter method of moments (MoM) Weibull model.Comparatively, the LonoWe model yielded a slightly lower (2.7%) value, while for the 2-parameter and 3-parameter maximum likelihood estimation (MLE) Weibull models, the same sea state (  = 2.66 m -  = 5.54 s) along the contours resulted in largest load  90% , which was 6.9% lower than the value obtained from MoM model.Despite the 10 MW wind turbine having a relatively short natural period, and significant aerodynamic damping due to the operational rotor, the effect of a steep sea state near the natural period was notable.For instance, wave events of the load case with   = 2.66 m and   = 5.54 s (Fig. 13(c)) resulted in similar extreme responses ( 90% ≈ 378 MNm -Fig.14) compared to wave events for sea states with significantly larger   but less steep and not in the vicinity of the natural period (Fig. 13(d)).The influence of wave steepness is also evident in the lower responses observed for sea states with similar   but longer periods.
A similar behaviour along the contours was observed for the 15 MW OWT, as shown in Fig. 14(a).However, the highest response was not found at the top of the contour.Due to its longer natural period (5.5 s), the 15 MW OWT is primarily affected by moderate but sufficiently large and steep wave events around the natural period, which induce high dynamic response, see e.g Fig. 14(b), rather than the largest waves located at the top of the contours.Similarly to the 10 MW  OWT, the choice of conditional   model influenced the estimated extreme responses.The largest extreme response was estimated using the 3-parameter MoM model, with  90% reaching 616 MNm.The LonoWe model resulted in a slightly lower (3.4%)extreme, while the 3-parameter MLE and 2-parameter Weibull models resulted in lower values of  90% by 7.3% and 8.7%, respectively.Wave steepness for some events exceeded theoretical limits and, in reality, wave breaking may occur; however, wave breaking was not considered in this study.

Cut-out speed (25 m/s)
At the cut-out speed, similar behaviour was observed for both the 10 MW and 15 MW OWTs, as depicted in Fig. 15.Extreme responses gradually increase along the contours, reaching the maximum values either at the top of the contour or for slightly smaller but steeper waves.For the 10 MW OWT, the maximum response ( 90% ≈ 457 MNm) was obtained for a sea state with   = 8.55 m and   = 14.94 s, which was 11.1% larger than the largest  90% observed at the rated speed.Similarly, for the 15 MW wind turbine, the load case with   = 8.81 m and   = 14.5 s resulted in the largest  90% ≈ 611 MNm, which was lower than the largest  90% at the rated speed (616 MNm).This suggests the significance of load cases at the rated speed for larger capacity wind turbines with increasing natural periods.The choice of conditional   model also influenced the estimated extreme responses at the cut-out speed as shown in Fig. 15 and Table 3.For the 10 MW turbine, the largest  90% was obtained using the LonoWe model.The 3-parameter model (MoM) resulted in a slightly lower (0.4%) extreme response, indicating negligible differences between the contours of the two models.On the other hand, the 3parameter (MLE) and 2-parameter Weibull models resulted in notably lower extreme responses, approximately 11.1% and 12.7% respectively, due to slightly less severe sea states for the two models.Similarly, for the 15 MW wind turbine, the LonoWe model resulted in the largest  90% ≈ 611 MNm, with the 3-parameter model (MoM) resulting in a slightly lower (2.2%) estimate.Meanwhile, the 3-parameter (MLE) and 2-parameter Weibull models resulted in similar  90% (8.5% lower) for two different sea states, respectively.

Effect of seed variability
The 90th percentile of the response for each sea state in Section 5.3 was based on the global maxima from 20 simulations with random wind and wave seeds.Although 20 simulations can provide a reasonable estimate of maximum response, there is still some uncertainty in quantile estimates as a function of number of samples due to seed variability.To evaluate the effect, the most severe load cases for the two OWTs for rated and cut-out wind speed were analysed with 100 random seeds.Then, using the sample size of 100 global maxima, a Gumbel distribution was fitted and scale and location parameters were estimated, as described in Section 5.3.Table 4 summarizes the load cases, the bending moment 90 th percentile  90% , and the corresponding Gumbel scale (  ) and location (  ) parameters.
The uncertainty of quantile estimation was evaluated using Monte Carlo simulations.Synthetic data samples of different sizes (between 3 and 100) were generated for a random variable  that is assumed to follow a Gumbel distribution, with location and scale parameters   = 0 and   = 1, respectively.Then, the 90 th quantile estimate for each set of the synthetic samples was determined.The procedure was repeated ten thousand times to create a distribution of quantile estimates, and the 2.5.%-97.5% confidence intervals were estimated.Then the confidence intervals of the distributed variable    for each load case were simply scaled following    =   +   ⋅ .Fig. 16 shows the results from the study.The top plots illustrate the 100 global maxima from the simulations, and the corresponding Gumbel fits.In the plots below, the blue/green markers for N = 1 show the global maxima from each simulation, and the red cross marker shows the  90% as estimated using all (N = 100) global maxima, which is used as the reference value.The shaded area indicates the 2.5% and 97.5% confidence intervals.As Fig. 16 shows, short-term maxima from 1-hr simulations have a notable spread.For both OWTs and for both load cases, the uncertainty gradually decreased as more global maxima were used to estimate  90% , which is expected.Higher variability was observed for the more severe sea states (LC2), at the cut-out speed for both turbines.This is also supported by the ratio between the scale (  ) and location (  ) parameters, which is related to the coefficient of variation (CV).The relatively lower CV for the load cases at rated speed, indicates that   is relatively smaller than   , leading to a narrower spread relatively to the mode of the distribution.It is worth noting that for both wind turbines, the  90% values at the rated and the cut-out speeds were comparable, with a difference ranging from 3.7% to 5.5%.Nevertheless, the results for the 15 MW OWT indicated the increasing importance of evaluating load cases at the rated speed for larger turbines, as it led to larger extreme values compared to the load case at cut-out speed.The results here are based on estimating the short-term extreme response distribution using only using the global maximum (single value) from each simulation.Other methods can be used such as peaks-over-threshold (POT) (Velarde et al., 2019) or the Naess-Gaidai ACER method (Naess and Gaidai, 2009), which consider a larger number of extremes from a given time series, resulting in better definition of distribution tails, and therefore smaller uncertainties.

Effect of soil modelling on extreme responses
The total damping in the system consists of aerodynamic, hydrodynamic, structural, and soil damping.The time-domain simulation tool inherently considers the effects of aerodynamic and hydrodynamic damping.To account for structural damping, stiffness-proportional Rayleigh damping was applied along the support structure, with a damping ratio () of approximately 1.0-1.1% of critical damping for the first natural frequency of each turbine.These values were based on published results (Shirzadeh et al., 2013) and typical values used for reference wind turbines (Bak et al., 2013;Gaertner et al., 2020).In the case of the  −  curves, soil damping was also incorporated using stiffness-proportional Rayleigh damping at four different levels.The amount of soil damping, expressed as a percentage of critical damping, was chosen to ensure that the total global damping ratio of the  −  curves closely matched the damping ratio of the macroelement model at the mudline, for low, intermediate, and high mean load levels as obtained from free-vibration tests.Fig. 17 illustrates the global damping ratio for the macro-element model and the  −  curves with different levels of soil damping applied in the analyses for each turbine.The presence of a non-zero slope in the global damping ratio for response amplitude reflects the hysteretic soil damping captured  by the macro-element model, which cannot be reproduced by Rayleigh damping.
To compare the macro-element model with the − curves, the most severe sea states for each wind turbine at rated and cut-out speeds were considered.As explained in the section on extreme statistical models (Section 5.3), twenty simulations were conducted for each sea state, and the 90 th percentile of the bending moment at the mudline was estimated.The effect of foundation modelling was evaluated for both the operational state, where aerodynamic damping dominates, and the parked state, where soil damping becomes the main damping contributor in the system.The same seeds were used for the foundation models and states to avoid stochastic variation in the extremes due to seed variability.Tables 5 and 6 summarize the extreme responses obtained from both the macro-element model and the  −  curves for different load cases, operational states, and damping levels.Additionally, Fig. 18 depicts the relative differences between the macro-element model (used as reference) and  −  models and for the various damping cases.
In operational states where aerodynamic damping dominates, the use of  −  curves led to only slight differences in extreme responses compared to the macro-element model.Furthermore, varying the soil damping in the  −  curves had a minor effect on the results.For the 10 MW turbine, the  −  curves overestimated  90% by approximately 1.6-2.6%for LC1 and 2.4-2.8% for LC2, when comparing the highest and lowest soil damping levels.Similarly, for the 15 MW wind turbine, the relative differences between the − curves and the macro-element model were around 1.0% for LC1 and below 0.5% for LC2.
In parked conditions, where the minimum values of soil damping in the − curves were used (corresponding to approximately 1.2% global G. Katsikogiannis et al. Fig. 18.Relative difference of  90% between  −  curves and macro-element (used as reference value) for the load cases shown in Tables 5 and 6. damping for both turbines), the extreme responses were overestimated.The overestimation ranged from 6.3% (LC2) to 10.4% (LC1) for the 10 MW turbine, and around 5.1% for both LCs for the 15 MW wind turbine.Increasing the amount of soil damping in the  −  curves gradually decreased the extreme responses for both load cases.For both turbines, a more significant effect was observed for LC1 (rated), as these load cases include wave events closer to the resonance period, making the responses more sensitive to damping variations.

Effect of diffraction on extreme responses
As mentioned in Section 4.3, diffraction effects can be important when the dimension of the structure is large compared to the wavelength, i.e.,  > 0.2.Therefore, evaluating their effect on extreme load and response estimates is important.To evaluate diffraction effects, the sea states for each OWT at rated and cut-out speed (see Section 5.3) were evaluated, with the turbines in operational state.For each sea state, twenty simulations were conducted, and the 90 th percentile of the mudline bending moment was compared.The load cases were analysed using the macro-element model and turbulent wind.For each 1-hour time-domain simulation with 2 nd order wave kinematics, the velocity, and acceleration at each time instant were extracted from mudline to instantaneous surface elevation , at nodal positions with dz = 0.5 m.Dynamic nodal forces were applied using the Morison equation with constant and frequency-dependent   (Section 4.3), using the extracted kinematics.The wave kinematics were assumed to be valid ± 0.25 m for each nodal position (half of the element below and above).The nodal positions along the monopile were compared to the instantaneous surface elevation at each instant, and the nodal forces were applied accordingly, considering if the node -corresponding element -is fully or partially submerged.A simplified illustration is shown in Fig. 19.
For a consistent comparison, and to avoid differences in results due to interpolation of wave kinematics within the software, the global  maxima using the nodal forces approach with constant and frequencydependent   were compared.Table 7 summarizes the  90% for the load cases analysed.For both OWTs, the load cases at the rated speed are mostly dominated by waves with shorter periods, where frequency-dependent   is lower than 2, resulting in lower extreme loads, and therefore in lower responses.For the 10 MW OWT,  90% was reduced by 1.8% while for the 15 MW OWT by 4.3%.The effect is more apparent for the 15 MW OWT G. Katsikogiannis et al. Fig. 20.Mudline response spectra from one simulation representative for LC1 (rated) and LC2 (cut-out) for 15 MW.Frequency-dependent (from MacCamy-Fuchs theory for diameter 11 m) inertia load coefficient is also plotted.due to the combined effect of its larger diameter, and the sea state with a lower peak period.For LC2 (cut-out), longer waves dominate extreme responses, where the theoretical   obtained from MacCamy-Fuchs is close to 2. Consequently, extreme responses were similar using both methods, with differences 0.5-1.6%, with frequency-dependent   model resulting in slightly larger response for the 10 MW OWT.
Fig. 20 shows an example of mudline response power spectral density (PSD) from one simulation for each of the load cases for the 15 MW OWT, using constant and frequency-dependent   .Together   is plotted over the frequency range.As shown, for the load case at rated speed, relatively short waves can result in large responses, where   values are expected to be lower than 2. As a result more notable effect is observed also in the response spectra between the two approaches.The same effect is also seen for the load case at the cut-out speed, nevertheless in that case lower frequency waves dominate, where   from MacCamy-Fuchs theory is similar to 2, and consequently diffraction is less important.

Conclusion
Estimating long-term extremes for OWTs involves uncertainties originating from different sources.This paper addressed some aspects of statistical and physical modelling uncertainties associated with extreme fore-aft bending moment responses at the mudline of monopile-based DTU 10 MW and IEA 15 MW reference offshore wind turbines.The environmental contour method using the global hierarchical approach and IFORM was used to establish contours for representative wind classes (rated and cut-out).Initially, the impact of different conditional probability distributions for significant wave height (  ) and peak period (  ), as well as methods for determining distribution parameters on the resulting contours were investigated.Then, the relative importance of using contours from different statistical models and seed variability on extreme responses were compared to the impact of different foundation and hydrodynamic load models, with a short summary shown in Table 8.
Several statistical models were employed for the conditional distribution of   given wind speed.Despite the challenge of capturing the tail behaviour when fitting the models to all observations (discussed in Section 2), for the site-specific dataset considered in the study, the 3-parameter Weibull model using the method of moments (MoM) to estimate distribution parameters, provided the best fit for   across all wind classes as it focuses on the upper tail of the data.The 3-parameter Weibull model using MLE resulted in less conservative contours compared to MoM model, particularly in moderate wind classes, leading to less severe sea states.The 2-parameter Weibull model yielded significant discrepancies compared to the observations, while the LonoWe model proved challenging to fit and was sensitive to the transition point between Log-normal and Weibull distributions.These effects were prominent in the upper tail of metocean contours, with differences between statistical models being more apparent in intermediate wind classes, attributed to the their inadequacy to capture the combined nature of moderate wind-sea states and a few severe sea states originating from swells.
For the conditional distribution of   given wind speed and significant wave height, the Log-normal and 3-parameter Weibull models were adopted.The Log-normal model deviated from the data trend near the wave steepness limit, particularly in intermediate wind classes, resulting in relatively high   values in the range of the turbines' natural periods.The conditional   Weibull model proved sensitive to the   class size, affecting   distributions and the resulting contours.Nonetheless, it provided a more accurate fit near the steepness limit, and compared to the commonly used Log-normal model, it can be considered as an alternative.It should be emphasized that other methods exist and this is the authors' recommendation given the two alternatives and the location considered in the study.To establish general applicability, more locations should be considered.The importance of using the appropriate steepness curve (deep or intermediate dispersion relation) is also evident to avoid unrealistic severe sea states.
The 3-parameter Weibull (MoM) and LonoWe models provided more reasonable contours, capturing the behaviour of the hindcast data in the upper tail for both wind classes, and resulting in higher extreme responses.in the largest  90% values between these were to range from approximately 0.5% to 3.5%, reflecting slight variations in the contours.Conversely, the and 3-parameter Weibull (MLE) models yielded the least conservative contours, with  90% values being approximately 7.0% to 12.7% lower, depending on the OWT and wind class.Furthermore, uncertainty in extreme responses was investigated, by analysing the effect of seed variability and the number of simulations on  90% values.The results demonstrated a higher variability on  90% for the sea states at cut-out speed, with variations up to approximately ±15% when using 20 1-hr simulations.For the larger turbine, it was found that evaluating load cases at the rated speed became increasingly important, as it resulted in larger  90% values.
A comparison was also made between a nonlinear elasto-plastic model (macro-element) incorporating damping and nonlinear stiffness during load reversals, and the conventional  −  curves.Variations up to 0.5% to 3.5% between the two models.LonoWe sensitive to transition .
Conditional distribution of   3-parameter Weibull/Log-Normal models were compared.3-parameter Weibull provided a more accurate fit near the steepness limit.
Seed variability/Number of simulations Significant variation in global maxima for each simulation and in quantile estimates  90% mainly at cut-out speed.

Foundation modelling
Nonlinear elasto-plastic model (macro-element) that inherently models hysteretic soil damping was compared to  −  curves.

Hydrodynamic modelling
Frequency-dependent and constant   were compared, using Morison and 2 nd order kinematics.
Up to 4.5% variation for the larger diameter OWT.
The representation of low, intermediate, and high levels of soil damping in the  −  curves employed stiffness-proportional Rayleigh damping.
The differences in extreme responses between the foundation models were within the range of approximately ±2% in operational states.However, in parked states, variations between −5% and 11% were observed, highlighting the importance of the foundation model in the absence of aerodynamic damping.Accounting for diffraction effects using frequency-dependent   becomes increasingly important for larger diameter and shorter sea states.Nevertheless, variations in estimated extremes compares to conventional Morison with constant   did not exceed 4.5%.
In conclusion, the estimated extremes can vary considerably based on the selected contour and stochastic variation due to seed variability.The foundation modelling was found to be important for parked conditions and sea states close to the natural period of the OWT, while accounting for diffraction effects became more crucial for larger diameters.Furthermore, the study emphasizes the increasing significance of load cases at the rated speed as wind turbine sizes and natural periods increase.The current study used 2 nd order wave kinematics with the Morison equation, without considering wave breaking.It is recommended that further studies, whether numerical or experimental, should be conducted, taking into account wave breaking, slamming loads, and additional factors such as wave directionality, wave spreading, and turbulence modelling in the incoming wind, to enhance the understanding of extreme response behaviour.Finally, in the present study, the percentile (90%) was chosen based on previous studies and established experience from other types of offshore structures.Further work including full long-term assessment is required to verify the appropriate percentile that corresponds to 50-year extreme responses of monopile-based OWTs for different wind classes and operational conditions, when the contour method is used.

Fig. 3 .
Fig. 3. 50-year contour sphere with radius  and contour circles of radius   for various wind classes  in Gaussian U-space.

Fig. 5 .
Fig. 5. Behaviour of  −  curves and macro-element models from the same decay test for a 15 MW OWT.

Fig. 6 .
Fig. 6.Comparison of inertia load coefficient for different monopile diameters at 30 m depth.Vertical lines represent the natural frequencies of the 15 MW and 10 MW OWTs, respectively.

Fig. 8 .
Fig. 8. 50-year contours using LonoWe model with various values of the transition point , between the Log-normal and Weibull models (15 MW OWT).

Fig. 9 .
Fig. 9. Hindcast data and fitted models for different   distributions (top) and resultant 50-year contours (bottom) for three wind classes (Left:10-12 m/s Middle:18-20 m/s Right:24-26 m/s) for 10 MW OWT.Colour in the hindcast data in the contour plots indicates the density of points and the horizontal lines the 10-and 50-year exceedance probability.

Fig. 10 .
Fig. 10.Wind-sea and swell components and contours for different wind classes (10 MW OWT).Colour in the hindcast data indicates the density of points.

Fig. 12 .
Fig. 12. 50-year contour for three wind classes using a Log-normal and a Weibull model for the   distribution.The wave steepness criterion is plotted.Colour in the hindcast data indicates the density of points.(Left:10-12 m/s Middle:24-26 m/s Right:26-38 m/s).

Fig. 13 .
Fig. 13.(a) Metocean contours at rated speed for different conditional   models for 10 MW.The colour scale of each dot (sea state) represents the 90 th percentile of the bending moment at mudline, obtained from 20 simulations.The same set of (20) wind and wave seeds were used for each sea state.(b) Individual wave events from 20 simulations for the most severe sea state found along the contour.Colorbar indicates the steepness of wave events.(c-d)Wave events from 20 1-hr analyses for two sea states that resulted in the same  90% (10 MW).Black crosses indicate the events that resulted in the highest mudline bending moment.Red vertical line indicates the natural period of the 10 MW OWT.

Fig. 14 .
Fig. 14.(a) Metocean contours at rated speed for different conditional   models for 15 MW.The colour scale of each dot (sea state) represents the 90 th percentile of the bending moment at mudline, obtained from 20 simulations.The same set of (20) wind and wave seeds were used for each sea state.(b) Individual wave events from 20 simulations for the most severe sea state found along the contour.Colorbar indicates the steepness of wave events.Black crosses indicate the events that resulted in the highest mudline bending moment.Red vertical line indicates the natural period of the 15 MW OWT.

Fig. 15 .
Fig.15.Metocean contours at cut-out speed for different conditional   models for the two OWTs.The colour scale of each dot (sea state) represents the 90 th percentile of the bending moment at mudline, obtained from 20 simulations.The same set of (20) wind and wave seeds were used for each sea state.

Fig.
Fig. Variation of extreme response due to seed variability.The analysed sea states are summarized in Table 4. Confidence intervals are scaled from Monte Carlo simulations of a synthetic Gumbel distribution.

Fig. 17 .
Fig. 17.Global damping ratio for  −  curves using different levels of soil damping (stiffness proportional), and for macro-element model.

Fig. 19 .
Fig.19.Illustration of nodal forces approach with dry, partially and fully submerged nodes with respect to instantaneous surface elevation .

Table 1
(Cheng et al., 2019) of the two monopile-based OWTs.Another reason is the uncertainty in parameter estimates due to finite sample size, particularly for cases with limited amount of data(Cheng et al., 2019).IEC(IEC, 2019)states that the joint environmental model should include the influence of possible upper limitations on   .

Table 2
Worst sea state along the contours for different conditional   models and corresponding  90% (rated speed).

Table 3
Sea state with largest load along the contours for different conditional   models and corresponding  90% (cut-out speed).

Table 4
Gumbel fit properties for the load cases analysed to evaluate seed variability.

Table 8
Summary of main observations in extreme responses from different statistical and physical load modelling aspects. 90% ranges Conditional distribution of   3-parameter Weibull (MoM)/LonoWe both give reasonable contours in the upper tail and higher extreme responses.