Assessment of Practical Methods to Predict Accumulated Rotations of Monopile-Supported O ﬀ shore Wind Turbines in Cohesionless Ground Proﬁles

: Monopiles supporting o ﬀ shore wind turbines can experience permanent non-recoverable rotations (or displacements) during their lifetime due to the cyclic nature of hydrodynamic and aerodynamic loading exerted on them. Recent studies in the literature have demonstrated that conventional cyclic p–y curves recommended in di ﬀ erent codes of practice (API-RP-2GEO and DNVGL-RP-C212) may not capture the e ﬀ ects of long-term cyclic loads as they are independent of the loading proﬁle and the number of applied cycles. Several published methodologies based on laboratory scaled model tests (on sands) exist to determine the e ﬀ ect of cyclic lateral loads on the long-term behaviour of piles. The tests vary in terms of the pile behaviour (rigid or ﬂexible pile), number of applied loading cycles, and the load proﬁle (one-way or two-way loading). The best-ﬁt curves provided by these tests o ﬀ er practical and cost-e ﬃ cient methods to quantify the accumulated rotations when compared to Finite Element Method. It is therefore desirable that such methods are further developed to take into account di ﬀ erent soil types and the complex nature of the loading. under the actions of a typical 35-h (hour) storm as per the Bundesamt für Seeschi ﬀ fahrt und Hydrographie (BSH) recommendations. Using classical rain ﬂow counting, the loading time-history is discretized into load packets where each packet has a loading proﬁle and number of cycles, which then enables the computation of an equivalent number of cycles of the largest load packet. The results show that the loading proﬁle plays a detrimental role in the result of the accumulated rotation. Furthermore, ﬂexibility of the pile also has an important e ﬀ ect on the response of the pile where predictions obtained from formulations based on ﬂexible piles resulted in a much lower accumulated rotation than tests based on rigid piles. Finally, the ﬁndings of this paper are expected to contribute in the design and interpretation of future experimental frameworks for O ﬀ shore Wind Turbine (OWT) monopiles in sands, which will include a more realistic loading proﬁle, number of cycles, and relative soil to pile sti ﬀ ness.


