Critical assessment of hydrodynamic load models for a monopile structure in finite water depth

The response from steep and breaking waves for a monopile structure is investigated by analysis of experimental results and by application of numerical models for six irregular sea states. The experimental monopile was designed to reproduce the first and second natural frequencies of the NREL 5 MW reference monopile wind turbine. The measured response is reproduced in a finite element model using the Morison equation extended to full Lagrangian acceleration with second-order wave kinematics, and with fully non-linear kinematics and the axial divergence term. For the fully nonlinear wave kinematics, the additional point forces of Rainey and Kristiansen & Faltinsen (KF) are further added. In the latter model, the kinematics at the still water level are obtained by Taylor expansion of the kinematics from the free surface. The shear force at the sea bed and the structural accelerations are next compared between the force models and the experimental data. Among the findings are that the extreme force events are generally smallest for the second-order Morison approach, followed by the extended Morison model, and then the Rainey and KF model which produce similar results. While the total accelerations are found to generally match the measurements with fair accuracy, a modal decomposition shows that all models overpredict the response at the first eigenfrequency and underpredict it at the second eigenfrequency for extreme events. The latter is linked to the missing description of slamming loads in the modelling approach. The point force models of Rainey and KF are found to give quite similar results for the extreme events, the reasons for which are demonstrated by regular wave analysis.


Introduction
The number of bottom-fixed offshore wind turbines around the world has been increasing for the past decade and is expected to continue this trend in the years to come [1].In Europe, 81.7% of all substructures for offshore wind turbines are currently monopiles [2].Over their lifetime, many of these structures will encounter storms with steep and breaking waves which may produce large responses thus threatening the structural integrity.In particular, two load effects have been the focus of several research projects over the last years, namely ringing and slamming.
Ringing responses are an intermittent excitation of the first mode of the structure by a steep wave, not necessarily breaking, whose fundamental frequency lies far from the first eigenfrequency of the structure.Fig. 1 illustrates such an event: large first mode response is triggered by a very steep wave.The response builds up over one wave period and slowly decays.Ringing started gaining attention in the 1990s when it was observed during experimental campaigns on offshore oil and gas platforms [3], and has recently been observed in model tests of bottom-fixed offshore wind turbines as well [4][5][6][7].Although the mechanisms producing ringing responses are not yet fully understood, all research agrees on the fact that it is a nonlinear phenomenon [8][9][10][11], and that linear wave load models cannot reliably reproduce the complete ringing response [12][13][14][15].Slamming events occur when a breaking or near-breaking wave impacts the monopile, producing an impulse load on the structure.The duration of slamming loads is usually very small compared to the first eigenfrequency of the structure [17] and it has therefore been argued that slamming loads alone are less important for the first-mode response, since this is to a large extent driven by ringing loads [6].However, large responses around the second mode of the structure have been observed when slamming events occur, both in full-scale measurements [18] and during experiments [5,7].Suja-Thauvin et al. demonstrated that the dynamic amplification of the second mode of the structure can account for up to 20% of the maximum total response, thus showing that models that do not include slamming are likely to underestimate the response [16].
This discussion highlights the importance of having models for external forcing that can drive the dynamic response of at least the first and second modes when designing a bottom-fixed offshore wind turbine.In the North Sea, the assessment of the response of the structure under extreme environmental conditions -so-called Ultimate Limit State (ULS) analysis-is commonly performed following the DNV-OS-J101, DNV-RP-C205 and IEC-61400-3 standards [19][20][21].These standards suggest that the hydrodynamic loads produced by extreme waves are calculated with the Morison equation [22] applied with stream function theory wave kinematics [23].In cases where the wave is breaking, a slamming model must be applied.For plunging breaking waves, the method derived by Wienke and Oumeraci [24] is usually applied; whereas for spilling breakers, the method developed by Nestegård et al. [25] is generally used (an overview of the different types of breaking waves is given by Galvin [26]).
The standard hydrodynamic load models used to predict ringing and slamming were initially developed for oil and gas platforms in deep waters and present some limitations regarding application to offshore wind turbines.These models fail to accurately depict the physics, as they generally do not reproduce the balance between the measured response of the first mode of the structure and that of the second mode [27].In the present paper, four load models, briefly described hereafter, are analysed.Experimental data produced by DHI Denmark as part of the Wave Loads project [28] are compared to the responses modelled with those numerical models.
The first model uses wave kinematics that include components up to second order in terms of wave steepness (hereafter referred to as second-order wave kinematics).The hydrodynamic loads are then computed applying those kinematics to the well-known Morison equation [22].Compared to a traditional application of the Morison equation, the load model applied here includes convective acceleration terms.
The second model uses fully non-linear wave kinematics applied to the Morison equation and also includes the convective acceleration terms.In addition, a correction term presented by Manners and Rainey [29] is added to account for the fact that the cylinder is not slender in its axial direction.
The third model uses the same fully non-linear wave kinematics, and the force model developed by Kristiansen and Faltinsen [10], hereafter referred to as KF model.The point force of the KF model requires the kinematics at the mean water level.This is not trivial because when a wave trough passes the cylinder, the area around the mean water level is dry and no kinematics are available.A solution using Taylor expansion is suggested here.
The last model uses the same fully non-linear wave kinematics, and the so-called Rainey force model which is derived from energy balance arguments.When used with a circular cylinder, the model reduces to the Morison force model plus two additional terms.The first additional term is the previously mentioned axial divergence term [29].The second one is a point force at the surface which accounts for the change of kinetic energy in the fluid associated with the time variation of the wet portion of the cylinder [30].For Fig. 1.Illustration of a ringing event, recorded in the experiments presented by Suja-Thauvin et al. [16].In this figure the quasi-static contribution to the bending moment has been filtered out to keep only responses near the first eigenfrequency of the structure.

L. Suja-Thauvin et al.
long-crested waves, the Rainey and KF models are identical apart from the point force terms that act at the free surface elevation (Rainey) and at the mean water level (KF).The difference between these models is analysed in detail in Section 6.
The aim of the present paper is to compare the predicted responses in strong sea states at intermediate depths for the four models with experimental data.Special focus is given to first-and second-mode response, and the role of slamming.
The paper is organised as follows: Section 2 presents the experimental set up used by DHI, and Section 3 analyses the experimental data.In Section 4, we describe the numerical models in detail and some of the limitations in the way we apply them.Section 5 compares the responses obtained with the numerical models to those measured during the experimental campaign, and Section 6 analyses the differences between the Rainey model and the KF model.Conclusions are drawn in Section 7.

