Harmonic structure of wave loads on a surface piercing column in directionally spread and unidirectional random seas

We investigate the harmonic structure of wave-induced loads on vertical cylinders with dimensions comparable to those used to support offshore wind turbines. Many offshore wind turbines are designed, so the structural resonance is two or three times the fundamental wave loading period, which may be excited by higher harmonics of wave loading. The suitability of a ‘Stokes-type’ model for force is examined. We analyse experimental data for unidirectional and directionally spread sea-states. We demonstrate that approximate harmonic components may be extracted from a random time series, using a novel signal processing method. The extracted harmonics are shown to follow a ‘Stokes-like’ model, where the higher harmonics in frequency are proportional to the n-th power of the linear inline force component and its Hilbert transform. Results show that the third harmonic component of force is fitted less well by this model, which is consistent with the literature. A key new result is that the harmonic coefficients for spread and unidirectional seas are nearly identical when fitted as powers of the linear inline force and its Hilbert transform. This finding has the potential to greatly simplify the process of generating harmonic data for new wavelength to cylinder and wavelength to depth ratios. This also implies that the size of the inline component of the nth\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n\textrm{th}$$\end{document} harmonic for force in a directional sea relative to the same component in a unidirectional sea scales as cosnσθ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\cos ^n \sigma _\theta $$\end{document}, where σθ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma _\theta $$\end{document} is the rms directional spreading angle of the incident wave field. Hence, higher harmonics will be of less importance in directionally spread seas than unidirectional ones. We show that the computationally fast harmonic decomposition approach taken here can reproduce the shape and magnitude of the loading from non-breaking waves waves in a wide range of realistic unidirectional and directionally spread sea-states. The proposed force model has potential as an engineering tool for computationally fast prediction of nonlinear wave loads.


Introduction
Offshore wind turbine monopiles are exposed to the hostile maritime environment and must be designed to survive Ultimate Limit State and Fatigue Limit State Loads. Nonlinear waves and their interaction with a structure results in forces which are more nonlinear than the waves themselves.