Introduction
Offshore Wind Turbines (OWTs) are currently one of the main sources of renewable energy in the United Kingdom and other parts of Northern Europe. Recently, it has been announced that the price of electricity generated from OWTs has hit a record low and several new projects have been announced, which will enable the grid to power more than seven million homes in the near future (Harrabin [1]). Similar trends are observed in other European countries and are emerging in other countries across the world, such as the United States, China, Taiwan, and others. These reduced cost installations are the result of extensive academic and industrial efforts in research and development in the past decade.
With regard to foundations, monopiles are the dominant types of foundations for water depths up to 40 m and are even considered for deeper waters. In addition to the limit state design criteria, these foundations must satisfy certain limits of deformation during their design life. However, due to the cyclic nature of the hydrodynamic and aerodynamic loading exerted on them, the monopiles may experience permanent non-recoverable rotations, which may result in excessive tilt of the turbine, see Figure 1. If certain limits specified by building codes or the turbine manufacturer are exceeded, the warranty from the manufacturers may be lost, leading to possible financial risks and implications on wind farm developers. Therefore, ongoing research is directed towards the development of methods to understand the mechanism of this accumulation in order to develop reliable design methods soon. Due to the complex nature of the process, current design standards such as the API [2] and DNVGL [3] codes of standards do not provide a comprehensive method to calculate the accumulated rotation due to cyclic loading. Instead, current guidelines only provide p-y curves that have adjusted soil stiffness values, which are independent of the monopile size and the number of applied loading cycles. There is a consensus in the literature, such as the work of Kallehave at al. [4] and Thieken et al. [5] that the recommended curves do not represent the response of the monopile subjected to both static and/or cyclic loading as the p-y curve formulations are calibrated against slender-small diameter piles with a low number of load cycles.
Energies 2020, 13, x FOR PEER REVIEW 2 of 24 price of electricity generated from OWTs has hit a record low and several new projects have been announced, which will enable the grid to power more than seven million homes in the near future (Harrabin [1]). Similar trends are observed in other European countries and are emerging in other countries across the world, such as the United States, China, Taiwan, and others. These reduced cost installations are the result of extensive academic and industrial efforts in research and development in the past decade.
With regard to foundations, monopiles are the dominant types of foundations for water depths up to 40 m and are even considered for deeper waters. In addition to the limit state design criteria, these foundations must satisfy certain limits of deformation during their design life. However, due to the cyclic nature of the hydrodynamic and aerodynamic loading exerted on them, the monopiles may experience permanent non-recoverable rotations, which may result in excessive tilt of the turbine, see Figure 1. If certain limits specified by building codes or the turbine manufacturer are exceeded, the warranty from the manufacturers may be lost, leading to possible financial risks and implications on wind farm developers. Therefore, ongoing research is directed towards the development of methods to understand the mechanism of this accumulation in order to develop reliable design methods soon. Due to the complex nature of the process, current design standards such as the API [2] and DNVGL [3] codes of standards do not provide a comprehensive method to calculate the accumulated rotation due to cyclic loading. Instead, current guidelines only provide py curves that have adjusted soil stiffness values, which are independent of the monopile size and the number of applied loading cycles. There is a consensus in the literature, such as the work of Kallehave at al. [4] and Thieken et al. [5] that the recommended curves do not represent the response of the monopile subjected to both static and/or cyclic loading as the p-y curve formulations are calibrated against slender-small diameter piles with a low number of load cycles. To have a better understanding of the monopile response under cyclic loading, numerous scaled model tests have been carried out in the last decade for establishing methodologies for the assessment of the long-term performance (deflection and tilting) for monopiles. The tests vary in terms of the pile behaviour (rigid or flexible pile), number of applied loading cycles, and the load profile (one-way or two-way loading). The results of these tests are also presented along with best-fit curves or formulations as a means to offer practical and cost-efficient methods to quantify the accumulated rotations when compared to more complex methods such as the Finite Element Methods. It is therefore desirable that such methods are further developed to take into account different soil types and the complex nature of the loading. To have a better understanding of the monopile response under cyclic loading, numerous scaled model tests have been carried out in the last decade for establishing methodologies for the assessment of the long-term performance (deflection and tilting) for monopiles. The tests vary in terms of the pile behaviour (rigid or flexible pile), number of applied loading cycles, and the load profile (one-way or two-way loading). The results of these tests are also presented along with best-fit curves or formulations as a means to offer practical and cost-efficient methods to quantify the accumulated rotations when compared to more complex methods such as the Finite Element Methods. It is therefore desirable that such methods are further developed to take into account different soil types and the complex nature of the loading.
The main understanding that emerged from the tests, numerical work, and relevant codes of practice are as follows: Energies 2020, 13, 3915 3 of 22 (1) The parameters that control the accumulation of tilt are the loading ratios: (M min /M max ) and (M max /M R ) where M min and M max are the minimum and maximum mudline bending moments corresponding to the appropriate load combination of the wind and the wave. On the other hand, M R is the Moment of Resistance of the pile i.e., ultimate static moment capacity of the pile considering soil failure around the pile (i.e., geotechnical failure) and strictly not the structural failure of the pile through buckling or yielding. (2) Jalbi et al. [6] assessed load ratios (M min /M max and M max /M R ) for 15 operating wind turbines in European waters. The study showed that the load ratio varies widely ranging from extreme one-way loading to extreme two-way loading. There is clear mismatch of the actual load scenario and the parameters in the model tests. (3) Bundesamt für Seeschifffahrt und Hydrographie (BSH [7]) is the only code that provides guidance on the applicability of a typical OWT load time series arising from a typical storm or other loading condition such as the 35-h (hour) storm.
The main aim of this paper is to compare the performance of the available formulations under the actions of a typical 35-h storm as per the BSH [7] recommendations and the objectives are summarised as follows: (1) Present and outline the current available methods to calculate the accumulated rotation from scaled tests; (2) Present an overview of how to process a time series from the BSH [7] storm through classical rain flow counting and how these loads can be post-processed and used to calculate accumulated rotation; (3) Compare the outcome of the accumulated rotation from different sources and discuss the range of applicability of each method; (4) Compare the accumulated rotation during different periods of a storm; (5) Provide general guidance on future work in this area of research.
It is hoped that the findings from this study can be used primarily to assist in the design and interpretation of future experimental frameworks for OWT monopiles in sands that include more realistic loading profiles, numbers of cycles, and relative soil to pile stiffness.
As previously stated, the formulations presented in the published methods are typically best fit curves obtained from scaled tests of single piles under cyclic lateral loads. The tests vary in terms of the pile behavior (rigid or flexible pile), number of applied loading cycles, and the load profile (one-way or two-way loading). In reference to these properties, the following definitions pertain to this paper: (a) Rigid piles are short and bulky enough to undergo rigid body rotation in the soil under operational loads instead of bending like a clamped beam; the shear strength of the soil governs the design (i.e., the soil fails before the pile). Monopiles supporting offshore wind turbines are more likely to show rigid rather than flexible behaviour due to their large diameters. However, in very stiff soils and weak rock, flexible behaviour may be observed. (b) Flexible (or long) piles undergo deflection under operating loads and fail typically through plastic hinge formation. (c) One-way loading refers to when the ratio of the minimum to maximum load amplitude in each cycle is between 0 and 1, whilst two-way loading occurs when the ratio of the minimum to maximum load amplitude is between 0 and −1. A ratio of −1 is indicative of a full reversal in the load amplitude and sometimes termed as complete "two-way loading" in the literature.
Due to the inherent differences in the tests undertaken, careful examination of the assumptions and the applicability of each test is required prior to the selection of the appropriate formulation.

Methods to Predict Accumulated Rotation
The performance of monopiles under cyclic loads has been studied and reported in numerous research items through full-scale tests, scaled model tests, element testing, and finite element analysis. The objective of this section is to provide a brief overview of the studies on accumulated rotation and summarise the main conclusions and driving parameters of the different methods. As previously stated, the mechanism of this phenomenon is a complex process and therefore disagreements in results are evident.

Monopiles in Sands
For piles installed in sands, early tests were performed on flexible piles, which were loaded up to a small number of cycles that were not representative of recent monopiles. More recent scaled tests considered dimensions and loading conditions relative to monopile foundations and their results may be deemed representative for current applications.

The Methods of Long & Vanneste and Lin & Liao
Long and Vanneste [8] performed 34 full-scale tests applying lateral loading in sands and studied the effect of load ratio (M min /M max ), installation method, and the effect of soil density on the modulus of subgrade reaction, which is a measure of how much a pile deforms under cyclic loads. The study considered different pile sizes that exhibit mostly flexible behaviour. Based on the test results, degradation factors for calculating the soil resistance were derived, which can be used to produce adjusted non-linear p-y curves as shown in Equations (1) and (2). (2) where N is the number of loading cycles and t is the degradation parameter, which can be calculated by Equation (3).
• F L is dependent on the cyclic load ratio and was recommended to be 1.0 for one-way loading (M min /M max = 0) and 0.2 for two-way loading (M min /M max = −1), F I is dependent on the installation method and ranges from 0.9 to 1.4, F D is dependent on the soil density and ranges from 0.8 to 1.1.
Hence, it may be concluded that the main degradation parameter is the load ratio where one-way loading results in a higher degradation factor than two-way loading. It is important to note that these tests were performed for up to 50 loading cycles only, which is significantly less than the cycles experienced by a typical offshore wind foundation.
Following a similar approach, Lin and Liao [9] derived degradation parameters from 26 full-scale tests for both slender and rigid piles, most of which were slender. The degradation parameter used was dependent on the load ratio, installation method, and soil density. This best fit curve suggested by the authors is shown in Equations (4) and (5) where the degradation parameter t can be calculated using Equation (5) Energies 2020, 13, 3915 5 of 22 t =0.0032 × L P × 5 n h E p I P × F L × F I × F D (5) where: L p is the length of the monopile, E p I p is the bending stiffness of the monopile, N is the number of cycles, n h is the modulus of subgrade reaction, F L is dependent on the cyclic load ratio and is recommended to be 1.0 for one-way loading (M min /M max = 0) and 0.09 for two-way loading (M min /M max = −1), F I is dependent on the installation method and ranges from 0.3 to 1.0, F D is dependent on the soil density and ranges from 1.0 to 1.3.
The factors of Equation (4) have been updated by Li et al. [10] who suggested a t value of 0.125 (0.08 for rotation) for monopiles in sand.