Experimental set-up
The experiments analysed in this paper were carried out at DHI Denmark in collaboration with the Technical University of Denmark (DTU) in the Wave Loads project [28].A more detailed description of the experiments is given in the report of the Wave Loads Project [28] and in the work by Bredmose et al. [7]; only the most relevant points are given here.The experiments were designed at a scale of 1:80 following Froude scaling.In the following, all dimensions are given at full scale with the model scale in parentheses and italic font, unless otherwise specified.For the tests presented here, one model was tested, which was a flexible circular cylinder of constant diameter and constant thickness of 6 m (0.075 m) and 0.144 m (1.8 mm), respectively, whose bottom was fixed onto a force transducer on a 1:25 slope.Note that the model does not extend below the seabed as is the case for full-scale wind turbines.Three water depths were tested, 20.8 m, 30.8 m and 40.8 m (0.26 m, 0.39 m and 0.51 m).Temporary walls and a rock berm were placed to create an enclosed area 17 m long and 18 m wide in the basin.The model was placed 7.75 m (model scale) away from the wave maker in the longitudinal direction.At the other end of the enclosed area, the rock berm was used to reduce wave reflection.
Two point masses of 937 t and 936 t (1.786 kg and 1.784 kg) were placed 128.6 m and 87.0 m (160.75 cm and 108.75 cm) above the sea bed, respectively, to match the first two fore-aft eigenfrequencies of the NREL 5 MW wind turbine [31], at f 1 ¼ 0:28 Hz and f 2 ¼ 2:0 Hz (2.5 Hz and 18 Hz), as shown in Table 1.Fig. 2 shows the dimensions of the experimental set up.The position ðx; yÞ ¼ ð0; 0Þ corresponds to the axis of the cylinder, and z ¼ 0 corresponds to the still water level.
The wave conditions correspond to situations where the turbine is expected to be idling, therefore, no aerodynamics were modelled.Wet decay tests were carried out and the damping for the first and second mode was found to be 1.7% and 2.7% of the critical damping, see Table 1.For the first mode, this is within the range of damping ratios measured on similar idling full-scale wind turbines (1.05%-2.8%depending on the wind speed and inclusion of tuned mass damper [32][33][34][35][36]).
Several wave gauges were placed in the wave basin at the centreline of the wave maker and towards the structure, and additional gauges were placed around the structure.One of the wave gauges was placed 53 diameters away from the cylinder in the transversal direction.This wave gauge is used in Section 4.1 to compute the second order wave kinematics.A three-component force transducer at the bottom of the structure and five accelerometers distributed along the cylinder measured the response of the structure.One of the accelerometers was placed at the level of the highest point mass and is used in Section 5.1.All data were sampled at 200 Hz (model scale) giving a time step of 0.005 s (0.0447 s in full scale).

Environmental conditions
Two different irregular sea states based on the JONSWAP spectrum [37] were considered.Combined with the three tested water depths, this gave six wave conditions, which are referred to hereafter as load cases (see Table 2).Some of these load cases have been used by Bredmose et al. [7] and Robertson et al. [38].In the following, each load case is defined by its targeted (as opposed to measured) significant wave height H S , targeted spectral peak period T p and water depth h.All sea states have a target spectral peakedness factor of γ ¼ 3:3.Each load case is named with a letter corresponding to the depth (S-shallowest, M-medium and D-deepest) and a number corresponding to H S .Note that the waves producing the largest responses of each case correspond to waves in hydrodynamically intermediate water [39].
For each sea condition, we calculate an irregular Ursell parameter using the method suggested by Stansberg [40]: where k p is the spectral peak wave number, calculated with the spectral peak period and the second order dispersion relation given in equation (1) of Kirby and Dalrymple [41].The Ursell parameter is typically defined for regular waves, and it has been showed that for U r > 0:33, second-order wave theory is not valid to represent the regular wave kinematics [42].Stansberg [40] extended the analysis to irregular waves and showed that this limit is also applicable.As shown in Table 2, load cases S8.3 and S11, which correspond to the shallowest water depth, are above this limit.
The maximum height H b before wave breaking is expected is also included in Table 2.In load cases S8.3 and S11, this wave height was reached.The breaking wave height is calculated following the equation suggested in the standards [19,20]: All sea states were run for 3 h (full scale) and produced about 1000 waves.During the experiments and especially for the most extreme cases, several waves broke at the cylinder.

Analysis of the measured data
We now analyse the data obtained during the experimental campaign.An analysis of the measured shear force at the seabed and acceleration at the highest point mass has also been performed by Bredmose et al. [7].A positive shear force corresponds to an external load in the direction of the waves (i.e.positive x-direction).Fig. 3. Exceedance probabilities for the measured wave heights.

Statistical analysis
Fig. 3 shows the exceedance probability plots for the measured wave heights for all six load cases.The wave height is measured between two zero-downcrossing values of the free surface elevation time series.As expected, the load cases with H s ¼ 8:3 m show lower wave heights than the ones with H s ¼ 11 m.For the H s ¼ 8:3 m load cases, non-linear shoaling effects induce an increase in wave heights for decreasing water depth.For the H s ¼ 11 m load cases, this effect occurs for the main part of the population.For the most extreme events however, the shallowest water case (S11) produces lower wave heights than the other two cases.This is due to wave breaking, which occurs more frequently in shallow waters.
Fig. 4 shows the exceedance probability plots for the seabed shear force.To obtain Fig. 4, two consecutive zero-downcrossings in the free surface elevation signal are found, and the maximum positive value in the force for each wave is identified.This procedure is applied throughout the paper for all exceedance probability plots.The measured shear force was filtered to remove frequencies above 4 Hz (full scale) by applying a Butterworth filter, which removes most of the contribution of the higher modes of the structure.The third eigenfrequency of the structure lies at 5.6 Hz [28].Excitation of modes higher than second was observed for breaking wave events, but this response is not representative of a full-scale turbine as these modes were not tuned to fit any full-scale turbine.
As for the wave heights, larger H s implies larger shear forces for all three depths.Considering the H s ¼ 11 m cases, for 80% of the population (P exceedance > 0:2), larger water depth implies larger shear forces.This is due to i) a larger water column acting on the structure and ii) a larger mode shape displacement further away from the sea bed.For deeper water, a load at the surface will induce higher modal loads than the same load at the surface for a shallower case.At P exceedance � 0:12, the shear force in S11 becomes larger than that of D11 and M11.The shallower water depth produces more non-linear (steeper) waves which then have the necessary frequency content to excite the first mode of the structure.At P exceedance � 0:03, the same phenomenon happens for M11, which then shows larger shear forces than D11.Similar trends can be observed when comparing load cases corresponding to H s ¼ 8:3 m.