B Mark L. McAllister
Wave-induced forces may be capable of inducing a dynamic response which may be detrimental to structural and foundation integrity if not properly accounted for in the design process, as a dynamic response can increase the effects of fatigue and accumulated foundation displacement (Richards et al. 2020). Monopiles for offshore wind turbines are typically designed, so that the structural natural frequency is around two or three times the frequency of severe storm waves (Kallehave et al. 2015) and in between the rotor-and blade-passing frequencies of the turbine. Any loading close to this natural frequency, such as from the higher harmonics of wave forces, may contribute significantly to the structural fatigue and the ultimate limit state loads, as analysed in Schløer et al. (2016) and Wang et al. (2021), respectively.
Various theoretical models based on perturbation expansions in wave steepness have been developed in an effort to understand the physical mechanisms and quantify the nonlinear components of wave loads. Most of the diffraction literature is limited to second order in wave steepness (see, e.g., Eatock Taylor and Hung 1987;Molin 1979;Newman 1996), and is based on classical Stokes-type expansions. Beyond this, there have been some efforts to extend diffraction theory, e.g., Malenica and Molin (1995) who follow the classical wave steepness assumption, and the FNV model (Faltinsen et al. 1995) (see also Kristiansen and Faltinsen (2017) for compact cylinders). Numerical methods have also been used, although the challenge of resolving loads associated with short components of the wave field (requiring high wavenumbers for their description) remains significant. There are many published experimental works on higher harmonic forces in the literature, e.g., Krokstad and Stansberg (1995); Huseby and Grue (2000); Riise et al. (2018), which report parametric analysis of higher harmonics from experimental tests. The origin of such harmonics was investigated in Grue and Huseby (2002) and Chaplin et al. (1997), with particular focus on ringing effects. Grue et al. (1993) report a secondary loading cycle that occurs after the wave crest has passed the centre of cylinder. The timing of this loading is around one-quarter to one third of the local wave period after the main peak of the force.
Recent work on wavegroups has used experimental and numerical methods to examine higher order harmonics (Chen et al. 2018(Chen et al. , 2019Feng et al. 2020;Fitzgerald et al. 2014) and fit these to a 'Stokes-like' model where the time-history of the total force, F(t), can be written at least conceptually as the sum of harmonics of some underlying linear force time-history F 1 cos(φ), which scale as F n 1 F(t) = F 1 cos(φ) + F 2 1 Γ 2 cos(2φ + ψ 2 ) + F 3 1 Γ 3 cos(3φ + ψ 3 ) +F 4 1 Γ 4 cos(4φ + ψ 4 ) + O(F 5 1 ), where Γ n is the nth harmonic force coefficient, and φ = ωt + ψ 1 is the phase angle. Unlike the harmonic components in a conventional Stokes' expansion for surface waves (e.g., Fenton 1985) which are all in phase, ψ n is the phase shift associated with each harmonic force component. Here, we only consider superharmonic terms and do not include components such third-order frequency difference terms, that scale as the cube of the wave amplitude but have a frequency close to linear. We also omit sub-harmonic terms which occur at lower frequencies than linear force time-history (we discuss these briefly in §4.4). This paper has two primary objectives. The first is to extract and explore the harmonic structure of wave-induced loading of a surface piercing column for a range of representative sea-states, using a novel signal processing method. Extracting harmonic components of measured wave-induced loads has not previously been achievable for broadband irregular sea-states. The second is to illustrate a method that may be used to generate synthetic loading data for offshore wind turbine design. Previous work has primarily involved uni-directional waves and sea-states with realistic amounts of directional spreading have not been considered. As such, there remains an open question as to how the higher harmonics of the loading in a directionally spread sea relate to those in unidirectional waves and indeed whether the hypothesis that the load can be separated into harmonics is still valid, and if so how can this be done. The present paper addresses this issue as well as reporting results from an extensive set of laboratory measurements of forces on a vertical cylinder in random seas from the DeRisk project ).

Experimental setup
The data we use herein are taken from an extensive experimental campaign undertaken at DHI in Denmark, with the aim of better understanding the nature of extreme wave loading on offshore wind turbine foundations in directionally spread and unidirectional seas (see Bredmose et al. 2016). The experiments were performed in the shallow water basin at DHI Hørsholm. The results of the experiments are Froude scaled to the field using a geometric length scaling ratio of 1:50. The full-scale description is representative of wind turbine structures in the North Sea. Hereafter, all values are presented at field scale unless stated otherwise.
The radius of the cylinder was 3.5 m (0.07 m at lab scale), which is representative of an offshore wind turbine monopile support structure. In all the tests, the diameter of the monopile is sufficiently large, such that the wave-induced forces are dominated by inertial loading. Two water depths were considered for the unidirectional wave tests, 33 and 20 m (0.66 and 0.4 m at lab scale). We only present directionally spread experiments performed with 33 m water depth. The vertical cylinder was installed a distance of 365 m (7.3 m at lab scale) from the wave paddles (see Fig. 1). The uniform cylindrical form extends throughout the entire water depth.
The monopile is connected to force transducers, two at the top and two at the bottom, measuring inline (F x ) and lateral forces (F y ). The sampling frequency for the force measurements was 51.56 Hz (256 Hz at lab scale). In the present paper, we consider inline forces, i.e., those in the mean wave direction (positive x-direction). The monopile structure comprised a near-rigid simply supported cylinder and was intended to exhibit a static response to the incident waves. With this in mind, the support structure was designed, so that the natural frequency was high enough to ensure that only slamming loads generated any significant dynamic structural response, with the more slowly varying load components inducing quasi-static behaviour. Post-processing of the force measurements found the natural frequency to be ≈ 3 − 4 Hz (33 and 20 m water depth, respectively). We point out that these dynamic properties were not intended to be consistent with those of a full-scale turbine structure.  Each case corresponds to a particular wave spectrum with H s the significant wave height, T p is the peak period and γ the peak enhancement factor of the JONSWAP spectrum, and duration of 6 and 70 h for unidirectional and directionally spread experiments, respectively. Also, a indicates the cylinder radius and d indicates the water depth The raw force signal is low-pass filtered at ∼ 7 to 8 times the frequency 1/T p (so, a high-frequency cut-off of 0.76 Hz at full-scale), for all the cases, to remove noise when searching for force peaks. We will subsequently refer to this filtered signal as the measured signal throughout the paper.

Analysis methods
Our aim is to investigate the harmonic structure of waveinduced loads, and examine the suitability of Stokes' type force model, for a surface piercing column in irregular waves. This process involves four steps, the methods of which we explain in the following section. First, we explain how we extract the average harmonic components of extreme loading from the random force time-histories (Sect. 3.1); second, how we fit a Stokes' type model to the extracted harmonic components (Sect. 3.2); third, how we estimate the added mass coefficient using undisturbed surface elevation and force measurements (Sect. 3.3); fourth, how we nondimensionalize force, so that the fitted model may be used to predict wave loading (Sect. 3.4).

Extraction of harmonic components from random signal
Previous studies, which have examined focused wave groups, have used 'phase separation' to extract the harmonic components of wave loading. The 'four phase' separation approach is outlined in Fitzgerald et al. (2014), and the method is introduced in terms of a Stokes expansion for surface waves and then applied to hydrodynamic forces. This method requires controlled phase shifted repeats of the same experiment (φ, φ + π/2, φ + π , φ + π 3/2) which are combined to extract harmonics. This approach works extremely well for wave groups (Zhao et al. 2017), and random sea-states (Kristoffersen et al. 2021;Pierella et al. 2021) but has been shown to be less reliable in certain conditions (Adcock et al. 2019), owing to wave breaking and reflections moving components out of alignment. In the present work, we use the technique described in Zhao et al. (2018a, b) to generate average time series of extreme events with different phase shifts, for random wave experiments which were not repeated with four phase shifts. Using the same approach that is applied to phase shifted repeated experiments in Fitzgerald et al. (2014), these average time series can be combined and filtered to extract average harmonic components F n . The black lines in Fig. 2 illustrate the first five harmonic components extracted using this process for Exp. 4. These extracted force time-histories represent the average nonlinear structure of the individual harmonics contributing to extreme wave loading events. Full details of the steps involved in this process including illustrative examples are given in Appendix A.

Fitting of Stokes' type force model
We aim is to fit a Stokes' type force model, based on the extracted average linear force time-history, to the average harmonics extracted from the irregular force times series. Through appropriate normalisation, such a model may then be used to include higher harmonic force components in the first-order reliability models to predict wave loading distributions. Our model is based on the assumption that the harmonic components of wave loading follow the form shown in (1) and each term may be expressed as Γ n F n 1 cos(nφ − ψ n ), where F 1 is the envelope of the extracted average linear extreme force time-history. Γ n a scaling coefficient and ψ n is a phase coefficient. Unlike for Stokes terms for surface waves, higher harmonic forces exhibit a phase shift relative to the fundamental linear harmonic. Therefore, both the amplitude, Γ n , of the harmonic, and the phase shift relative to the fundamental, ψ n , need to be estimated for each term F n in our model. In practice, we fit cos(nφ) and sin(nφ) and then convert to amplitude Γ n and phase ψ n where α n and β n are unknown coefficients determined using a least-squares approach. When fitting each term, the in-phase cos(nφ) and out-of-phase sin(nφ) force component of each harmonic may be approximated using the linear signal and its Hilbert transform making use of the multiple angle formulae. For example, the second harmonic may be expressed as and, assuming we have a narrow-banded signal, F 1 sin(φ) may be approximated using the Hilbert transform of the linear signal F H 1 , and F 1 cos(φ) ≈ F 1 which gives from which the coefficients α 2 and β 2 may be fitted based on information contained within the linear component alone. The same process is used for each harmonic to determine the values of α n and β n that produce the best model fit F n of the extracted harmonic component F n .

Estimating inertial load coefficient C M
Having defined a model for the higher harmonics of the inline force in terms of the linear component of force, we now turn to the relationship between the incident waves and the linear force component. For unidirectional waves, the first-order inertial wave loading may be expressed as where the coefficient C M combines the Froude-Krylov and added mass forces, ρ is the water density, g is the acceleration due to gravity, A is the linear wave amplitude, a is the cylinder radius, k the wave number, and d is the water depth. For directionally spread waves, it is common practice is to reduce unidirectional force using an inline kinematics factor, cos(σ θ ). We can estimate the value of C M using the spectra of the inline force and the quasi-undisturbed wave elevation at WG4 over the linear frequency range (see Fig. 1) measured in Exp1 -11. This process yielded values ranging from 1.95 to 2.3. This is slightly larger than values presented using the same data in Pierella et al. (2020), as they also include the small contribution from drag in their estimation of C M . The larger values we estimated corresponded to less steep waves and for experiments carried out at shallower water depths. For the majority of waves, the values we estimate lie around 2 (close to the theoretical value for spatially uniform flow around a compact cylinder).

Force non-dimensionalisation
The n-th harmonic component of force F n may be nondimensionalized as This is effectively the scaling factor of Huseby and Grue (2000) with an additional factor of 2π . Using (6), and assuming a value of C M = 2, we may estimate the wave amplitude associated with large linear wave loads (unidirectional deep water waves) where max(F 1 ) is the maximum of the averaged profile of F 1 , which is the same as the maximum of the envelope of the linear component max( F 2 1 + F 2 1H ). On substitution of the expression for the peak linear amplitude A (9), the nondimensional force becomes or which may be rewritten in a more general form as where C F is a normalisation coefficient defined as results in a normalized force that accounts for the effects of water depth and directional spreading of the incident waves.

Quality of coefficient extraction
Before considering values of the fitted force coefficients, we evaluate the quality of the fits achieved for the different experiments. Disagreement between the extracted harmonic components and the model fits to describe these may arise if the approximate extraction method we present fails to accurately decompose the signal, or if the assumptions we make about the harmonic structure of force signal are incorrect. There is also the issue of statistical variability, whether enough individual examples are included in the data. As an example, Fig. 2 plots the the extracted harmonic components and the harmonic model based on fitted coefficients for Exp4. Overall, there is good agreement between the profile of the extracted harmonics and the model fit achieved based on only the linear force component.
To quantify the quality of the fitted force model, we use the measure R 2 which is calculated as where F is the extracted profile, F the mean, and F is the fitted model. A value of 1 indicates a perfect fit. Table 2 presents the R 2 values of the harmonic fits for the all experiments. For each harmonic, the R 2 values were computed over a range where the envelope is more than 5% of its maximum. This is done to exclude the effects of noise at times far away from the centre of the linear force envelope in time.
On average, the second, fourth, and fifth harmonic fits have R 2 values over 0.9, indicating that the assumed force structure in time for these harmonics fits the experimental data well. In contrast, the third harmonic on average has lower R 2 values. This discrepancy can also be seen in Fig. 2. Intriguingly, for some cases, the third harmonic does appear to fit better than others. For the deeper water cases, the largest mismatch is for Exp1 and Exp2 of the unidirectional cases and Exp102 and Exp103 of the directionally spread cases. These with the fitsF n based on the linear force F n 1 using the ordinary leastsquares method (red dashed lines) for Exp4. The grey shade area represents the envelope of the extracted harmonics, and the red shaded area represents the envelope of the linear signal raised to the appropriate power and scaled to the peak of the group (color figure online) cases have similar sea-state parameters. However, we have no suggestions as to why these particular sea-states would be different to the rest. For the directionally spread experiments, which are longer in duration, the average R 2 values for the third harmonic are higher than the corresponding unidirec-tional experiments. Therefore, some of the issues in fitting F 3 may be overcome when there are more large events to average over.A more systematic variation of spectral parameters may also help reveal if any specific physical processes are causing the issues we observe. The discrepancies we observe may suggest that what is extracted as the third harmonic component in these particular sea-states is not completely bound to the fundamental harmonic, indicating the possibility of additional phenomena to the Stokes-like harmonic components. The third harmonic may be influenced by a local physical process with its own dynamical time-scale.
Inconsistencies in the higher harmonics have previously been attributed to various mechanisms which occur on a different time-scale to the linear waves. In particular, the work of Grue et al. (1993) stresses the contribution of a secondary load cycle to higher harmonic forces, while the study of Swan and Sheikh (2015) ascribes the nonlinearities to the generation of two types of scattered waves resulting from the interaction with the cylindrical structure. Both works note the occurrence of fluid circulating around the cylinder and giving rise to nonlinear effects at an additional time-scale.
A recent paper by Ghadirian and Bredmose (2019) presents analysis of the interaction of a focused wave group with a vertical cylinder. They performed simulations in Open-FOAM and concluded that the secondary load cycle is caused by suction when the piled-up water at the back of the cylinder falls under gravity. This accelerates downwards as the local wave elevation drops after the peak of the crest has passed. Paulsen et al. (2014) found that the harmonic content of the secondary load cycle inline force is around the sixth harmonic of the incident wave. It is therefore not expected to contribute to classical ringing loads in the third harmonic range, or to any of the 2nd-5th harmonics that we analyse here.
In our analysis, we do not observe a significant secondary load cycle in the averaged profile of the extreme wave loadings. Such effects may be lost in the averaging process, based on the timing of the peak of the linear signal, used to extract the harmonic components. However, the secondary cycles are present in the average of the forces following the top wave crests, i.e., on the averaged force profiles conditioned on maximum crests. The observation suggests that the secondary cycles are primarily linked to the interaction of large wave crests with the structure, and they do not generally lead to the largest peaks in the force time-histories. We note that the passage of wave crests past the column is associated with the zero-crossing of the linear inertia load rather than its peak. The averaging used in this work is based on the properties of 100s of force peaks. If the 2nd load cycle is not phase locked to these peaks, then the averaging process will remove it. This is a possible explanation as to why it appears in individual force time histories, and in the results of Schløer et al. (2017) where only a few examples were averaged, but not in our analysis.
Our analysis assumes that the harmonic structure of the loading has a form consistent with the sum harmonics of a Stokes expansion. Therefore, the n-th harmonic in frequency scales as the n-th power of the linear force component (and through this with the n-th power of the wave amplitude). One might also expect to find load components arising from drag. At the Keulegan-Carpenter number and Reynolds number regime for an offshore wind turbine column, the load linear in wave frequency is dominated by the inertial (potential-flow) contribution. However, we would expect there to be some small drag contribution. Assuming the Morison u|u|-form, there will also be Fourier harmonics within this drag force at n-odd multiples of the wave frequencies. It is possible that the relatively poor performance of our Stokes-type model for the third harmonic may be due to a drag harmonic. We would not expect any significant drag harmonics for n=2 or 4, and we note that for n=5 our model fits the experimental results rather well. Therefore, our analysis is used to build a force representation under the assumption that drag loading is small.

Magnitude and phase of harmonic coefficients
We now turn to the magnitude of the fitted coefficients. Table 3 presents these for the different experiments carried out in this study. There is clearly consistency between dif- ferent experiments and some identifiable trends in the data. The harmonic coefficients, particularly for higher harmonics, are significantly larger in the shallower water depth (Exp6 to Exp11). Whilst there is some variation in the phase data, this still is reasonably consistent across all the data. For example, the phase of the second harmonic is consistently close to the theoretical result for a compact vertical cylinder in deep water of π/2 relative to the crest of the linear force component. The magnitudes of the harmonics in unidirectional and spread seas cases for the same water depth match very well; this aspect of the results is examined in more detail below.

Comparison of unidirectional and directionally spread results (d = 33 m)
The duration of each of the unidirectional experiments is too short for us to look robustly at how estimates of the different harmonics varied with the steepness of the corresponding waves. However, for the very long (roughly 25000 waves) directionally spread wave experiments, it is possible to examine the fitted coefficients in greater detail. We first compare the coefficients in the same probability of exceedance range as in the 2D tests; so, the top 100 (largest to 100th largest) peak forces approximately correspond to the 10th to 999th largest cases for the 3D seas. These are the data shown in Table 3. Next, we investigate the variability of the coefficients over different ranges of exceedance probability. Figure 3 depicts the variation of the amplitude and phase coefficients, respectively, for all six of the directionally spread seas. Each of the points in the plot corresponds to the average of 500 peak forces. The variation of the coefficients are made versus k p A where A is the average of the amplitude for blocks of 500 waves, sorted in the same way as the peak forces. This is done by first sorting the amplitude of forces and waves in decreasing order. The variables are then grouped into separate bands, each comprising of 500 cases (i.e., 1-500, 501-1000,....).
There are 20 bands for each sea-state in the plots, and this required the consideration of the top 10000 cases of force and elevation amplitudes. The coefficients are almost consistent over a wide range of parameters. The largest variations occur in the third harmonic, which as discussed above, we suspect has components not directly bound to the fundamental mode, and in the fifth harmonic at lower exceedance probabilities where its magnitude (Γ 5 ) itself is small and close to the level of background noise in the measured force time history. Figure 3 demonstrates good agreement between the coefficients found in unidirectional and directional sea-states with otherwise similar parameters. The colors in this figure denotes pairs of directional spread and unidirectional sea-states with otherwise comparable properties. This suggests that unidirectional tests may be used to make an engineering estimate of the harmonics of wave loading in directionally spread seastates, by including the effect of wave spreading into the force coefficients. This potentially greatly simplifies the process of generating harmonic data for new wavelength to cylinder and wavelength to depth ratios.
The key difference between the 2D and 3D tests is in the reduction of first harmonic inline force contribution due to directional spreading effects. The higher harmonic forces can then be computed using the same (similar) coefficients with respect to the linear component, so the directional effects can be accounted for through a single spreading factor. We make a comparison of the peak of the linear force between unidirectional and spread seas (see Table 4), corresponding to the average of the largest events lying within the same exceedance probability range (as described earlier). In each case, we identify the 2D and 3D sea-states which are closest in description-similar H s , T p , γ , and k p d. The magnitude of the peaks from the unidirectional tests exceeds those from spread sea tests in all the five comparisons. The third column shows the ratio of the force peaks normalized by the significant wave height, and on an average, it has a value of around 0.9-this represents the effective inline kinematics factor in the spread seas and is close to what one would expect theoretically including capturing the trend that the inline kinematics factor should be smaller for greater directional spread. The normalization is done to account to difference between the 2D and 3D sea-states with respect to H s . The inline kine- matics factor cos σ θ derived using the RMS wrapped normal spreading angle supplied for the tests by DHI Denmark (from cross-correlation of wave gauge records, using data from three wave gauges to compute spreading angle at peak frequency with the extended maximum entropy method in the WAFO toolbox (Brodtkorb et al. 2000)) is shown in the righthand column of Table 4. The spreading estimates made using the wave gauge array and the inline (along the mean wave direction) component of the linear force component are consistent. With a measure of the effect of directional spreading absorbed in the definition of the linear force component F 1 , we return to the question of the effect of spreading on the higher harmonics. In Figure 3, we observe that, at upper end of the range of wave steepness, the non-dimensional amplitude and phase coefficients for the 2D and 3D spread seas become very similar in value. This implies that the effect of realistic wave directional spreading can be accounted for in the definition of the linear force F 1 and to higher force harmonics (2nd-5th) through F n 1 . These results imply that higher harmonics will be relatively less important in directionally spread seas compared to unidirectional ones. If we assume the magnitude of the coefficient of each harmonic remains the same, then the size of the nth harmonic will be reduced by f n = cos n (σ θ ). Thus, one would expect that the peak values of the largest total loads would would be reduced by a factor greater than that predicted merely by the inline kinematics factor, because the nonlinear harmonic contributions are reduced by directional spreading more than the linear component.

Sub-harmonic components
So far, we have not considered sub-harmonic components, e.g., the second-order difference force F − 2 . The magnitude of its maximum is very small compared to the total force peak, and it is reasonable to neglect this contribution from engineering analysis for bottom fixed structures. In random wave experiments, it is hard to absorb long waves and these lead to significant noise in this part of the spectrum. Our analysis of these is therefore problematic, and for brevity, we do not present details in the present paper. However, we note that our analysis suggests that the inertia loading from the 'return current' under energetic parts of the random wave train dominates the loading relative to the sub-harmonic mean load calculated for regular waves by Eatock Taylor and Hung (1985). An example of this load component extracted from potential-flow simulations is shown in Fig. 3 of Fitzgerald et al. (2014). The time-history is close to skew symmetric in time, consistent with the temporal derivative of the varying return current beneath the compact wave group used in that study, and hence with a slow Morison inertia loading. The lack of perfect skew symmetry is consistent with a much smaller mean force contribution than in 2nd-order regular waves.

Impulsive loads
When sufficiently steep wave crests are incident upon a surface piercing column with the correct phase, slamming can occur. When a slamming event occurs, the free surface col- The brackets beside each numerical value indicate the test to which it corresponds, along with the estimated spreading angle for spread sea cases. The last column shows the ratio of the force peaks (directional to unidirectional). The water depth is 33 m lides with the column causing impulsive forces in addition to the potential-flow wave-induced force components that we estimate using our method. Therefore, without explicitly accounting for this additional form of loading, some values of peak force will be underestimated. Various methods exist to identify and model slamming loads on monopile structures (see Maes et al. 2018;Ghadirian and Bredmose 2019). However, slamming loads are not something we aim to account for herein. We propose an ad hoc method to establish when a slamming event may have occurred using force measurements alone. This will allow us to identify if mismatches between our predictions and experiments are occurring as a result of slamming.
When an impulsive force is applied to the surface column, this may induce a response at the natural frequency of the structure. Therefore, we identify force peaks where a significant dynamic structural response is excited as potential slamming events. Each dynamic structural response was detected by searching for a peak at the natural frequency in the spectral representation of extreme forces. This was achieved by performing a short-duration windowed FFT centred around each extreme force. We examine experiments carried out at two different water depths, which causes the natural frequency of the model structure to change as the added mass of the systems also changes as a function of water depth. This information is used in the following comparison of measured and modelled exceedance probabilities of the inline force.

Exceedance probability
The coefficients from the harmonic decomposition can also be used to compute the maximum force for a given exceedance probability. This is achieved by first defining linear force amplitude using a Rayleigh distribution, and then using this as an input for our fitted force model to provide nonlinear corrections to fifth-order for a given sea-state and column. For each linear amplitude, the phase shift of the linear component ψ 1 is selected to maximise the total force. A comparison of the exceedance probability of the crests of force from experimental data and force computed using estimated coefficients and linear amplitude based on a Rayleigh distribution is presented in Fig. 4. Six cases are shown-the left-hand column (Exp2 and Exp4) shows unidirectional deeper water tests, the middle column (Exp9 and Exp10) shows unidirectional shallower water tests, and the right-hand column (Exp102 and Exp105) shows directionally spread cases in deeper water. The top and bottom rows show sea-states with lower and higher steepness, respectively. Force reconstructions using coefficients up to second order (F2) and fifth order (F5) are shown in the plots. There is good agreement except for differences in the largest loads (∼ top 0.5% of events for deeper water). In panel f, force peaks where a potential slamming event has been detected via structural response are indicated by a red marker. Based on this method of identification, it appears that slamming events may be the main cause of this departure.
This method of predicting forces works reasonably well even in shallow water despite the shortcomings of the Stokes-type formulation in such depth regimes. The results demonstrate the importance of the higher harmonic components in the description of extreme wave forces.

Conclusions
This paper presents an analysis of experimentally measured wave induced forces on a cylinder representative of the support structures typically used for offshore wind turbines. We analyse a range of unidirectional and directionally spread random sea-states. We analyse the data under the hypothesis that the wave-induced loads on the turbine column can be modelled using a 'Stokes-like' expansion of the linear force.
We extract amplitudes and phases of the higher force harmonics up to the fifth order. These are generally cleanly extracted with the largest uncertainty found in the third harmonic. This is consistent with previous studies which have found that localized fluid structure interaction (e.g., scattered waves or suction), can cause the extracted third harmonic to fit less well to the 'Stokes-like' model assumed here. In  (   Fig. 4 Comparison of the exceedance probability of the force crests in experiments Exp2, Exp4, Exp9, Exp10, Exp102, and Exp105. Note that Exp2 and Exp4 (left column) correspond to the unidirectional deeper water case, and Exp9 and Exp10 are in the shallower water with unidirectional waves, while the Exp102 and Exp105 (right column) are for the deeper water case with directionally spread waves. The black dots correspond to force maximums in the measured force data. The black dash-dot line shows the Rayleigh distribution. F2 is the maximum force reconstructed using linear NewWave and derived coefficients up to second order, while F5 uses coefficients up to fifth order. Red circles indicate force peaks where potential slamming events were identified by a structural response (color figure online) doing so, we demonstrate the effectiveness of the approximate method we present for the extraction of harmonic components from random signals, when applied to force measurements.
As expected, the amplitudes and phases of the different harmonics are dependent on (non-dimensional) water depth and cylinder radius. Several longer duration experiments show little dependence of the coefficients on wave steepness, consistent with the assumed model. We find that there also appears to be little dependence on the directional spread (i.e., amplitude and phases for the higher harmonics for similar unidirectional and spread sea cases are close) when the harmonic structure is written in terms of the linear component of force resolved along the mean wave direction. This suggests that for other combinations of depth, cylinder, and wave properties, unidirectional simulations or experiments would suffice, though we recommend that the duration of such tests should be > 10 4 waves. We show that, as expected from theory, the linear part of the force signal scales with the inline wave kinematic factor. However, the higher harmonics will scale by this factor to the power n and thus will be relatively less important in directionally spread seas compared to unidirectional ones, at the same return period. Our 'Stokeslike' model captures the statistics of the forces on the cylinder for all but the largest loads. The large forces not captured are likely caused by slamming events and these would need to be modelled separately.
The experiments we have used constitute a representative range of sea-state conditions taken from a pre-existing series of experiments. To build a more comprehensive force model and to fully understand the sensitivity of the derived coefficients to given sea-state parameters, a wider range of experiments is required. Multiple repeat experiments may also be needed to understand the uncertainty of the derive coefficients. For full-scale structures in-situ that have instrumentation, the method we present may also be used to extract nonlinear force components and build force models.
At present, there is no computationally fast engineering, model available to calculated the higher harmonics of wave loading. The simple approach we present here may be used to calculate very fast random time-histories, approximate loads of the correct magnitude and frequency and has the potential to be a valuable engineering tool. First, the Harmonic extraction method we present may be used to decompose the nonlinear structure of measurements of force or moment recorded on any structure, be it at full scale or in the lab, assuming the symmetry arguments consistent with a Stokes-like expansion of the response variable are appropriate. Second, decomposed measured forces for a given structure may be used to build a 'Stokes like' nonlinear force model. The normalisation we present can be used to produce force coefficients, that, for the limited range of conditions we have presented, may be used in scenarios of varying water depth and directional spreading. Accurate reconstruction of the total force time-history from the harmonic expansions validates the approach. Thus, by running unidirectional experiments, or collecting data in-situ an engineer may build a predictive nonlinear force model using the methods we have presented for a given structure. The alternative approaches to calculating nonlinear force statistics require ensembles of high-fidelity fluid structure interaction simulations, such as fully nonlinear potential flow, multi-phase CFD or SPH simulations, which are all very computationally expensive to run for many cases.

Conflict of interest
The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indi-cate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.

A Method
The method used in this work to extract nonlinear components from a random time-series is outlined in this appendix. There are three main processing steps, and they are described below along with illustrative examples. Here, we use the example of a measured force time-series; however, the method may be used to decompose any nonlinear time-series data that are expected to exhibit harmonic structure. STEP 1: Linear signal estimation The first step is to estimate the linear component of the measured force time-series. First, the force time-series is band-pass filtered to remove bound wave components. Lowand high-frequency cut-offs of 0.4 f p and 2.5 f p were used, where f p is the peak frequency. For a JONSWAP spectrum, the majority of free wave energy will lie within this range of frequencies. This filtered signal is our initial estimate of the linear force (F 1,0 ). However, the filtered signal may still contain significant second-order superharmonic bound wave components at frequencies greater than ∼1.3 f p .
To remove the second-order superharmonics that remain in the filtered signal, we make use of the following approximation. We assume that the linear signal can be written as the product of a slowly changing amplitude function (envelope) F 1 and a fast changing phase function φ as F 1 = F 1 cos φ. The second-harmonic sum contribution, comprising of an in-phase component (cos 2φ) and an out-of-phase component (sin 2φ) can be estimated from the linear signal and its Hilbert transform where α 2 and β 2 second-order coefficients, and the superscript H denotes the Hilbert transform of the linear signal F 1 . The in-phase term affects the vertical asymmetry, while the out of phase term affects the horizontal asymmetry. By removing this approximated second-order signal from the filtered time-series, we can obtain a signal which is symmetric both vertically and horizontally and hence approximately linearised. The values of coefficients α 2 and β 2 are estimated by minimising skewness of the resulting linearised signals The accuracy of this process can be assessed using the crest height statistics of the resulting time series. Figure 5 shows the probability of exceedance P exceed and ordered crests versus ordered troughs of the linearised signal (panels a and b) and its Hilbert transform (panels c and d) for Exp2. The black solid line in the figures on the left side shows the exceedance probability of the signals, with a Rayleigh distribution shown by dashed-dotted lines. The grey shaded band indicates error bounds of two the standard deviations, calculated using one thousand ensembles of crest height through bootstraping. The plots on the right show the ordered crests versus troughs (black solid line). A linear signal of sufficiently long duration will follow a line of unit slope (shown in dotted line shifted repeats. The following sub-steps are used to determine the phase shifted average extreme waveform: - A visual illustration of these steps is provided in Fig. 6. It is important that the peaks are selected from the linearized time-series instead of original time-series. The peaks in the original time-series can be influenced by non-harmonic effects and other contamination, while the phase-based separation technique is suited for applications to peaks with a fixed phase. The resultant signals are statistical averages, and unlike deterministically phased shifted experiments, they will be contaminated by some noise and energy leakages across frequency bands.

STEP 3: Extraction of harmonics
Once four phases shifted extreme waveforms are approximated, four sets of harmonic components can be extracted through the following four linear combinations of the waveforms and their Hilbert transforms: where the subscripts of the extracted force components f refer to the power in amplitude and the harmonic, respectively. Our analysis focuses only on superharmonic terms F n = f n,n . The first three combinations extract the first three harmonic components plus higher order terms, while the last combination gives the second-order difference and Fig. 8 The time-series of the different harmonics reconstructed from the spectra (shaded) in Fig. 7. The top plot shows the first harmonic, followed by the second-order sum, third-order sum, fourth-order sum, and fifth-order sum. The dashed lines are the scaled envelopes A n 1 . This figure shows data from Exp1 fourth-order superharmonic which are separable by frequency filtering.
This technique makes use of the structure of a Stokes expansion for free surface elevation. Herein, we assume that the same description may also be applied to wave loadings; this is a reasonable assumption if the transfer functions between force and surface elevation are nearly flat versus frequency, and that we can average enough extremes to remove noise. Figure 7 shows the power spectral density of the force time-series from the averaged four phase separation method described above, with the vertical axis plotted on a log scale. The harmonic structure is clearly visible in the figure, with each harmonic exhibiting distinguishable humps in the power spectrum across its frequency range of influence. We also observe a significant fifth harmonic contribution in all the experiments. The presence of higher harmonics is largely due to the nonlinear nature of the wave loadings. Note, these plots are the force power spectral density for the averaged large events, not of the underlying random process. Given the statistical nature of the method, energy leakage can be noticed around the peak frequency for the higher order contributions to the spectra, although these are negligible in magnitude compared to the linear component. Figure 8 shows the time-history of the various reconstructed harmonics from the filtered spectra in case of Exp5. The dashed-dotted lines show the wave envelopes scaled with respect to the linear wave group and they give an indication of how the higher order profiles are bound to the first order, e.g., estimate of the n-th harmonic envelope in time is derived from F n 1 = ( F 2 1 + ( F H 1 ) 2 ) n/2 , and then scaled to Reconstructed extreme force crest profiles (in red) and average of the top 100 force crests (in black) for experiments-a Exp1, b Exp2, c Exp3, d Exp4, and e Exp5. For these experiments, the water depth was 33m, and the waves were unidirectional (color figure online) the peak value of that harmonic. Note, F 1 is the linear component corresponding to the averaged profile of the extreme wave loadings, while F 1 , described previously, is the linear component of the whole time-series. Also, F H 1 is the Hilbert transform of F 1 . Figure 8, shows that the extracted sum harmonics may be considered as bound harmonics of the linear signal as they fit well within the envelopes estimated using the linear signal. We note the oscillations in the second harmonic to the right of the main wave group. These variations may be due to the linear interaction of the second-order sum free wave (an error wave generated at the wavemaker) with the structure.

B Force Reconstructions
Having extracted the average nonlinear harmonics from a random signal, and fitted coefficients to the harmonic components, it is then possible to reconstruct or predict extreme loadings using the fitted model. In this section, we provide additional examples of the force time-histories reconstructed using the fitted force coefficients and compare them to the experimental observations. Figure 9 plots the force profiles of the reconstructed force maximum (shown in red), comparisons are made to the average of the top 100 force maxima extracted from the original data. The diagram shows the cases corresponding to a water depth of 33m in full scale. There is good agreement between the data and the reconstructions across all the sea-states.

Fig. 10
Reconstructed extreme force trough profiles (in red) and average of the top 100 force minima (in black) for experiments-a Exp1, b Exp2, c Exp3, d Exp4, and e Exp5. For these experiments, the water depth was 33m, and the waves were unidirectional (color figure online) Figure 10 shows the time-history of the reconstructed force troughs (in red) and the average profile of the top 100 force troughs (in black) from the original force data, which also show good agreement.