The Method of LeBlanc, Houlsby & Byrne
The main drawbacks regarding the applicability of the previous methods on monopile foundations were the flexible behaviour of the tested piles and the low number of cycles. Researchers recently undertook scaled model tests to study the response of rigid piles under a higher number of loading cycles.
LeBlanc et al. [11] performed scaled model tests on rigid piles in dry sands of different relative densities. The tests were performed for up to 65,000 cycles and for both one-way and two-way loading. Based on the results, formulations were provided to predict long term rotation accumulation considering scaling between the tests and the prototypes. The accumulated rotation based on this method can be calculated from Equation (6). where: • θ s is the static rotation, which can be obtained from a non-linear analysis using standard p-y curves, T b is a function of the ratio of the maximum moment M max to the static resisting moment M R and the relative density of the sand, T c is a function of M min /M max of each loading cycle and accounts for the load profile, θ N is the rotation after N cycles of loading , α is an empirical factor that ranges between 0.3 to 0.31. Methods to calculate M R and other parameters are further elaborated in Appendix A.
Based on the original paper by LeBlanc et al. [11], T b and T c are estimated from charts as the equations are not provided. In order to include these parameters in a numerical analysis, linear piecewise functions can be derived and are further discussed in Appendix A. From the equations, it is evident that two-way loads (of about M min /M max = −0.6) resulted in high rotation accumulation of the pile head.
In this paper, the work of LeBlanc et al. [12] is extended and demonstrates how the cumulative strains due to a typical stochastic load process can be processed into Markov matrices and assessed using the "superposition" approach proposed in Stewart [13]. An equivalent number of cycles can then be obtained considering the loading history of the loads as shown in Figure 2. It is important to note that this method requires reordering the load packets in ascending order (lowest to highest accumulated rotation) and therefore it is assumed that the order of the loads has minimal effect on the final outcome. In addition, the equivalent number of cycles is calculated for each net direction independently, which is a conservative assumption. Another important aspect is that the proposed formulations are dependent on the static rotation (θ s ) and different methods will result in a different prediction of the static rotation. For instance, using Finite Element Analysis or PISA type p-y springs will result in higher foundation stiffness and produce lower values of θ s . In this paper, θ s is calculated using typical API [2] springs for sands using commercial software package OPILE [14]. For instance, assuming the first load cycle has an amplitude "a" and number of cycles Na, and the second loading cycle has an amplitude "b" and number of cycles Nb, the equivalent number of cycles of "a" and "b" can be found using Equation (7).
• θsa is the static rotation due to load packet "a", θsb is the static rotation due to load packet "b", Na number of cycles of load packet "a", Tba and Tca are the Tc and Tb factors in Equation (6) due to load packet "a", Tbb and Tcb are the Tc and Tb factors in Equation (7) due to load packet "b" The procedure may then be repeated for all other successive load levels, extracted from the rain flow count, which is elaborated upon in subsequent sections of the paper.
Thus, assuming that the largest θ is due to a cycle of amplitude "x" and number of cycles Nx, the total rotation, including the effect of accumulation of rotations during a storm, can be obtained using Equation (8).
It should be noted that Equation (8) incorporates the effects of accumulation on top of the static rotation, which is slightly different than what is proposed in Equation (6), which accumulates over the initial rotation θ1. The authors do distinguish between the initial θ1 and static rotation θs where θs was obtained from separate static tests and θ1 comes from much more rapid cyclic loading and is typically much less than θs. However, the authors use θs as it is more conservative and is easier to compute using conventional p-y analysis. It is later shown that the effect of using θs rather than θ1 has detrimental effects when considering the accumulated rotation under the action of a 35-h storm.
Another important aspect discussed in the aforementioned paper is the effect of load reversals. The experimental results showed that 2.51 N to 1.96 N cycles are required to neutralize the effect of loading reversal. In other words, if a cyclic load with certain values Mmin/Mmax, Mmax/MR profile is applied to the monopile head in the positive X direction for N cycles, this would result in the rotation  [13] applied on the accumulation model of LeBlanc et al. [12].
For instance, assuming the first load cycle has an amplitude "a" and number of cycles N a , and the second loading cycle has an amplitude "b" and number of cycles N b , the equivalent number of cycles of "a" and "b" can be found using Equation (7).
• θ sa is the static rotation due to load packet "a", θ sb is the static rotation due to load packet "b", N a number of cycles of load packet "a", T ba and T ca are the T c and T b factors in Equation (6) due to load packet "a", T bb and T cb are the T c and T b factors in Equation (7) due to load packet "b" The procedure may then be repeated for all other successive load levels, extracted from the rain flow count, which is elaborated upon in subsequent sections of the paper.
Thus, assuming that the largest θ S is due to a cycle of amplitude "x" and number of cycles N x , the total rotation, including the effect of accumulation of rotations during a storm, can be obtained using Equation (8).
It should be noted that Equation (8) incorporates the effects of accumulation on top of the static rotation, which is slightly different than what is proposed in Equation (6), which accumulates over the initial rotation θ 1 . The authors do distinguish between the initial θ 1 and static rotation θ s where θ s was obtained from separate static tests and θ 1 comes from much more rapid cyclic loading and is typically much less than θ s . However, the authors use θ s as it is more conservative and is easier to compute using conventional p-y analysis. It is later shown that the effect of using θ s rather than θ 1 has detrimental effects when considering the accumulated rotation under the action of a 35-h storm.
Another important aspect discussed in the aforementioned paper is the effect of load reversals. The experimental results showed that 2.51 N to 1.96 N cycles are required to neutralize the effect of loading reversal. In other words, if a cyclic load with certain values M min /M max , M max /M R profile is applied to the monopile head in the positive X direction for N cycles, this would result in the rotation of θ 1 . Consequently, if the same load is applied for N cycles in the negative X direction, this would result in a rotation of θ 2 , which is in the opposite direction of θ 1 . The model tests applied by the This effect of load reversals is shown in Figure 3, which is extracted from the references and shows how 2.51 × N cycles are required to restore the rotation. A conservative way to derive the applied net cycles is therefore to assume that one cycle of a specific load packet with the same magnitude and opposite direction is enough to reverse the original same magnitude load packet. of θ1. Consequently, if the same load is applied for N cycles in the negative X direction, this would result in a rotation of θ2, which is in the opposite direction of θ1. The model tests applied by the authors have shown that it would require 2.51N to 1.96N cycles in the positive direction to bring the monopile back to its original of position of θ1 prior to load reversal.
This effect of load reversals is shown in Figure 3, which is extracted from the references and shows how 2.51 x N cycles are required to restore the rotation. A conservative way to derive the applied net cycles is therefore to assume that one cycle of a specific load packet with the same magnitude and opposite direction is enough to reverse the original same magnitude load packet.