Spectral analysis
Fig. 5 and Fig. 6 show the power spectral density of the free surface elevation for ðH S ; T P Þ ¼ ð8:3 m; 12:6 sÞ and ð11 m; 14 sÞ, respectively.For the milder sea state of H s ¼ 8:3 m shown in Fig. 5, the energy content around the peak frequency f p is very similar for all three load cases.Differences appear in the difference-frequency region (between 0 and 0.05 Hz) and the sum-frequency region (around 0.16 Hz), where the load cases with shallower water have higher energy content.This is consistent with the fact that waves in shallow waters are more non-linear than in deep waters.
For the extreme sea state of H s ¼ 11 m shown in Fig. 6, load cases in shallower water contain less energy around f p .This is due to wave breaking, which is induced by the smaller depth and shoaling on the slope towards the structure.Also, energy is transferred from the spectral peak into sub-and super harmonics.As for the previous sea state, the shallower cases show more energy in the differenceand sum-frequency region.Fig. 7 and Fig. 8 show the Fourier components of the seabed shear force of the structure for ðH S ; T P Þ ¼ ð8:3 m; 12:6 sÞ and ð11 m; 14 Fig. 9. Exceedance probability plots for the wave height for all load cases.Comparison between the measured free surface elevation and the free surface elevations obtained with the FNL and second-order kinematics models.
sÞ, respectively.From these figures it appears that deeper waters induce larger shear forces around f p .As explained in Section 3.1, this is due to a larger water column and to the loads acting at a level with a larger mode shape displacement.Around 2f p , the shear force in shallower waters becomes larger, due to the waves being more non-linear and therefore containing more energy at this frequency (as shown in Figs. 5 and 6).
The dynamic amplification of the first mode (0.28 Hz) is clearly visible in both Figs.7 and 8.For the H s ¼ 11 m cases (Fig. 8), the amplitude of the shear force around the first eigenfrequency is similar for all depths: the relative contribution of the response at the first natural frequency increases with decreasing water depth.This is consistent with Fig. 6, as shallower depths induce steeper waves that contain energy at higher frequencies.This effect is not visible in Fig. 7, where the amplitude of the shear force around the first eigenfrequency for the deepest case (D8.3) is larger than for the other cases.A possible explanation is that for the H s ¼ 8:3 m and T p ¼ 12:6 s cases, the sea state is milder and therefore produces less steep waves.Except for the most extreme events (as shown in Fig. 4), the effect of non-linearity of the waves is less important for the mild sea states than for the stronger sea states of the H s ¼ 11 m and T p ¼ 14 s cases.

Numerical models
In Section 5, we reproduce the measured shear forces with different numerical models.This is done in three steps: first the wave kinematics are reproduced, then the hydrodynamic loads are computed and finally the shear force of the structure at the seabed is calculated.
The models used to perform these three steps are presented in this section.

Wave kinematics
In this study, we compare models using two different sets of wave kinematics.The first set of wave kinematics includes components up to second order in terms of wave steepness.This set of kinematics is obtained by first linearizing the free surface elevation measured by the wave gauge placed at the side of the cylinder, see Section 2.1.The linearization is carried out by first removing the differencefrequency wave components, and then iteratively selecting the cut-off frequency of a low-pass filter such that the reconstructed wave spectrum (including first and second order terms) gives the best possible match to the spectrum of the measured free surface elevation.Further details of the linearization method are given by Bachynski et al. [4].The reconstructed set of kinematics is hence referred to as 'second-order'.
The second set of wave kinematics is produced with the fully non-linear potential flow solver OceanWave3D, presented by Engsig-Karup et al. [43].The code discretizes the domain by means of finite differences, and integrates the Laplace equation with fully nonlinear boundary conditions at the water surface.The Laplace problem is solved in an ðx; y; σÞ coordinate system, where Here η is the free surface elevation, h is the water depth and ðx; y; z; tÞ are the spatial and temporal coordinates.Using the nondimensional σ coordinate, the Laplace equation is solved in a time-invariant grid, increasing the computational efficiency of the solution.The kinematics in the physical coordinate system are calculated from the potential Φðx; y; σÞ via a chain rule derivation.For example, the horizontal particle velocity is In the computations, the waves were not induced via a wave paddle generator as in the experiments, but were rather enforced in a so-called relaxation zone.First, a linear representation of the incident waves was obtained from the experiments by the analysis of four wave gauges, located further than 76 cylinder diameters from the model in the offshore part of the domain, by applying a variant of the reflection analysis of Goda and Suzuki [44] as presented by Bredmose et al. [45].This linear wave field was next imposed in a 640.0 m (8.0 m) long flat-bedded relaxation zone, located between the domain offshore boundary and the beginning of the slope.The approach is similar to the one used in Bredmose et al. [7] and other OceanWave3D studies.In the flat-bedded relaxation zone, the numerical solution from the wave model is in each time step replaced by a weighted average of the numerical solution and the prescribed linear wave field.The weighting function undergoes a smooth transition from the prescribed wave field to the numerical wave field through the relaxation zone.Beyond the relaxation zone, only the internally-calculated potential is retained.In the experiments, the distance between the domain boundary and the slope start was only 80 m (1 m), due to the different wave generation technique (see Fig. 2).
The waves are then allowed to propagate through the shoaling numerical domain fulfilling the Laplace equation for the velocity potential and the fully non-linear kinematic and dynamic free surface boundary conditions.To represent the effect of wave breaking, which is beyond the physics of a potential-flow model, an ad-hoc dissipative filter was applied in each time step, whenever the downward local vertical particle acceleration exceeded 0.5 times gravity.This method is well suited for spilling type of breakers and the cutoff at 0.5 g, as suggested for example by Longuet-Higgins [46], has provided good results in earlier studies, e.g.Schløer et al. [47].Other research suggest that different cutoff values can provide accurate results (for example Dawson et al. [48] recommend 0.3 g), however, recent work by Pierella et al. [49] with the same data set as this paper shows that 0.5 g gives the best agreement across the whole range of investigated sea states.
The kinematics are then sampled in pre-defined positions in space and time and saved to memory.The fully non-linear kinematics are referred to as 'FNL' in the following.
Fig. 9 shows the exceedance probability plots of the wave heights.For each load case, the exceedance probability was computed for the measured and modelled waves.The FNL wave heights generally agree with the measured wave heights.For the most nonlinear sea state (S11), the 2% largest wave heights are over-predicted by the FNL model.This is expected to be due to breaking in the experiment, which is not accurately reproduced by the 0.5 value of the ad-hoc breaking filter.As expected, although the second-order wave field is produced from a wave gauge at the same position as the structure, the kinematics from the second-order model produced wave heights slightly lower than those from the FNL model but still in good agreement with the measurements.As noted in Table 2, S8.3 and S11 have an Ursell number higher than 0.33 and therefore exceed the limit of validity of second-order theory.The second-order wave kinematic model is therefore not studied for these two cases.Note that the second-order kinematics are obtained assuming a flat bottom, which is not consistent with the test conditions.