Additional references in sand
This section briefly presents other recent references for the calculation of the accumulated rotation in sands, which are expected to further modify and influence the widely used methods discussed above.
Rudolph et al. [15] studied the effect of multi-directionality under one-way loads in centrifuge tests in sand. The study concluded that multi-directional one-way loads increased the accumulation of displacements.
Nanda et al. [16] performed scaled model tests for monopiles in sand for loading of up to 1000 cycles. The experimental setup contained a gearbox that changed the directionality (changing the load application angle) of both one-way and two-way cyclic loads. The results showed that both

Additional References in Sand
This section briefly presents other recent references for the calculation of the accumulated rotation in sands, which are expected to further modify and influence the widely used methods discussed above.
Rudolph et al. [15] studied the effect of multi-directionality under one-way loads in centrifuge tests in sand. The study concluded that multi-directional one-way loads increased the accumulation of displacements.
Nanda et al. [16] performed scaled model tests for monopiles in sand for loading of up to 1000 cycles. The experimental setup contained a gearbox that changed the directionality (changing the load application angle) of both one-way and two-way cyclic loads. The results showed that both multi-directional and uni-directional one-way cyclic loads produced a higher horizontal displacement at the monopile head. It was also found that multi-directional loads in both cases resulted in higher accumulated horizontal deformations.
Considering the above, it is evident that there is no clear conclusion in the literature on whether one-way or two-way loading is more critical for the accumulated rotation of a monopile. The differences in the results might be due to a number of reasons such as the size and scale of each test, the number of loading cycles, the rate (frequency) of loading, the amplitude (intensity) of loading, the installation method, and the total execution time of each test. Furthermore, the model testing setup is usually solely Energies 2020, 13, 3915 8 of 22 either one-way or two-way, whereas the loading on an offshore wind foundation is a combination of the two with cyclic load packets of different means and ranges. However, recent scaled tests by Richards et al. [17], Truong et al. [18], and Abiker and Achmus [19] show that the effect of two-way loading is less severe than the ones previously reported in LeBlanc et al. [11]. It is shown that the load factor for partial two-way loading factor T c can be reduced from 4, as originally presented, to a number in the range of 2-3, which has a significant impact on the effect of accumulated rotation considering the total number of cycles in a given storm. However, in this study, the T c factors from the latter are used to compare the effect of including different parts of the BSH storm.
There is a consensus in the literature that rigid piles cause higher degradation (more accumulated rotation) than flexible piles because soil deformation and degradation of soil stiffness occur over the entire length of the pile (Abiker and Achmus [19]). Hence, tests based on flexible piles such as Lin and Liao [9] and Long and Vanneste [8] result in a low value of accumulation considering that the piles tested were flexible and under a small number of cycles. Li et al. [10] suggested Equation (9) as an adjustment for the Lin and Liao [9] formulation as: Hence Nb* between packets "a" and "b" can be found using Equation (10) Nb * = exp 1 0.08

Modification of the LeBlanc, Houlsby & Byrne Method
The method by LeBlanc et al. [11] accumulates the rotation above the static rotation θ s , which is stated to be a conservative approach. It has been of interest in this study to further investigate the sensitivity of the results when the accumulation is undertaken over θ 1 rather than θ s . This was performed in order to clarify whether a better understanding of the rotation of the first cycle is required or whether the static rotation is sufficient for the analysis. Therefore, slight modifications to the formulation are suggested, which are summarised below: • Use of the superposition method on the total rotation θ rather than the accumulated rotation ∆θ. This is following the work of Lin and Liao [9] and Li et al. [10], which used this approach to find the total strain on a laterally loaded pile. • Accumulation over the initial rotation θ 1 rather than the static rotation θ s .
The same methodology presented in Section 2 was adopted to calculate the equivalent number of cycles. Specifically, the equivalent number of cycles between the two load packets can be calculated as: where θ 1a and θ 1b are the total permanent rotations due to the first load cycles "a" and "b" and the other parameters have the same meaning as in Equation (7). The procedure may then be repeated for all other successive load levels extracted from the rain flow count.
Currently, there is no standard method to calculate the rotation due to the first cycle and the accumulation equation is implemented assuming n = 1 as shown in Equation (12). For the purpose of an initial understanding of the effects, this assumption is reasonable since all the results are based on these best fit equations.
Thus, if the largest θ 1 is due to a cycle of amplitude "x" and number of cycles N x , two calculations can be made: Maximum rotation experienced by the monopile (note the rotation is added to the static): Rotation after restoration (note that the accumulation is added to θ 1 ):