Load models
Four load models are used to compute the hydrodynamic loads on the structure: i) the Morison equation presented by Morison [22] with second-order wave kinematics, hereafter 'Morison second-order', ii) the Morison equation with FNL kinematics, hereafter 'Morison FNL', iii) the model proposed by Kristiansen and Faltinsen [10], with FNL kinematics, hereafter 'KF' and iv) the model proposed by Rainey [50], with FNL kinematics, hereafter 'Rainey'.
Note that the motion of the structure was not taken into consideration for any of the models.It is assumed that the motion of the structure was small enough for relative viscous damping to be insignificant.This assumption was verified by modelling all load cases with and without the motion of the structure, where differences of less than 1% were found.The main damping of the response was due to structural damping.

Morison second-order
The first model is based on the well-known Morison equation [22], often used due to its simplicity.The difficulty of using the Morison equation resides in the selection of the wave kinematics and coefficients.As mentioned in the introduction, the simplest solution of applying linear wave kinematics does not produce accurate results in steep waves.The Morison equation is composed of a drag term and an inertia term with and where η is the incident wave elevation, ρ is the water density, a is the cylinder radius, u is the horizontal wave particle velocity and C D and C M are the drag and inertia coefficients, respectively.In the present analysis, the values of C D and C M are obtained by calculating an irregular version of the KC and Reynolds numbers [51] based on the measured T P and the standard deviation of the particle velocity from the FNL kinematics.Table 4.11 in Sumer and Fredsøe [52] is then used to determine C D and C M .In the present study, for all load cases, C D ¼ 1:1 and C M ¼ 2 were chosen.The derivative of the particle velocity is taken as the Langrangian derivative in 2D: where w is the vertical particle velocity.The total load can thus be expressed as In the Morison second-order model, we use the second-order kinematics to compute the loads.The loads are thus integrated up to the second-order wave elevation.

Morison FNL
The Morison FNL model is the same as the Morison second-order model, except that 1) the FNL kinematics are used and, 2) to account for the fact that the cylinder is not slender in its axial direction, the term ðC M À 1Þρπa 2 u∂w=∂z is added [29].As in the previous model, the Lagrangian acceleration is applied.The total load can then be expressed as: As for the Morison second-order model, values of C D ¼ 1:1 and C M ¼ 2 are used.

KF model
The first version of the so-called FNV (short for Faltinsen-Newman-Vinje) model was developed by Faltinsen et al. [9] in order to predict ringing responses of oil and gas platforms observed on model tests in the 90s [3].This model was initially developed for deep water regular waves based on perturbation theory and includes the effect of long wave diffraction and the waves scattered by the cylinder.It is consistent up to third order in terms of wave steepness.Krokstad et al. [53] presented a validation of the model in deep water, and Paulsen et al. [54] showed that the third-order load matched calculations computed with CFD.This model was further extended to irregular waves by Newman [55] and eventually to finite water depth by Kristiansen and Faltinsen [10].The latter implementation, referred to as KF model, is the one used in the present analysis, and the load is calculated as follows F'ðz; tÞ is a distributed load given by where m a ¼ ρπa 2 is the added mass in surge.F ψ is the force due to the scattered potential and is applied at z ¼ 0.
where g is the acceleration due to gravity.In the original paper from Kristiansen and Faltinsen, the kinematics used to calculate F ψ are to be taken at z ¼ 0 [10].However, in the current analysis, no FNL kinematics are available in the dry parts of the structure: whenever a trough passes the structure, there are no kinematics at z ¼ 0. Therefore, the kinematics at z ¼ η are Taylor expanded to z ¼ 0. This approach is assessed in section 4.4.1.
In the present implementation, we add a drag term based on the Morison equation with C D ¼ 1:1.The total force calculated with the KF model is then with C M ¼ 2.
In the present analysis we use the FNL kinematics to compute the KF loads.To be consistent to third order, contributions higher than third order should be removed from the hydrodynamic loads thus calculated.This is not trivial and was not carried out in this analysis.This implies, for example, that load components of fourth order are included due to the multiplication of contributions of second-order kinematics in the convective acceleration, whereas fourth order terms due to the scattered potential will be missing since the scattered potential is consistent to third order only.The effects of this shortcoming are analysed in section 4.4.2.

Rainey model
The Rainey model is derived from energy balance arguments.When applied to circular cylinders, it reduces to the Morison equation as presented in section 4.2.1 with two additions.The first one accounts for the fact that the cylinder is not slender in its axial direction by adding the term m a u ∂w ∂z [29].The second one is a point force at the free surface to account for the change in fluid kinetic energy due to the variation in time of the submerged portion of the cylinder (equation 7.4 in the original paper [30]).This force is given by where η is the instantaneous free surface elevation at the centre of the column.The Rainey model is then

Table 3
Summary of the load models.'X' means that the term is included in the model, '-' means that it is not.All models include a drag and an inertia term, and use the Lagrangian definition of the acceleration.
In the present implementation, the FNL kinematics are used for the Rainey model, C D ¼ 1:1 and C M ¼ 2. It should be noted that the integrated force of the Rainey model and the KF model in the present implementation are the same and equal to the force computed in the Morison with FNL kinematics model.The only difference between the Rainey and KF models lies in the point force and is further explored in section 6.
Table 3 is a summary of the four load models implemented in this analysis.Note that for all models, the term in the integral is vertically distributed along the cylinder whereas loads F ψ and F η are applied at one point.