Monopiles in Clay
Though the main research focus of this paper is in sands, this section provides a brief insight of the research in cohesive ground profiles. It is noted that tests of piles under cyclic loading in clay are less common, possibly due to the difficulty of modelling undrained conditions in laboratory scaled model tests. It is common practice to ensure that a certain level of threshold shear strain is not exceeded, which can lead to conservativism in the assessment. In other words, clays are known to have unrecoverable degradation after a certain strain called the threshold shear strain and conservative designs in clayey ground profiles ensure that this threshold strain is not reached. Other methods include the work of Carswell et al. [20], which presented a numerical model to adjust the p-y curves of clay to account for cyclic loading. The ultimate resistance p u is reduced based on the number of cycles and the displacement caused by each load packet through the application of a degradation factor λ N . This factor was used to reduce p u as shown in Equation (15) and (16). where: • p uN is the ultimate resistance after N load cycles, D is the diameter of the monopile, y static is the static deformation of the monopile It should be noted that this model does not explicitly consider the load profile, i.e., whether the load is one-way or two-way. The reference mentions that the p-y models on which this method was based were derived for one-way testing and the applicability of this method in two-way loading was not discussed. A similar approach has been adopted by Zhang et al. [21] for the clay layers to calculate the damage to p u by the incorporation of the results of Direct Simple Shear (DSS) and Cyclic Triaxial Tests and deriving cyclic p-y springs that are dependent on the number of cycles. Lombardi et al. [22] performed a scaled test of a steel pipe in clay supporting a wind turbine generator. Subsequently, cyclic loads at certain loading amplitudes were applied using an actuator mounted on the tower. It was shown that the natural frequency decreased with increasing number of cycles especially when the loading amplitudes were high (imposing higher shear strains on surrounding soils), which is indicative of the cyclic degradation of the surrounding clay. This section emphasises the need of further research and model testing in clays, especially with the emergence of offshore wind sites in east of China and Taiwan, where soft clayey deposits are commonly found.

Summary of Literature
Judging from the above, there has been an increasing number of tests in the literature regarding the accumulated rotation. These tests have enhanced the knowledge of this engineering problem and provide useful best-fit formulations for preliminary design. However, it is evident that there is a disagreement in the literature on whether one-way or two-way loading is more detrimental. Even in more recent literature, which predominantly focused on rigid piles, agreeing that two-way loads are more detrimental than one-way loads, disagreements still arise on the intensity of the damage of the two-way loads. This is further clarified by the summary provided in Table 1.

The differences in results
Energies 2020, 13, 3915 10 of 22 might be due to multiple reasons, including the size and scale of each test, number of loading cycles, the rate of loading, the amplitude (intensity) of loading, the installation method, and the total time of each test. One-way loading results in a higher degradation parameter to the p-y formulation Tests were performed on flexible piles for less than 100 loading cycles, which is less than what is experienced by an offshore WTG Lin & Liao [9] Full scale tests in sands One-way loading results in a higher degradation parameter, which leads to a higher pile head deformation Tests were mostly on flexible piles and loading for less than 100 loading cycles, which is less than what is experienced by an offshore WTG LeBlanc et al. [11] Laboratory scaled tests for rigid monopiles in sand Partial two-way loading results in a much higher degradation and leads to more rotation accumulation than one-way loads Tests were performed on a rigid pile in sands for up to 65,000 loading cycles, which is more representative of an offshore WTG environment Li et al. [10] Field tests on small scale to intermediate rigid piles in sand Only performed tests for one-way load and did not assess the damage for two-way loads Tests were performed on piles in sands for 5000 cycles Albiker & Achmus [19] Presented a model for stiffness degradation in a finite element model Applicable to one-way loading in drained conditions Method compared to and validated against Li et al. [10] Richards et al. [17] Laboratory scaled tests for rigid monopiles in sand Partial two-way loading results in a much higher degradation and leads to more rotation accumulation than one-way loads, however with a smaller factor than LeBlanc et al. [11] Tests were performed on a rigid pile in sands for up to 10,000 loading cycles. The response of the monopile to multi-directional storm loading was also studied

Computation of Mudline Forces and Moments
This section briefly describes the computation of the load required for the analysis. The primary loads at the mudline due to the wave and wind loading are the base shear and the corresponding overturning moment. For use in the accumulated rotation assessment, the time-series loading history has to be converted to equivalent discrete load packets of shear force and overturning moment. For each combination of force and moment load packet, a corresponding number of cycles is required.
In order to determine the number of cycles, a Rainflow algorithm was used. Rainflow counting is a popular method of analysis typically used in fatigue analysis in order to decompose a time-series into simplified load packets with a mean and range, see Figure 4 for a schematic. It was primarily developed by Matsuishi and Endo [23], with numerous modifications such as the one proposed by Rychlik [24]. The outcome of the rainflow count is Mmin, Mmax, and number of cycles. The method can be extended to compute the accumulated rotation as previously discussed in LeBlanc et al. [11] with the equivalent number of cycles being calculated as per Figure 4, where the resulting matrix is generally called a Markov matrix. There are also numerous algorithms available in the public domain that are programmable in open source software, which makes the implementation of these methods relatively more cost effective when compared to a time-series simulation on a finite element solver.
has to be converted to equivalent discrete load packets of shear force and overturning moment. For each combination of force and moment load packet, a corresponding number of cycles is required.
In order to determine the number of cycles, a Rainflow algorithm was used. Rainflow counting is a popular method of analysis typically used in fatigue analysis in order to decompose a time-series into simplified load packets with a mean and range, see Figure 4 for a schematic. It was primarily developed by Matsuishi and Endo [23], with numerous modifications such as the one proposed by Rychlik [24]. The outcome of the rainflow count is Mmin, Mmax, and number of cycles. The method can be extended to compute the accumulated rotation as previously discussed in LeBlanc [11] et al. with the equivalent number of cycles being calculated as per Figure 4, where the resulting matrix is generally called a Markov matrix. There are also numerous algorithms available in the public domain that are programmable in open source software, which makes the implementation of these methods relatively more cost effective when compared to a time-series simulation on a finite element solver. The algorithm decomposes a time-series to a number of discrete load packets with different mean values, ranges and numbers of associated cycles. When considering mudline shear force and overturning moment, there is no direct correlation between a base shear load range of a certain number of cycles, and an overturning moment load range of different number of cycles. Therefore, an assumption needs to be made in order to establish a correlation between the base shear and the overturning moment using their extreme values across the simulations. The algorithm decomposes a time-series to a number of discrete load packets with different mean values, ranges and numbers of associated cycles. When considering mudline shear force and overturning moment, there is no direct correlation between a base shear load range of a certain number of cycles, and an overturning moment load range of different number of cycles. Therefore, an assumption needs to be made in order to establish a correlation between the base shear and the overturning moment using their extreme values across the simulations.
Two different approaches are commonly used to establish this force-moment correlation: (1) Determine an "eccentricity" value to convert the mudline force to moment. This is calculated as the maximum value of the ratios of the maximum moment over the maximum shear force within each design event duration. (2) A regression analysis of the maxima of forces and moments for the different loading conditions and the establishment of a relationship between the two based on a "best-fit" approach. This is similar to the approach adopted in Carswell et al. [20] and Wang and Larsen [25].
Both approaches have been used in this assessment and their influence on the results is discussed in subsequent sections of this paper.