Structural model
The hydrodynamic loads calculated in the previous section are used as input to a finite element representation of the cylinder described in section 2.1 with 160 elements.The finite element software applied here is Ashes, which uses the Newmark-Beta integration method to solve for the deflections of the structure [56].The added mass coefficient used is C a ¼ C M À 1 ¼ 1, and the first and second eigenfrequencies obtained with the finite element model match perfectly those given in section 2.1.The damping applied is Rayleigh damping, with the damping ratios tuned to match that of the first and second modes of the physical model.Fig. 10 shows the mode shapes of the structure obtained with Ashes for h ¼ 30:8 m.The mode shapes and the eigenfrequencies show slight variations depending on the water depth.These variations are accounted for in the numerical simulations.

Limitations of the current method to compute the F ψ term 4.4.1. Selection of the kinematics
As explained in section 4.2.3, the theory developed by Kristiansen and Faltinsen implies that F ψ must be computed with kinematics at z ¼ 0 and applied at z ¼ 0 [10].In the current work, since no FNL kinematics are available at z ¼ 0 when η < 0 at the structure, the kinematics at z ¼ 0 are obtained through first-order Taylor expansion: uðz; tÞj z¼0 ¼ uðz; tÞj z¼η À η ∂uðz; tÞ ∂z L. Suja-Thauvin et al.
∂uðz; tÞ ∂t In the present study, ∂ 2 u=∂z∂t was not output by the fully non-linear solver, so it was calculated by numerical differentiation using the available grid points.
To assess the impact of this method for calculating F ψ , the response shear force of the structure is calculated in regular waves with kinematics calculated with the stream function wave theory, following the method presented by Rienecker and Fenton [23], which solves for a velocity potential in the whole domain and therefore provides kinematics in the dry parts of the domain also.Three procedures are studied: � z ¼ 0 based F ψ : the kinematics are taken at z ¼ 0, as proposed in the original method � Taylor expansion based F ψ : the kinematics are obtained from first order Taylor expansion (equation ( 17) and ( 18)).� z ¼ η based F ψ : the kinematics are taken at z ¼ η and the force is applied at z ¼ η To assess the differences among these formulations, the zero-downcrossing characteristics of the wave that produced the largest shear force of the S11 case are selected: the wave height is H ¼ 13:4 m and the wave period is T ¼ 15:2 s.The water depth for this case is h ¼ 20:8 m.As recommended in the standard IEC-61400-3 [21] for a wave of these characteristics, the order of the stream function is taken as N ¼ 11.Fig. 11 shows the shear force as calculated with the three procedures.The approach with the kinematics taken at z ¼ η produces an overestimation of the response of about 25%.Using the Taylor expansion approach produces an underestimation of the shear force of 3%.This procedure is therefore kept in the rest of the analysis.
Adding more terms to the Taylor expansion produced good results in the simulations with the regular waves but gave clearly unphysical loads when applied to irregular waves.Therefore, it was decided to keep the first-order Taylor expansion.

Contributions of higher than third order
In our implementation, we consider all orders in the incoming wave kinematics, as explained in section 4.2.3.To assess this limitation, the incoming kinematics must be separated into contributions to different orders.Since it is not straightforward to perform such a separation on a stream function wave, a Stokes fifth order wave is used instead, implemented following the work by Fenton [57].The Stokes expansion does not converge for the regular wave previously analysed in section 4.4.1, therefore a larger depth is considered.The zero-downcrossing characteristics of the wave producing the largest response of the M11 case are selected: the wave height is H ¼ 13:7 m and the wave period is T ¼ 11:2 s.The water depth for this case is h ¼ 30:8 m.Two procedures are analysed here: � KF3: only kinematics that produce loads up to third order are kept, both for the integrated load and for the F ψ term, as proposed in the original method.� KF5: all kinematics are kept.The integrated load is integrated to the instantaneous fifth-order free surface.Fig. 12 shows the shear force of the structure for these two procedures.The KF5 procedure produces a maximum shear force about 15% lower than the KF3 procedure.Investigating the kinematics, it is found that higher-order kinematics reduce the excitation load at lower harmonics.For instance, the fourth-order kinematics reduce the second harmonic of the load, and the fifth-order kinematics have a similar effect on the third harmonic of the load.As a consequence, the excitation load obtained with the KF5 procedure is lower than that obtained with the KF3 procedure.This is consistent with the work of Kristiansen and Faltinsen, who showed that the KF5 Fig. 11.Shear force of the structure in regular waves with the KF model.F ψ is computed following three different procedures.

L. Suja-Thauvin et al.
procedure generally produces lower loads than the KF3 procedure at the first, second and third harmonic of the force [10].They also point out that the KF3 procedure tends to overestimate the excitation loads when compared to measurements, and finally show that using the KF5 procedure, while still overpredicting with respect to the measurements, gives a closer match [10].
In the case of fully nonlinear kinematics, it is most natural to apply the full kinematics with no attempt of decomposition into orders.While the above analysis is restricted to regular waves and fifth order, the results support that this approach is applicable.This will be further tested in Section 5 where comparisons with measurements are presented.

Numerical responses
In this section, we assess how the different numerical models perform by comparison to the measurements.As pointed out in section 4.1, for S8.3 and S11, the Ursell parameter shown in Table 2 is above the classical limit of 0.33, above which second-order theory is not sufficient to predict the wave kinematics [40].The shear forces produced by the Morison second-order model therefore largely underpredicted the measurements and are not shown in this section.