The Time Series
For this assessment, the mudline loads were divided into four design event durations as per the BSH, which typically include the following and are shown in Figure 5: have higher lever arm than an extreme storm as the predominant load (wind) is applied at a greater distance from the mudline. This depends on the governing loading source and the point of rotation of the monopile. In Figure 5, the wind bins represent the wind velocity vw (in m/s) in the time frame of the bin itself (they are also referred to as Wind Steps). Similarly, the Hs bin represents the significant wave height in the same time frame.

Analysis Example
An analysis example is taken to allow a comparison between the methods discussed in Section 2. A monopile with embedment length Lp = 30.5m and diameter Dp = 8m and a uniform wall thickness tf = 70mm is taken. The monopile is assumed to be installed in sandy soil with φ' = 31.5° and the relative density is taken as 75%, which is the equivalent of 38% of the tests discussed in LeBlanc et al. [11]. This separation was undertaken because, depending on the loading case, the eccentricity and the correlation between the force and the moment are different. For example, the fatigue loading may have higher lever arm than an extreme storm as the predominant load (wind) is applied at a greater distance from the mudline. This depends on the governing loading source and the point of rotation of the monopile. In Figure 5, the wind bins represent the wind velocity vw (in m/s) in the time frame of the bin itself (they are also referred to as Wind Steps). Similarly, the Hs bin represents the significant wave height in the same time frame.

Analysis Example
An analysis example is taken to allow a comparison between the methods discussed in Section 2. A monopile with embedment length Lp = 30.5 m and diameter Dp = 8 m and a uniform wall thickness tf = 70 mm is taken. The monopile is assumed to be installed in sandy soil with ϕ' = 31.5 • and the relative density is taken as 75%, which is the equivalent of 38% of the tests discussed in LeBlanc et al. [11].
These are preliminary sizes that are usually found in the German North Sea installations. The wind turbine used for the analysis is an 8 MW wind turbine. Typically, wind interface loads are provided by turbine manufacturers and are project-specific/confidential. A generic wind load has been utilised for this assessment, which represents the location of the foundation discussed above. Markov Matrices were generated, and both the maximum eccentricity and regression models were used as discussed above. It is reminded that the main objective of this manuscript is to establish the areas of research required for simplified methods of quantifying the accumulated rotation for design. Figure 6 shows a description of the methodology required to calculate the accumulated rotation for the monopile foundation presented above. The following methods have been considered: • A simplified method of assessing long-term tilt using a loading/unloading model based on cyclic p-y curves obtained from the API [2] and DNVGL [3] codes of practice. This is termed Method 1. The unloading stiffness is assumed to have the same slope as the loading stiffness, which is conservative considering the studies of Achmus et al. [26] (see, for example, Figure 7 for a schematic of the method).

•
The LeBlanc et al. [11] method, which considers random two-way lateral loads in sands. This is termed Method 2.

•
The modified LeBlanc et al. [11] method (the use of θ 1 ) to investigate the sensitivity of the results to the parameters used. This is termed Method 3.

•
The LeBlanc et al. [11] method, which considers random two-way lateral loads in sands. This is termed Method 2.

•
The modified LeBlanc et al. [11] method (the use of θ1) to investigate the sensitivity of the results to the parameters used. This is termed Method 3.
The methods to calculate the factors of Method 2 and Method 3 are found in Appendix A. Figure 6. Description of the methodology used to calculate the accumulated rotation. Figure 6. Description of the methodology used to calculate the accumulated rotation.

Comparison Between Methods 1,2, and 3
The accumulated rotation was calculated based on the presented methods. Loading packets were derived from the Markov matrices and applied to the updated non-linear soil model described above. Subsequently, the static rotations were obtained using OPILE and implemented to the formulations described in Section 2 to obtain the equivalent number of cycles and accumulated rotation. The loading packets were obtained using both the maximum eccentricity method and the regression method.