Statistical analysis
Fig. 13 shows the exceedance probability plots for the shear force at the sea bed for the six load cases and four load models.The numerical results are compared against the experimental data.
Below shear forces of about 3.5 MN, all models compare well with the measurements, with a slight underprediction for the 8.3 cases and a slight overprediction for S11.For extreme forces, the models that use FNL kinematics overpredict the D8.3 and D11 cases and underpredict the M-S cases.The Morison second-order model generally produces smaller results than the other models, except for the M8.3 case.Generally, the Morison second-order model produces the smallest forces, followed by the Morison FNL model and then the Rainey and the KF models, which produce similar results.
For load cases D-M-S8.3 and D-S11, and with the exception of the largest force in M8.3, the forces calculated by the KF and the Rainey model lie within 12% of the measured forces.However, none of these two models can be said to generally over-or underpredict the results.For M11, all models underpredict the extreme peaks.This analysis was also performed with a less strict wave breaking filter and produced better results for M11.The sensitivity of the analysis to the wave breaking filter was not studied in further detail.
The Morison FNL model is identical to the Rainey and KF models without their respective point forces.It can be seen in Fig. 13 that this model underpredicts most of the extreme events, which shows the importance of the point loads.The KF model has been shown to overpredict the third harmonic loads for steep waves [10].In this analysis, the steepness for the waves producing the largest shear forces is calculated with the following formula: where k is the wave number of each wave calculated from the second-order dispersion relationship (eq (1) given by Kirby and Dalrymple [41]) based on the downcrossing period and H is the crest-to-trough height.The steepnesses thus calculated are well above the 1=40 limit for which Kristiansen and Faltinsen showed that the third harmonic load is overpredicted by the KF model [10].Despite this overprediction, the modelled extreme forces agree reasonably well with the measurements.
To examine the dynamic load effects in detail, the structural accelerations are now analysed.Fig. 14 shows the exceedance probability plots for the absolute value of the acceleration at the highest point mass of the structure.The absolute value of the acceleration is used because the largest acceleration in a steep wave load event with dynamic response can be in both the positive and the negative direction.For the measured acceleration, Bredmose et al. found that increasing depth led to an increased acceleration level due to the larger moment arm.Extreme accelerations, though, were strongest at the smallest depths due to the impulsive loads from very steep and breaking waves [7].
In general, the models using the FNL kinematics produce similar results, with Morison FNL producing lower accelerations than the KF and Rainey models.
With the exception of D11 and D8.3, all models match the main part of the population fairly well.The effect of the Ursell number is illustrated by the differences between the Morison second-order model and the models using FNL kinematics: in the D8.3 case, the nonlinearities in the waves are not as significant as in the other cases (as indicated by the low Ursell number) and therefore the models using FNL kinematics do not produce larger accelerations than the Morison second-order model.As the Ursell number increases, the accelerations produced by the models with FNL kinematics become larger than those produced by the Morison second-order model.
For the deep water cases, the load models quite well match the experimental data in D8.3 while D11 shows an overprediction.For the M and S cases, the load models produce results close to the measurements for the main population (about P exceendance > 3%) while To understand this effect, the accelerations are decomposed around the first and second eigenfrequencies of the structure.The decomposition is performed by applying a Butterworth band-pass filter.For the first mode, the filter was of ninth order and a frequency band of f 1 � 0:022 Hz was used.For the second mode, the filter was of sixth order and a frequency band of f 2 � 0:2 Hz was used.Different frequency bands and filter orders were tested.Although the amplitude of the filtered acceleration showed a slight sensitivity to those parameters, the qualitative results did not exhibit significant variations.Note that the accelerations that did not occur around the first or second eigenfrequency are not considered in this analysis.This implies that the sum of the first and second mode accelerations is not equal to the total acceleration.Fig. 15 illustrates such a decomposition for the largest measured acceleration of M8.3.Fig. 16 shows the exceedance probability plots for the structural acceleration filtered around the first eigenfrequency.Note that the order of the events shown in this figure might be different than that of the unfiltered acceleration shown in Fig. 14, as the largest unfiltered acceleration does not necessarily exhibit the largest first mode acceleration.
For the main population of events, the results follow a similar trend as the unfiltered acceleration shown in Fig. 14.This is consistent since non-extreme events are mainly the results of first mode forcing.Except for D8.3 and a few extreme events in M-S 8.3, all the models overpredict the first-mode acceleration.As for the shear force and total acceleration, the accelerations are generally smallest for the Morison FNL model while Rainey and KF produce similar results.The second order Morison results overpredict largely for D8.3 and M8.3.For M11 it gives the smallest accelerations and thus closest to the measurements, while for D11 it is generally in between the Rainey and KF results.It is not included for the shallow cases of S8.3 and S11.Fig. 17 shows the acceleration around the second mode of the structure.The load models driven by fully non-linear kinematics, provide remarkably similar results.Relative to these results, the second-order Morison model provides lower accelerations.For all load cases, the models largely underpredict extreme events.This is due to the absence of a slamming model in the numerical models: large events are produced by large waves which in the shallow water cases are breaking and thus will excite the second mode of the structure [18].de Ridder et al. [5] showed based on similar experimental data that triggering the second mode of the structure has a large influence on the acceleration at the top.Suja-Thauvin et al. [27] compared different numerical models to calculate loads on offshore wind turbines and demonstrated that slamming models were required to produce impulse loads with a duration short enough to trigger second mode response.Since none of the models analysed in this study include a slamming model, the underestimation of the second mode acceleration is expected.
The decomposition of the acceleration around the first and second eigenfrequencies of the structure shows that the models tend to overpredict the response around the first mode but underpredict the response around the second mode.