Comparison Between Methods 1, 2, and 3
The accumulated rotation was calculated based on the presented methods. Loading packets were derived from the Markov matrices and applied to the updated non-linear soil model described above. Subsequently, the static rotations were obtained using OPILE and implemented to the formulations described in Section 2 to obtain the equivalent number of cycles and accumulated rotation. The loading packets were obtained using both the maximum eccentricity method and the regression method.
As previously stated, it was assumed that loads leading to deformations in one direction are independent of those in the opposite direction i.e., positive and negative deformations/rotations were assessed independently, which is conservative. Table 2 summarises the results for a 35-h BSH storm considering the maximum eccentricity. Figures 8 and 9 depict the implementation of the Superposition Model ( Figure 2) to obtain the equivalent number of cycles resulting from using Method 2 and Method 3, respectively.    From Table 2, it is observed that the cyclic p-y springs provided in Method 1 underpredict the accumulated rotation when compared to Method 2 and Method 3. Methods 2 and 3 predict a significantly higher maximum rotation after cycling as each load packet uses the maximum  Table 2, it is observed that the cyclic p-y springs provided in Method 1 underpredict the accumulated rotation when compared to Method 2 and Method 3. Methods 2 and 3 predict a significantly higher maximum rotation after cycling as each load packet uses the maximum eccentricity between the base shear, H and mudline overturning moment, M. The results from this analysis can be considered a highly conservative "upper-bound". In addition, Method 3 predicts a lower value than Method 2 as it has been modified to accumulate over the initial rotation θ 1 rather than the static rotation θ s , which slightly reduces the conservatism of Method 2. Table 3 below summarises the results for the same process but using a regression "trend" between base shear and overturning moment. Comparing the two tables, the maximum accumulated rotation is reduced along with the accumulated and restored rotations. It is also worth noting that an analysis was undertaken using the "net" number of cycles (reducing the number of cycles for packets with exact but opposite loading patterns, as discussed in Section 2) and it was found that this had a negligible impact on the results. This is because only the packets with identical magnitudes, but opposite directions, are taken into account. Further research is needed for the net amount.

Comparison to Other References
The same superposition model (Figure 2) was applied to Li et al. [10] formulations presented in Equations (9) and (10) and resulted in a maximum rotation θ N = 0.56 • . This is lower than both Method 2 and Method 3 because Equation (10) was calibrated for piles under one-way loading, which has lower damage implications than the two-way loads in Method 2 and Method 3 (through the value of T c in Equation (6)). It may be noted that the highest load packet from the extreme wave (which results in the static rotation of 0.52 • ) is two-way as shown in Table 2. This means Method 2 and Method 3 will always result in higher rotations than Equation (10) and this is due to the disagreement in the literature when considering the severity of one-way and two-way loads. In addition, it is realised from a quick implementation of Equation (4) and (5) that the damage parameter t is even lower than the value of 0.08 suggested by Li et al. [10], which is an indicator that the formulations originally based on flexible piles may not be used on rigid piles as they result in much lower accumulation of the rotation.
Comparing the results from all the different methods, the following can be concluded: • Method 1 (direct implementation of cyclic p-y curves) does not capture the effect of the loading profile and is shown to result in a lower value of both the maximum and the accumulated rotation. • Method 2 will result in higher rotations than Method 3 and this is because the accumulation in Method 2 occurs over the static values θ s , whilst Method 3 considers the initial value θ 1 . Therefore, due to the noticeable difference between the two Methods, it is encouraged that further research is needed in addressing a practical method to calculate θ 1 . • Method 2 and Method 3 are likely to show higher rotations than other sources in the literature such as Li et al. [10]. This is because these methods consider two-way loading to be more onerous for accumulated rotation than one-way loading. • Methods 2, 3 and Li et al. [10] are expected to predict higher deformations than Lin and Liao [9] and Long and Vanneste [8] as the piles considered in these tests were rigid. • All methods considered the most "dominant loading profile", where in reality the loads in opposing directions will counteract the movement and will possibly result in a lower rotation

Comparison Between Design Event Durations
This section aims at demonstrating the difference between the results when considering the four design event durations as discussed in Section 3, namely: • A 10-min simulation with an extreme wave; A 3-h peak storm; A full 35-h BSH storm; A full 35-h BSH storm and 25 Years of fatigue loading.
The aim of the assessment is to demonstrate the influence of the design event duration on the accumulated rotation and help illustrate the level of computational complexity that is required for future assessment.
Method 3 has been selected as the basis of this comparison and both the "maximum eccentricity" and "regression analysis" have been performed. The results for using the maximum eccentricity are presented in Tables 4 and 5 for the regression model. For both eccentricity methods, it can be observed that the levels of accumulated rotation increase, with increasing design event durations.
When considering the maximum eccentricity in the analysis, it is observed that the increase in accumulated rotation for each increase in design event duration is significant. This is principally due to the over-conservatism in this method, particularly when considering the fatigue data. The levels of accumulated rotation are not considered to be representative and should only be considered as an upper bound.
When considering the regression model, the increase in accumulated rotation with the increase in design event duration is more reasonable. As expected, the majority (approximately 85% of the total) of the accumulated rotation comes from the 10-min simulation with an extreme wave. A further 11% of the total is added when considering the 3-h peak storm (i.e., 96% of the total). Considering the full BSH [7] storm and the 25 years of fatigue loading makes up the remaining 3% and 1%, respectively. This is similar to the results obtained by Wang and Larsen [25] when analysing monopiles under storm loading.

General Conclusions and Recommendations
The results of this study are predominantly based on tests available in academic references for monopiles in homogenous sand. The following can be concluded and recommended based on this study:

•
The loading profile (M min /M max ) plays a detrimental role in the result of the accumulated rotation. Different tests from different sources of literature debate the effect of the loading profile on the long term-tilt. Without other supporting data, it is recommended to follow the methods that would provide conservative estimates in terms of loading direction until further developments are available in the near future.

•
The flexibility of the pile also has an important effect on the response of the pile. Results obtained from formulations based on flexible piles resulted in a much lower accumulated rotation than tests based on rigid piles.

•
When selecting a design event duration, it can be observed that the 10-min time series with extreme wave predicts over 80% of the accumulated rotation. This is considered acceptable for early stages of the design. It may be necessary to consider at least the 3-h peak storm in detailed design to account for the majority, i.e. >96%, of the accumulated rotation. It is recommended that, if the turbine is sensitive to SLS displacements, then a small number of studies are carried out using the full 35-h BSH storm and/or the fatigue loading.