Single event analysis
We now analyse two characteristic events of the experimental campaign, marked with an 'X' in Figs. 4 and 13.For each event, four time series are shown comparing the measurements and the numerical models: � The free surface elevation.Note that the KF, the Rainey and the Morison FNL models use the same kinematics and therefore show the same free surface elevation.� The modal excitation, computed with the mode shape of the first mode following the equation where F ex is the hydrodynamic load calculated with the models described in section 4.2 and ψ 1 is the first mode shape.Note that no measurement of the excitation load is available.
� The shear force at the seabed.� The acceleration at the highest point mass.Fig. 18 illustrates an event where a mild shear force was recorded.No strong dynamic amplification appears in the measured or modelled shear forces.The lowest plot shows significant acceleration around the first mode and a lower acceleration at the second mode.This plot also shows that the Morison second-order model produces higher acceleration at the first natural frequency than the other models, which is consistent with Fig. 16.This is due to the different kinematics: the upper graph of Fig. 18 shows that the wave obtained with the second-order kinematics is steeper than the one obtained with the fully nonlinear kinematics.This produces a larger amount of high frequency content in the load, which leads to a larger excitation of the structure's first mode.
It can also be noted that the KF, the Rainey and the Morison FNL models produce very similar excitation loads and therefore very similar shear forces.As the waves composing the studied event have low steepness, the particle velocities and acceleration are also low.The F ψ and F η terms are of third order in terms of wave steepness, so these point forces become small compared to the F I term of the Morison equation defined in equation (5).Therefore, the excitation loads from the KF and the Rainey models are mainly due to the In Section 5.1 it was found that, only the models that use the FNL kinematics can correctly predict the most extreme events.To understand this, the event that produced the second largest measured response of M11 is now analysed.Fig. 19 shows the same analysis as in the previous example for this event.The upper plot shows that the wave modelled using the second-order kinematics model is less steep than the measured one.Combined with the absence of a point load, the lower steepness produces lower excitation loads than the other models and the Morison second-order model therefore underpredicts both the shear force and the acceleration, consistent with the statistical analysis of section 5.1.
The steepness of the measured wave and the large second mode energy observed both in the measured shear force and acceleration indicate that the wave is breaking.The wave associated with the FNL kinematics is significantly steeper than the one associated with the second-order kinematics but less steep than the measured one.The shear force and accelerations graphs show that the models using the FNL kinematics produce first mode response but fail to trigger the second mode.This suggests that the excitation forces calculated with these models do not have the required frequency content to significantly excite the second mode of the structure.
The KF and Rainey models generate similar modal forces, which translates into both models giving similar shear forces, as observed in Fig. 13.This is consistent with Section 5.1, where it was found that the KF and the Rainey models produced similar responses.
Note that for this event, there is a mismatch in the phasing of the numerical and the measured free surface elevation.This results in a similar mismatch in the shear force and acceleration and illustrates the need for an accurate wave profile when attempting to match the response at a deterministic level.13.This figure also shows that the Rainey model and the KF model produce similar results.However, it is known that F ψ ¼ 8F η in the small amplitude limit [8,9,58].
To understand how the two models can produce similar results despite this apparent eightfold ratio, the point forces are analysed in detail in the following by comparing F ψ and F η for regular waves.The waves analysed are shown in Fig. 20: the three water depths studied in this work (h ¼ 20:8 m, h ¼ 30:8 m and h ¼ 40:8 m) are tested, a period T ¼ 11:2 s is chosen for all waves and the wave height H is varied to go from small steepness waves up to the limiting wave height for regular waves.As in section 4.4.2, the stream function theory is used to model the kinematics of the waves.The order of the stream function is chosen according to the guidance lines included in Fig. 20.
We focus on the third harmonic of the F ψ and the F η forces, since this load component is important for ringing events.The third harmonic is extracted by band-pass filtering the time series for the two forces.Fig. 21 shows the amplitude of the third harmonic of the two point forces for the waves shown in Fig. 20.In this figure, for low wave height, the amplitude of the third harmonic predicted by the F ψ term is about eight times larger than that of the F η term.Note that this ratio is not visible in the total force (for example in Fig. 13) because the contribution of both point forces is not significant compared to the integrated load for low wave heights.
As the wave height increases, the ratios between amplitudes of the third harmonic decrease until, for the largest wave heights, the third harmonic of both point forces becomes similar.The result at small amplitudes is consistent with the previous research mentioned above [8,9,58]: the present analysis considers all orders, but at small amplitudes the terms beyond third order are not significant compared to third-order terms and F ψ � 8F η .For larger amplitudes, however, the higher orders play an important role and the third-harmonic load amplitudes become close as the waves approach the limiting wave height.This means that the apparent under-estimation of the third-harmonic load by the Rainey model disappears for large waves, where this force component is important.This can explain why the results of the Rainey and KF model in Section 5 have generally been found to be similar, despite the known difference for F ψ and and F η in the small-amplitude limit.
For small amplitude waves, Rainey [58] shows that the eightfold difference between F ψ and F η comes from neglecting the diffraction of the waves by the cylinder and suggests the addition of a so-called 'surface distortion' force to compensate for the difference.This force is given by With this addition, the point force in the Rainey model produces similar loads to the F ψ point force in the small amplitude limit.However, the derivation of the surface distortion force F dist assumes that the Froude number Fr ¼ c= ffi ffi ffi ffi ffi ag p , where c is the wave celerity, is small [58].This is not the case for waves producing large responses of the structure, and given the results of the present analysis, the missing forcing is compensated by higher-order contributions for large waves.The analysis performed in this section indicates that the eightfold difference between F ψ and F η present in the small amplitude limit Fig. 18.From top to bottom: free surface elevation, modal excitation, shear force and acceleration for an event recorded during D8.3.In the second subplot, the differences in modal force F m1 for the Morison FNL, the Rainey and the KF models are indistinguishable.
disappears for large waves because the difference between the wave kinematics taken at z ¼ 0 for F ψ and at z ¼ η for F η becomes significant.This also justifies the procedure of Taylor expanding the kinematics from z ¼ η around 0 to apply them to the KF model.We next analyse the difference between computation of F ψ at z ¼ 0 (by Taylor expansion from z ¼ η) and direct computation at z ¼ η which is formally consistent to leading order.The ratio of the two approaches is shown in Fig. 22, which presents the ratio of the third harmonic of the term u 2 du=dt taken at z ¼ 0 and taken at z ¼ η.For large waves, the kinematics at z ¼ η are significantly larger than at z ¼ 0. This implies that a direct calculation at z ¼ η of this third-order forcing term will give very large third-order forces.This approach was tested for irregular waves and led to strongly over-predicted extreme accelerations (not shown).The present analysis explains the reason for this and further suggests that the best way of using the KF model with fully nonlinear kinematics is to compute F ψ at z ¼ 0, in terms of a Taylor expansion of the kinematics from the free surface.