•
If fatigue loading is to be considered in the assessment, the regression model is recommended as the maximum eccentricity method is believed to over-predict the accumulated rotations.

•
Accumulating over θ 1 resulted in a significant reduction when compared with accumulating over θ s . However, this would need further validation through testing and advanced analysis.
As previously stated, the methods discussed in this paper present an efficient and cost-effective way to estimate the accumulated rotation. Therefore, it is hoped that the discussed conclusions pave the way of future research such that these methods are further enhanced for more accurate analyses. The following are some of the areas that require further research:

•
In general, there is more research directed to sands than clays. With the emergence of offshore markets in other areas around the world, it is recommended that some additional insight must be available on the performance of monopiles in clay. There are also no insights on the performance in two-layered or multi-layered ground profiles.

•
There has not been extended work on understanding how the order of loads affects the accumulated rotation of the monopile. The N eq method rearranges the load packets from largest to smallest, which is acceptable to storm type loading. However, it is recommended that a mixture of loads are applied in scaled tests to gain a further understanding of the response.

•
The load reversal of identical but opposite (in direction) loads has been touched upon in the literature. However, as previously discussed, it is current practice to analyse each loading profile separately and it is not clear how to incorporate the "net" effect of the positive and negative loads.

•
It is evident that using the static rotation for accumulated rotation calculations results in higher values than other methods. It is therefore recommended that guidance should be available on different methods to predict the value of θ 1 and how to include it in design calculations.

•
Loads applied in scaled tests available in the literature range from extremely two-way to extremely one-way. Analysis of the loads shows that the highest load packets of a storm are partially two-way. In addition, other parameters such as M max /M R should be more reflective of the actual loads applied in an offshore wind farm.
It is important to note that the density values in the formulations above are the relative densities in the scaled model tests rather than the in-situ relative density. The effective stress the authors propose is that at 0.8 x monopile embedment depth. Other similar formulations on this issue can be found in literature such as the work of Bolton [27].
It may be noted that the formulations provided by the authors of the publication only include laboratory relative densities of 4% and 38%, which correspond to 8% and 75% in terms of realistic in-situ conditions. However, it is evident that the formulations to compute Tb for both 4% and 38% are linear and, seemingly, linear interpolation may be used for values in between and linear extrapolation for values greater than 38%. This is a preliminary estimate based on the available data, which needs to be confirmed with further testing. ζ b is an important parameter to enable the calculation of T b . Consider Figure A1, where a certain load packet H max , M max is plotted. In addition, a failure line can be established between the maximum base shear and maximum overturning moment the monopile can take. As per the figure, ζ b can be calculated as In the figure, H R is the maximum lateral load the monopile can take without any moment. Similarly, M R is the maximum over-turning moment the monopile can take without a lateral load. M R , H R , and the failure line can be obtained through the static equilibrium of the ultimate force, as per Poulos and Davies [28], and is described below. Consider the monopile shown in Figure A2, which is embedded in a sandy ground profile. The Figure also shows the ultimate failure pressure variation with depth and assumes a linear idealization. However, it is important to note that the ultimate failure pressure shown in Figure A2 is a simplified idealisation of the actual distribution. The linear idealisation was implemented in this Appendix in order to simplify subsequent derivations of the failure envelope.
The failure envelopes from which MR is calculated from throughout this manuscript are based on non-linear failure pressure distributions as per the recommended formulations of API [2] and DNVGL [3]. Consider the monopile shown in Figure A2, which is embedded in a sandy ground profile. The Figure also shows the ultimate failure pressure variation with depth and assumes a linear idealization. However, it is important to note that the ultimate failure pressure shown in Figure A2 is a simplified idealisation of the actual distribution. The linear idealisation was implemented in this Appendix in order to simplify subsequent derivations of the failure envelope. The failure envelopes from which MR is calculated from throughout this manuscript are based on non-linear failure pressure distributions as per the recommended formulations of API [2] and DNVGL [3].
Consider the monopile shown in Figure A2, which is embedded in a sandy ground profile. The Figure also shows the ultimate failure pressure variation with depth and assumes a linear idealization. However, it is important to note that the ultimate failure pressure shown in Figure A2 is a simplified idealisation of the actual distribution. The linear idealisation was implemented in this Appendix in order to simplify subsequent derivations of the failure envelope.
The failure envelopes from which MR is calculated from throughout this manuscript are based on non-linear failure pressure distributions as per the recommended formulations of API [2] and DNVGL [3]. The following notations apply to the figure: • L p is the embedment depth of the monopile. • d is the point of rotation at a certain load level. • P 1 and P 2 are the soil failure pressures at the point of rotation "d". They are equal in magnitude and opposite in direction. • P 3 is the soil failure pressure at the pile toe. A B-y spring can also be added in this location if needed. • H is the external applied shear at the mudline and M is the external applied overturning moment.
The horizontal equilibrium of the system can be taken as: H − 1 2 ×P 1 ×d + 1 2 × (P 3 − P 2 ) × (L P − d) + (P 2 × (L P − d))= 0 Similarly, the moment equilibrium of the system can be calculated by taking the summation of the moments at the point of rotation d: It may be noted that to calculate M and H for the equilibrium equation above, the only unknown is the point of rotation "d". To plot the M-H envelope, d should be varied from 0 to Lp and the resulting pairs of H and M (H, M) should be recorded. Figure A3 is the result of plotting all the combinations of (H, M) at different d values.
A simpler method would be to load the p-y model to the points (0,H R ) and (M R ,0), respectively, and record the failure point. It may be noted that to calculate M and H for the equilibrium equation above, the only unknown is the point of rotation "d". To plot the M-H envelope, d should be varied from 0 to Lp and the resulting pairs of H and M (H, M) should be recorded. Figure A-3 is the result of plotting all the combinations of (H, M) at different d values. Figure A3. Failure Envelope.