Conclusions
Experimental data of an idling bottom-fixed offshore wind turbine under extreme irregular wave conditions have been analysed.Two point masses were placed on the scaled model so that its first and second natural frequencies fit those of the NREL 5 MW turbine mounted on a monopile [31,59].Three water depths (20.8 m, 30.8 m and 40.8 m) and two JONSWAP spectra (H s ¼ 8:4 m, T p ¼ 12:6 s and H s ¼ 11 m, T p ¼ 14 s) were considered, for a total of six load cases.The shear force at the seabed was measured and it was found that for the main population of events, the largest shear forces correspond to the deeper water cases.This is due to a larger water column acting on the cylinder and the loads acting at a point with a larger modal displacement.For the largest events however, the Fig. 19.From top to bottom: free surface elevation, modal excitation, shear force and acceleration for the event that produced the second largest shear force in M11.Note that the high frequency oscillation of the third and fourth plots correspond to the second eigenfrequency of the structure.trend was inverted, with the largest shear forces occurring for the shallower depths.In shallower water, the waves are more non-linear and carry energy at frequencies close to the eigenfrequencies of the structure and are breaking to a larger extent.This led to larger extreme shear forces at the smallest depth.
In the second part of the paper, four numerical models were compared against the experimental results: the classical Morison equation with second-order wave kinematics, the Morison equation with fully non-linear kinematics, the Rainey model presented by Rainey [30] and the KF method presented by Kristiansen and Faltinsen [10].The last three models use fully non-linear wave kinematics as input.Exceedance probability plots showed that the wave heights obtained with the fully non-linear wave model are generally in good agreement with the measurements.The KF model requires using kinematics at z ¼ 0, which are not available for the FNL kinematics when a trough passes the cylinder.Instead, the kinematics at z ¼ η were Taylor expanded to z ¼ 0.
The Morison equation with second-order kinematics matched the main population of events well but underestimated the largest measured shear forces.There are two reasons for this underprediction: first, second-order kinematics miss some higher harmonics in the loads that are necessary to excite the first eigenmode of the structure.Second, this model did not include a point force load.
The models using fully non-linear kinematics predicted the most extreme shear forces reasonably well, but they were not in general conservative.Generally, the smallest forces were obtained with the Morison FNL model.The forces obtained from the Rainey and KF (in Taylor expanded version) models were generally very close.Similar findings were achieved for the acceleration at the highest point mass close to the top of the structure, where again the smallest accelerations were obtained with the Morison FNL model, followed by the Rainey and KF models.These three force models generally matched the main population of measured acceleration reasonably well, with underprediction, however, of the extreme accelerations.
To understand where the inaccuracies of these models come from, the acceleration was decomposed into contributions around its first and second mode eigenfrequencies.It was found that the models generally overpredicted the first mode acceleration but underpredicted the second mode acceleration for large response events.This was confirmed by the analysis of the time series for two characteristic events, where it was seen that the load models could produce ringing responses but did not produce significant response at the second eigenfrequency.
The underprediction of the second mode shear force by all numerical models was expected since no slamming load model is implemented in the current analysis.As shown by Suja-Thauvin et al. [16], second mode response is triggered by the impulse load produced by waves breaking at the structure.Recommended future work is to include a slamming model consistent with the numerical models presented in this paper.
In the last part of this paper, the excitation from the point forces of the Rainey model, F η and of the KF model, F ψ were compared for regular waves of different steepnesses.In the small amplitude limit, F ψ is eight times larger that F η .This difference reduced to unity for large waves, mainly because the difference in the wave kinematics (taken at z ¼ 0 for F ψ and at z ¼ η for F η ), which is negligible in the small amplitude limit, becomes significant.
Several limitations to this work should be pointed out: the analysis was based on data obtained through an experimental campaign at a 1:80 scale, with all the limitations and uncertainties inherent to this type of testing.The study also dealt with only one pair of values for first and second eigenfrequencies, corresponding to a 5 MW wind turbine in idling conditions.Current wind turbines have exceeded this rated power, and it is expected that the second eigenfrequency of larger turbines will be lower than in this study.This will change the relative contribution of the first and second modes in the response and could allow for mechanisms other than wave breaking to trigger second mode response.A limitation inherent to the hydrodynamic load models used in this analysis is that none of them includes fourth-and higher-order loads due to the hydrodynamic-structure interaction.These could influence the overall excitation loads on the structure.In addition, it was observed that the response was dependent on the wave breaking filter applied to the fully non-linear kinematics, which shows the need to deal with breaking wave kinematics accurately.Finally, this analysis did not include responses at frequencies higher than the second eigenfrequency.The influence of such higher frequencies could have a relevant impact on the overall response of the turbine.
It should also be noted that the inputs to the load models applied were kinematics based on measured free surface elevation time series.The difficulties inherent in producing accurate kinematics play a role in how well the modelled responses match the measured ones.However, for designing substructures for offshore wind turbines in the industry, the common workflow implies generating random sea states rather than using measurements.This removes the need for kinematics that accurately fit measurements and demands that the models perform well on a statistical rather than a deterministic level.Therefore, the models based on fully non-linear wave kinematics analysed in this paper have the potential to be used in the design of bottom-fixed offshore wind turbines, and can be further improved by adding a consistent slamming model.Finally, to be relevant for the industry, a remaining challenge will be to reduce the computational cost of producing fully non-linear wave kinematics.A database of fully non-linear kinematics is under production within the DeRisk project [60].

Fig. 2 .
Fig. 2. Experimental set-up (not to scale).Model scale dimensions are shown in italic font.

Fig. 5 .
Fig. 5. PSD of the free surface elevation for the H s ¼ 8:3 m cases.

Fig. 4 .
Fig. 4. Exceedance probabilities for the measured shear force.The two events analysed in section 5.2 are marked with an X.

Fig. 6 .
Fig. 6.PSD of the free surface elevation for the H s ¼ 11 m cases.

Fig. 7 .
Fig. 7. Discrete Fourier transform of the shear force for the H s ¼ 8:3 m cases.

Fig. 8 .
Fig. 8. Discrete Fourier transform of the shear force for the H s ¼ 11 m cases.

Fig. 12 .
Fig. 12. Shear force for a fifth order Stokes waves with H ¼ 13.7 m, T ¼ 11.2 s in h ¼ 30.8 m depth.

Fig. 13 .
Fig. 13.Exceedance probability plots for the shear force.The two individual events analysed in section 5.2 (one in D8.3 and the other in M11) are marked with a black 'X'.(For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Fig. 14 .
Fig. 14.Exceedance probability plots for the absolute acceleration at the highest point mass.

Fig. 15 .
Fig. 15.Decomposition of the measured acceleration at the highest point mass for the largest acceleration during M8.3.

Fig. 17 .
Fig. 17.Exceedance probability plots for the absolute acceleration at the second eigenfrequency.

Fig. 20 .Fig. 21 .
Fig. 20.Selected waves for the comparison between F ψ and F η (taken from the IEC-61400-3 standard [21]).The water depth in this figure is denoted d

Table 1
Eigenfrequencies and damping ratios of the first two modes of the model.

Table 2
Characteristics of the analysed load cases.The load cases are named after the depth (Sshallowest, Mmedium or D -deepest) and the significant wave height (8.3 m, 11 m).