Evaluation of soil models for improved design of offshore wind turbine foundations in dense sand

In order to perform optimised and safe design of foundations for offshore wind turbines (OWT), it is important to have calculation tools that describe the key features of water-saturated soil subjected to complex and irregular loading over a wide range of strain levels. Soils subjected to cyclic loading are prone to strain accumulation. The accumulated (plastic) volumetric strain may result in excess pore pressure or stress relaxation, which will reduce the effective stresses, stiffness and strength of the material. Strain accumulation in dense sand is a complex mechanism of deformation and it is challenging to describe it properly. Four different soil models to describe the stress–strain relationships of dense sand are evaluated in this paper: two implicit models that follow the actual stress history and two explicit models that calculate the accumulated strains as a function of number of cycles. These models are first evaluated on the basis of their theoretical framework and back-calculations of laboratory tests specifically carried out for the design of OWT foundations in dense sand. Second, the models are implemented in finite-element analyses and evaluated on the basis of the analyses of an OWT monopile subjected to different loading conditions.


INTRODUCTION
Offshore wind turbines (OWTs) are dynamically sensitive structures subjected to complex cyclic loading conditions from wind, current and sea waves, as well as loads due to the operation of the turbine. The foundation is an essential part of the system as it has to resist the loads from the structure above and remain stable and functional during the entire lifetime of the OWT. In addition, the stiffness and damping characteristics of the foundation have a substantial impact on the global dynamics (Page et al., 2019).
The design of foundations for OWTs is driven by the requirements for capacity, installation and operational performance that follow the limit state design philosophy . These include: ultimate limit state (ULS) for the assessment of foundation capacity; serviceability limit state (SLS) for the prediction of the accumulated foundation displacement and tilt due to the cyclic loading experienced during the OWT lifetime; and fatigue limit state (FLS) for the assessment of structural fatigue. Predicting the foundation response in these three limit states requires that the soil model be used accurately to reproduce the cyclic soil behaviour. In particular, the conditions where the accumulation of strain and pore pressure can change the stiffness and strength of soil are crucial, because these changes will affect the foundation capacity in the ULS, the prediction of accumulated displacements and tilt in the SLS and the foundation stiffness in the FLS assessment.
The behaviour of piles under cyclic loading has been studied differently depending on the soil type. This is mostly due to assumptions considered in offshore oil and gas standards (e.g. API, 2010), where modelling the behaviour of clay as undrained and of sands as drained under cyclic loading conditions might be representative for the geometry of commonly used offshore piles.
The behaviour of saturated sand under cyclic loading conditions is typically neither fully drained nor perfectly undrained. The sand is therefore generally under partially drained conditions. This is especially true in the case of large-diameter monopiles in dense sand, as illustrated in Zhang et al. (2019). Based on a numerical study on the drainage conditions of a monopile in dense sand subjected to different cyclic loads, the dense sand is found to be close to an undrained condition during a single load cycle with a period of some few seconds. The results challenge the assumption that the behaviour of dense sand under cyclic loading is drained.
Several researchers have investigated the behaviour of monopiles in sand under cyclic loading conditions, either experimentally or numerically. Most experimental studies have focused on dry sand (Leblanc et al., 2010;Cuéllar et al., 2012;Klinkvort & Hededal, 2013;Bayton et al, 2018) or fully drained conditions (Burd et al., 2017). Numerical studies of monopiles in sand have been performed assuming both drained and undrained conditions (Achmus et al., 2009;Zdravković et al., 2015;Corciulo et al., 2017;Sheil & McCabe, 2017;Kementzetzidis et al., 2019). However, the results depend strongly on the constitutive model used to represent the sand behaviour. For example, Sheil & McCabe (2017) compared simulations of a monopile in sand utilising two different constitutive models, and found substantial differences in the computed results. These variations in the predicted foundation response highlight: (a) the need to understand better the differences among constitutive sand models implemented in numerical analyses to predict the cyclic behaviour of OWTs; (b) the need to assess the accuracy of the models in their description of dense sand behaviour observed in laboratory tests for drained, undrained and partially drained conditions. This paper discusses the laboratory behaviour of a dense sand from a real offshore wind site and studies the required features of the material models to be used in numerical analysis of foundations for OWTs. A reliable description of the behaviour of dense sand for the foundations of OWTs is a relevant topic because large offshore wind parks are planned at sites dominated by dense sandfor example, Dogger Bank (Fitch et al., 2005).
To evaluate the effect of different features in the modelling of dense sand, four material models were considered: two implicit models that can follow the actual irregular load history, and two explicit models that can calculate accumulated strain and pore pressure as a result of the number of load cycles around an average load level. The performance of each of the four models was compared with laboratory test results on the dense sand from the OWT location in the North Sea. The strengths and limitations of the models are highlighted, and the prediction performance of each of the models for monopile response in dense sand is illustrated.

DESIGN LOADS
The stiffness of the foundation may affect the structural loads and vice versa. It is therefore essential to identify the governing load cases for both structural and geotechnical design. For that purpose, standards like the International Electrotechnical Commission's IEC 61400-3 (IEC, 2019) define so-called design load cases (DLCs). The load cases are, however, focused on the structural design rather than the geotechnical design. DNV GL (2016) went a step further and categorised the DLCs in the ULS, FLS and SLS. The load effects from the wind and sea waves on the rotor nacelle assembly, tower, support structure and monopile are calculated using integrated coupled aerodynamic, hydrodynamic and structural codes. Input to these analyses are the DLCs and corresponding metocean data. These analyses also need to account for the stiffness and damping of the foundation resulting from the specific load condition. An example of a calculated 1 h storm load history of a monopile foundation for a 10 MW turbine is shown in Fig. 1. The figure shows the resulting overturning moment at the seabed. For this load case (idling situation), the maximum overturning moment is about 400 MNm and the average load is low compared to the cyclic load, resulting in a close to two-way cyclic loading condition of the foundation. Looking in detail around the peak load, it is seen that the peak load is applied within a period of around 1 s. In addition, the cyclic load period is about 4 s, which is in line with the natural period of the OWT considered. This load history will be used in the evaluation of the explicit soil models.

BEHAVIOUR OF SATURATED DENSE SAND
The deformation characteristics of sand are rather complex and depend on several index properties of the sand (e.g. grain size distribution, fines content, grain shape and angularity, mineralogy, etc.), void ratio, effective mean stress and stress path (compared to the actual stress state and the fabric of the grain skeleton), and the effect of stress and/or strain histories (e.g. maximum mean effective stress, cyclic loading, etc.). This was shown for instance by Wichtmann (2016) and Andersen (2015). An extensive laboratory test programme on a given sand is therefore required to determine the optimised sizes of an OWT foundation.
An extensive test programme was recently performed on the Dogger Bank sand (Blaker & Andersen, 2015. The Dogger Bank site is located in a shallow area of the North Sea between 125 and 290 km off the east coast of Yorkshire, UK, with water depth from 18 to 63 m. Three offshore wind farms are planned at this location. Figure 2 presents the grain size distribution of two mixed batches of Dogger Bank sand: a clean sand (batch A) and a sand with approximately 20% silt content (batch B). Only the batch A sand is considered in this paper; however, the silt content will affect the properties of the sand (e.g. reduced permeability). The sand is sub-angular to sub-rounded according to the classification by Russel & Taylor (1937). The grains are composed of 91 to 93% quartz, 2-4% carbonate and a small amount of other minerals. The maximum and minimum void ratios of batch A were e max = 0·865 and e min = 0·597, with a unit weight of the solid particles of γ s = 26·3 kN/m 3 . Fig. 2 also shows the grain size distribution of a rather similar North Sea sand, from the Siri field, also discussed later in this paper. All laboratory tests were performed on fully saturated specimens prepared by moist tamping to the target relative densities as described in Blaker & Andersen (2019). Figure 3 shows two oedometer curves, each including two unloading/reloading loops starting at different vertical effective stresses. These curves show how the oedometer stiffness increases with increasing effective vertical stress and increasing relative density. Furthermore, the initial unloading stiffness is significantly higher than the virgin loading stiffness at the same vertical effective stress. The reloading stiffness is slightly less stress-dependent than the unloading stiffness. During drained cyclic loading, the different behaviour of the sand during unloading and reloading may be one reason for the accumulated volumetric strain. The measured tangential oedometer stiffnesses can be described very accurately by the equations proposed in Andersen & Schjetne (2012).
The maximum shear modulus G max was measured by piezo-ceramic bender elements on direct simple shear (DSS) specimens before and after pre-shearing for different vertical stresses and relative densities. The measured values before pre-shearing agreed well with the equation of Hardin & Drnevich (1972); however, after pre-shearing the values were reduced, giving a minimum value for all tests of around 50% of the calculated value.
The stiffness and strength characteristics were determined from anisotropically consolidated drained and undrained triaxial compression and extension tests, starting at an axial effective stress σ′ a0 = 40 and 200 kPa with the coefficient of earth pressure at rest K′ 0 of 0·45 and a relative density D r of about 80 and 100%.
Figure 4(a) shows the shear stress τ = (σ a À σ r )/2 normalised with the effective consolidation stress σ′ a0 plotted against the axial strain ε a . In the undrained compression tests (TUC and TUE tests in the figure, where 'C' is for compression and 'E' is for extension), the shear stress ratio τ/σ′ a0 increases more or less linearly with axial strain ε a . Failure in these cases (as the sand dilates) is limited by the cavitation pressure in the pore water. From the drained compression tests (TDC and TDE tests in the figure), the peak strength increases with increasing D r and decreasing σ′ a0 . The peak friction values, which are mobilised at ε a = 0·3-0·5%, follow the same trends with respect to σ′ a0 and D r as shown in Andersen & Schjetne (2012) for a large number of different sands, but the strength values for the Dogger Bank sand are higher. After peak, the strength reduces with increasing strain. However, owing to non-uniform deformations (strain localisation), this part of the curves is not reliable. Fig. 4(b) shows the development of volumetric (dilation) strain ε vol plotted against ε a from the drained tests. These curves define the dilation properties of the sand as a function of void ratio, mean effective stress and shear strain. The dilation in extension is significantly smaller than the dilation in compression. Figure 5 highlights the different behaviour during undrained compression, undrained extension and constant-  volume DSS testing. The response is significantly stiffer in triaxial compression than in triaxial extension and DSS for a shear strain larger than 0·5%. The DSS specimen has the lowest stiffness. However, the stress paths considered are not necessarily representative of the response of the soil around a monopile. It is therefore important to have a model that can extrapolate the behaviour from the laboratory data to the behaviour under the actual stress paths around the monopile.
Based on the cyclic undrained triaxial and DSS tests with varying combinations of cyclic and average shear stress, the response of the soil (e.g. the accumulated pore pressure u acc and the accumulated shear strain γ acc as a function of the number of cycles) can be presented in so-called cyclic contour diagrams (Andersen, 2015). Examples of crosssections through these contour diagrams are shown in Fig. 6.
Additional cross-sections for the Dogger Bank sand are presented in Blaker & Andersen (2019). Figure 7(a) shows the effective stress paths during the first two cycles in a p′-q plot (where q = 2τ) for normalised cyclic shear stress ratios τ cy /σ′ a0 of 0·25, 0·5 and 0·75. Fig. 7(a) shows that there exists a line of maximum friction with slopes that are different in compression and extension. Fig. 7(a) also shows that the soil contracts immediately after stress reversal. The slope of the unloading stress paths becomes more inclined (less steep) with increasing dilation along the maximum mobilised friction lines. The unloading stress paths after stress reversal on the compression side are more or less linear until they reach the phase transformation line (PTL), where dp′ = 0, on the extension side. The PTL is unique on the extension side, while the slope is increasing with increasing amplitude on the compression side. The reason for this behaviour is not fully understood. According to the authors, it is unclear whether any models can account for this behaviour. The largest accumulated pore pressure u acc occurs during the first cycle. However, Δu acc is very small compared to the maximum change in p′ within the cycle. Therefore, the constitutive model needs to describe these stress paths very accurately in order to predict u acc properly for all combinations of average and cyclic shear stresses and number of cycles. Figure 7(b) shows the corresponding τ-γ curves. For τ cy /σ′ a0 = 0·25, the soil is significantly softer during the first quarter of a cycle. Thereafter, the secant cyclic stiffness is close to constant from one cycle to the other. For τ cy /σ′ a0 = 0·5, the stiffness is softest during dilation on the extension side, and the behaviour is significantly stiffer during reloading. For τ cy /σ′ a0 = 0·75, the characteristic 'boomerang' shape starts to appear, where the stiffness is lowest slightly after crossing the horizontal axis during reloading. The main reason for this low stiffness is seen from the p′-q plot, where p′ is lowest close to PTL on the compression side. This demonstrates the importance of capturing the effective stress path properly in order to be able to model the change in stiffness properly. In addition, the stiffness during dilation on the compression side is significantly higher than during dilation on the extension side. From laboratory tests on the Siri sand, it was found that the strain and pore pressure accumulation for the dense sand subjected to the same cyclic load history may show a wide spread in the results (Fig. 8), while the response within the undrained cycle is nearly uniquely defined by the effective stresses at the beginning of the cycle regardless of the previous load history (Fig. 9).
The rate of strain accumulation can also be found from cyclic drained tests. For instance, an extensive laboratory test programme of cyclic drained tests on sand was presented by Wichtmann & Triantafyllidis (2016a, 2016b. In cyclic drained tests, the accumulated strain increases with increasing number of cycles, while the rate of strain accumulation is decreasing. The advantage of cyclic drained tests is that the effective stress path is the same from cycle to cycle. It is then easier to isolate the effect of earlier stress or strain history in the measured response. The disadvantage is that it is difficult to apply as large cyclic shear stresses in drained tests as in undrained tests.
As an alternative to fully drained or fully undrained cyclic tests, a given package of undrained cycles could be applied several times with partial drainage and dissipation of the excess pore pressure between each package. This condition is more realistic for the actual behaviour around a foundation in sand with drainage during several cycles. This kind of test is exemplified in Fig. 10. The test shows the same trend as the drained cyclic teststhat is that the rate of pore pressure accumulation decreases with increasing volumetric strain, even though the volumetric strain is so small that it hardly changes the relative density.

HS-Small
The hardening soil (HS) model is an implicit model described in Schanz et al. (1999) and in the Plaxis manual . This model has been used to back-calculate model tests and to calibrate p-y springs for the PISA design method (Taborda et al., 2019). The PISA design method is implemented in the design tool Plaxis MoDeTo (Panagoulias et al., 2018). Key features of the HS model are listed below.
(a) Elasto-plastic formulation with isotropic hardening of two yield surfaces, a cone and a cap. (b) The expansion of the cone (mobilised friction) gives increased plastic shear strain following a hyperbolic hardening function. (c) The cap defines the pre-consolidation stress for a general three-dimensional (3D) stress state and the difference in bulk modulus during virgin loading and unloading/reloading. The same stress dependency (expressed by the same exponent) is used for the virgin compression, unloading and reloading stress paths. (d ) Dilatancy (the ratio between plastic volumetric strain increment and plastic shear strain increment) is based on Rowe's formulation (Rowe, 1962). (e) There is a dilatancy cut-off at a specified critical void ratio. ( f ) There is a tensile effective stress cut-off.
In addition, the model may be expanded by a small-strain stiffness degradation formulation as described by Benz (2007) and in the Plaxis manual . The hysteretic behaviour during unloading/reloading follows Masing's rule. The model requires a total of 13 input parameters. These parameters are listed in Table 1.
Interpretation of input parameters and back-calculations of laboratory tests. The material parameters were calibrated for the Dogger Bank sand with D r of 80%. An independent calibration is needed for D r of 100%. Without additional tests on other densities and since it is only the small-strain shear modulus that is a function of the void ratio, material properties for other densities need to be interpolated and extrapolated based on these two data sets.
Since HS-Small is based on isotropic hardening, the model can describe only small-strain cycles around a given effective stress level without strain accumulation. Therefore, only the monotonic tests were analysed. The effect of cyclic loading would need to be accounted for by another model.   The material properties were first calibrated by backcalculation of the oedometer test in Fig. 3. The virgin loading curve was duplicated with a stress exponent varying between m = 0·5 and 0·6. The unloading and reloading curves were more difficult to fit because the same m is used for virgin loading, unloading and reloading. For HS-Small, the unloading/reloading stiffness needs to be about 3 times the virgin loading stiffness at the same effective stress due to the mathematical formulation of the model. This means that a good fit can be obtained for a selected stress range only. With E oed,ref = 55 MPa, E ur,ref /E oed,ref = 3 and m = 0·5, a good fit is for instance obtained for an unloading/reloading range between σ′ a = 200 and 60 kPa.
The model parameters were also calibrated with the undrained triaxial tests. When using the parameters from the oedometer tests, it was not possible to fit the undrained triaxial tests: the slope of the undrained shear stress-shear strain curve was much lower than the experimental data. To overcome this problem, the pre-consolidation stress was set to an arbitrary value large enough not to be exceeded. A physical explanation for this may be that the pre-consolidation stress should not be updated without a change in volumetric strains. Furthermore, to obtain an improved fit of the non-linear undrained stress-strain curves, m had to be reduced from 0·5 to 0·35.
To fit the effective stress paths in the undrained triaxial tests, a drained friction angle of 43·6°was selected. However, it was not possible to fit the stress-dependent peak friction angle from the drained tests, or the dilatancy angle. The dilatancy angle used to fit the undrained shear stress-strain curves was found to be different in compression and extension (Figs 11 and 12, ψ = 30°and 12°), and even lower in DSS (ψ = 8°). Furthermore, HS-Small was not able to capture the initial contraction part in the extension tests, thus leading to deviation in the stress-strain curve at low strains.
It was therefore not possible to model adequately the drained and undrained compression, extension and DSS stress paths with the same set of parameters. Calibrating the model parameters based on the triaxial compression tests will give stiffness values that are too high for a general stress path.  To calibrate the small-strain stiffness for HS-Small, the shear stress-shear strain curves up to 0·1% shear strain were considered. A reference initial shear modulus G 0 ref of 200 MPa (minimum principal stress σ′ 3 of 100 kPa and shear strain of γ 0·7 of 0·02% at 28% reduction of the initial stiffness) was selected to fit this part of the curves. The results are shown in Fig. 13. The selected G max is then stiffer than obtained by the equation of Hardin & Drnevich (1972) that is G 0 ref = 115 MPa, and thus also stiffer than measured in the DSS apparatus.
To fit the stress-strain curves from the drained compression tests, the dilatancy angle ψ must be lower than that used to fit the undrained compression tests. The two different parameter sets that gave the best fit for the undrained and drained triaxial tests of Dogger Bank sand with D r of 80% are presented in Table 1. HS-Small is not able to describe the post-peak softening response observed in the drained tests and cannot be used to predict the effect of cyclic loading of an OWT foundation.

Sanisand
To overcome several of the limitations of HS-Small, a more advanced soil model should be used. There exists a family of models (e.g. Li & Dafalias, 2000;Dafalias & Manzari, 2004;Taiebat & Dafalias, 2008 etc.) based on the framework described in Manzari & Dafalias (1997). These models are often called Sanisand (simple anisotropic sand). The formulation from Dafalias & Manzari (2004) was used in this paper. It includes an effect of fabric change where the contraction after load reversal is enhanced by the effect of dilation before load reversal (as seen from the tests in Fig. 7).
A user-defined UMAT-subroutine in the finite-element product suite Abaqus (Abaqus, 2014), together with a Plaxis interface, are available at the SoilModels website (SoilModels, 2020). In this study, the authors' implementation of the model in Plaxis was used, which may give slightly different results than the UMAT-subroutine.
The key features of the implicit Sanisand model and its main input parameters are summarised below. The stress ratio (q/p′) based elasto-plastic bounding surface model is illustrated in Fig. 14(a), and the simplified equations for the triaxial stress state are summarised in Fig. 14(b).
(a) Mean effective stress-dependent critical state formulation (void ratio e and stress-path-dependent residual strength, M ). (b) Elastic shear modulus G defined by a dimensionless input parameter, G 0 , together with an equation dependent on p′ and e (Fig. 14(b)). (c) Elastic bulk modulus K given by the pressure-dependent G and a constant Poisson's ratio, ν.  (Fig. 14(b)). (e) The model is without a cap to control the difference in compressibility during virgin loading, unloading and reloading. A model with a cap surface is proposed in Taiebat & Dafalias (2008). ( f ) The plastic modulus depends on the distance from the current stress state to the bounding surface, the distance from the previous stress state of reversal, and is expressed with an equation of e, p′ and an input parameter h 0 . This means that the failure strain may vary both with e and p′, as observed in the tests in Fig. 4. (g) The dilatancy (ratio between plastic volumetric strain and plastic deviatoric strain) depends on the distance to the phase transformation surface (PTS) and an input parameter A 0 . The soil contracts when inside the PTS and dilates outside. The amount of contraction can be enhanced by an equation for fabric change due to dilation before stress reversal. (h) A PTS with inclination M d in the p′-q space depends also on ψ = e À e c and an input parameter n d . This means that the PTS may expand or shrink depending on the sign of ψ = e À e c . (i) The model has a small elastic regime (the size of this region is defined by the input parameter m) that moves together with the change in mobilised stress ratio, η = q/p′.
Interpretation of input parameters based on back-calculations of laboratory tests. The variation in the hardening parameter H and dilatancy parameter d can be calculated from the stresses and strains from the laboratory tests. The input parameters can then be determined by fitting the Sanisand functions of H and d to these curves. This type of interpretation was presented in Dahl et al. (2018). However, for an incremental model where the strain is obtained from the integration of incremental strains, this only works if the curves of H and d fit reasonably well with the equations used in the model, otherwise the accumulated error becomes too large.
In this paper, the input parameters were determined by trial and error to ascertain the effect of uncertainties in the selected parameters. The elastic parameters were initially estimated from the measured G max in the DSS apparatus (initial value of G 0 = 200 and ν = 0·05).
Since sands with relative densities D r of 80 and 100% only were tested, one challenge was to determine the p′-dependent critical void ratio e c on the critical state line (defined by three input parameters e c0 , λ c and ξ in Fig. 14(b)). Without tests on looser specimens, e c0 , λ c and ξ together with n b were simply calibrated to obtain a good fit with the measured drained peak strength M b for the two different stress levels (σ′ a0 = 40 and 200 kPa) and the two relative densities (D r = 80 and 100%). Fig. 15 shows the selected e c line where e c0 $e max = 0·87, λ c = 0·014, ξ = 1·0, n b = 1·33 and M c = 1·45 (i.e. ϕ c = 36°). However, other combinations, for instance starting with M c = 1·29 (ϕ c = 32°), also give a reasonably good fit to the measured drained peak strength M b with n b = 1·85-1·9. It was not possible to fit the post-peak softening response. This is most likely due to the non-uniform deformations observed from an examination of the specimens after the tests. For both the drained and undrained extension tests on the specimens with D r = 80%, a good fit was obtained with M e = 0·86 (ϕ c = 30°). However, this gave c = M e /M c = 0·86/1·45 = 0·59, which is a rather low value. One explanation could be that the peak extension strength is reached before failure due to necking. A ratio M e /M c of 0·59 may cause numerical problems, depending on the implementation of the model in the finite-element code.
The hardening parameter h 0 = 9 was found by fitting the calculated curves of the stress ratio q/p′ plotted against the deviatoric strain with the corresponding curves from the four drained compression tests.

EVALUATION OF SOIL MODELS FOR DESIGN OF OWT FOUNDATIONS
A good fit to the measured curves of volumetric strain plotted against axial strain (shown in Fig. 4(b)) was obtained with A 0 = 1·0 and n d = 3·0. However, a good fit can also be obtained with other combinations of these two parameters. To obtain a unique set of parameters, a well-defined phase transformation line is needed.
An improved fit with the drained tests and the initial part of the undrained tests could be obtained by increasing G 0 from 200 to 250. However, the stiffnesses for the undrained tests after about 1% shear strain then became too high.
The calculated responses using Sanisand with the set of monotonic material parameters in Table 2 fit reasonably well with all the measured responses from Fig. 4(a), as shown in Fig. 16.
Once a good fit of the static tests has been obtained, the model should ideally also be able to describe the soil response during complex cyclic load histories, as for instance shown in Fig. 1. Fig. 17 compares the measured laboratory response from Fig. 7 with the calculated response with the calibrated parameter set from the monotonic tests and with the set of adjusted parameters in Table 2 (denoted Sanisand-improved cyclic fit). The following observations are made: (a) the calculated cyclic shear strains and hysteretic damping with the parameter set from the monotonic tests were too large, except for the largest stress cycle; (b) the behaviour was improved significantly by increasing h 0 from 9 to 35; (c) the stress paths were slightly improved by using the equation for fabric change with the parameters in Table 2; (d) the model was not able to capture the high contraction immediately after stress reversal for the largest stress cycle, especially on the extension side; (e) the calculated accumulated shear strain and pore pressure did not agree with the measured response.
The present version of Sanisand is therefore not recommended to predict the effect of cyclic loading of an OWT foundation in dense sand. It would thus be interesting to use the same data set to calibrate the model including a memory-enhanced bounding surface proposed by Liu et al. (2019), which should improve the modelling of cyclic strain accumulation.

High-cycle accumulation model
Based on the evaluation above, neither HS-Small nor Sanisand appear to be able to predict the cyclic behaviour of a dense sand. Therefore, the high-cycle accumulation model (HCAM), as described in Niemunis et al. (2005), was considered. In this explicit model, the accumulated strain ε acc is calculated as a function of the number of cycles. The cyclic load history was idealised with packages of different constant cyclic load amplitudes that can increase or decrease from package to package.
In the present model, the changes in effective stresses were given by where K = Bp atm ( p av /p atm ) m , G = 1·5K(1 À 2ν)/(1 + ν), p atm = 100 kPa and the constants B, m and ν are input parameters. The quantity dε vol is the change in volumetric strain; dε q is the change in deviatoric strain; and m v and m q are the volumetric and deviatoric components of the normal vector to the modified Cam-Clay yield surface passing through the current stress point. The slope of the critical state line is used to fit the measured ratio between dε q and dε vol . The change in accumulated strain dε acc is defined by equation (3).
In the authors' user-defined model in Plaxis, plastic strains due to changes in effective stresses were not included, because    However, the performance of the model would most likely have been better if plastic strains due to changes in effective stresses had been included.
The key assumption in this model is that the rate of strain accumulation dε acc is given by the current average effective stresses, void ratio and an equivalent number of drained cycles N eqv of strain amplitude ε ampl . The rate of strain accumulated is given by a set of empirical functions where f ampl gives the effect of the cyclic strain amplitude; f N accounts for the cyclic history by the equivalent number of drained cycles; f e accounts for the effect of void ratio; f p accounts for the effect of effective mean stress level; f Y accounts for the effect of the average shear stress level; and f π accounts for the effect of change polarisation of the cycles. The actual functions used in HCAM are given Table 3.
In the model implementation, the strain amplitude ε ampl was, for simplicity, calculated by the following degradation function (Hardin & Drnevich, 1972) where γ cy is the cyclic shear strain; γ r is the reference shear strain; and G max is given by the equation of Hardin & Black (1966) where exponent α, constants g 0 , e g and n g are input parameters to the model. The γ cy is independent of the average shear stress and the cyclic load history. However, for undrained and partially drained conditions, p′ is reduced due to pore pressure accumulation. A more advanced model, as for instance Sanisand, could also be used to calculate γ cy .
Calibration of model parameters. The determination of the parameters to HCAM is usually based on a tailor-made set of drained cyclic triaxial tests with a large enough range of the governing variables (Wichtmann & Triantafyllidis, 2016a, 2016b. To demonstrate the calibration process, the parameters fitting the drained cyclic triaxial tests in Fig. 18 were determined. The test has an average effective stress state (σ′ ac = 200 kPa and K 0 = 0·45), D r of 80%, and a nearly constant cyclic strain amplitude ε ampl = 8 Â 10 À4 . Only the functionḟ N can be calibrated from this test. The product of the other functions was set to a constant C f = f ampl f e f p f Y f π . A very good fit was obtained with the parameters C f = 8·3, C N1 = 6 Â 10 À5 , C N2 = 1·4, C N3 = 1·0 Â 10 À4 and ϕ c = 25·6°a s illustrated in Fig. 18. The model was then calibrated to estimate the accumulated pore pressure around a monopile in dense sand during a ULS load history. The accumulated pore pressure during cyclic loading around the K 0 -stress state shown in Fig. 6(a) was used. Fig. 19 shows the accumulated pore pressure u acc during the first ten cycles for τ cy /σ′ a0 = 0·25, 0·32, 0·5 and 0·75.
The calibration of the HCAM parameters included the following steps (the model parameters are listed in Table 4).
(a) The model parameters that could not be obtained from the laboratory tests were estimated from the database presented in the paper by Wichtmann et al. (2015). Table 3. Functions in HCAM used in analyses (Niemunis et al., 2005) Influencing parameter Function Material constants Average mean pressure These are the parameters that account for the change in void ratio e and different average shear stress mobilisations (C e = 0·58 and C Y = 2·5). The parameter that accounts for different effective mean stresses was set to C p = 0 for simplicity. Cyclic tests consolidated to different p′ 0 would be required to estimate this parameter. However, the model accounts for different effective stress levels by G max and K that are dependent on p′. (b) K was estimated from the unloading/reloading loops in the oedometer test. A value of K of 235 MPa was used, obtained by B = 1162, ν = 0·25 and m = 0·5. (c) The measured γ cy during the first cycle (N = 1) from Fig. 6(b) was fitted by the cyclic degradation model (equations (4) and (5)) using g 0 = 300, γ r = 0·012% and α = 0·8 (as shown in Fig. 20).
A reasonably good fit to the measured accumulated pore pressure u acc was obtained using the parameters in Table 4 as illustrated in Fig. 19. However, fitting the accumulated shear strain γ acc from Fig. 6(c) with the same set of model parameters was not possible. Furthermore, HCAM was not able to back-calculate the response in Fig. 10 with partial drainage between the cycles. The significant u acc during the first cycle after pore pressure dissipation was not predicted. The calculated Δu acc in this cycle was even lower than in the last cycle before pore pressure dissipation due to increased N eqv given by the contribution of ε vol to ε acc and lower γ cy due to the increased p′. It may be that the first cycle should be calculated by another model, as recommended in Niemunis et al. (2005). However, this model would need to account for the previous history not explained by the small change in e from package to package (D r is only changing from 78·3 to 79·7% during the whole test). Based on this, one should be careful by using the present version of HCAM for calculations of accumulated pore pressure during ULS conditions. However, the model may be suitable for calculation of accumulated strains in cases with significant drainage during cyclic loading, for instance to predict the accumulated tilt of a monopile during the lifetime of the OWT.
Simulation of undrained and partially drained cyclic triaxial tests. The above parameter set was used to analyse the development of the accumulated pore pressure u acc during the ULS load history in Fig. 1. The idealised history of packages of constant cyclic load within each package (Table 5) was established by the procedure in Norén-Cosgriff et al. (2015). For simplicity, only the six largest load packages were modelled, starting from a load level that is 55% of the maximum load. The maximum cyclic shear stress applied was τ cy /σ′ a0 = 0·63 with σ′ a0 = 200 kPa. The cyclic load period T p was 4 s. The load history was applied τ cy /σ' a0 = 0·32 τ cy /σ' a0 = 0·5 τ cy /σ' a0 = 0·75 HCAM Lab Fig. 19. Accumulated pore pressure during first ten cycles from cyclic undrained triaxial tests, stress amplitudes of τ cy /σ′ a0 = 0·25, 0·3, 0·5 and 0·75 for σ′ ao = 200 kPa and D r = 80%, and as predicted by HCAM  for the simulation of both an undrained cyclic triaxial test and a partially drained test where the flow resistance of the filter at the top of the specimen was adjusted to represent a 9 m long drainage distance through the sand with a permeability k of 5 Â 10 À6 m/s. Figure 21 shows the calculated u acc /p′ 0 with the number of cycles for the two cases. For the partially drained case, u acc /p′ 0 is about 0·1 at the end of the history. This result is compared with a similar analysis using the partially drained cyclic accumulation model (PDCAM) in the next section.
Partially drained cyclic accumulation model The PDCAM was described in the paper by Jostad et al. (2015). The laboratory test programme carried out on the Dogger Bank sand was tailor made for this model, since the main input to the model is based on cyclic contour diagrams ( Fig. 6 and Blaker & Andersen, 2019). Compared to HCAM, one avoids the rather cumbersome calibration process to fit the laboratory data, since these data (u acc , γ cy and γ acc ) are already included in the contour diagrams as a function of τ cy /σ′ a0 , τ a /σ′ a0 and number of cycles N.
The key features and assumptions of the explicit PDCAM model include the following.
(a) The stresses are divided into average effective stresses and undrained cyclic stresses. Similarly to HCAM, cyclic and average calculation phases need to be coupled. (b) The incremental response within individual cycles is not modelled. Similarly to HCAM only τ cy (or γ cy ) needs to be calculated, γ acc and u acc are then given by the number of cycles of this amplitude. (c) Δp = p′ 0 Àp′ is a state parameter that accounts for the load history.
In addition to the cyclic contour diagrams, the inputs to the model are parameters for the stress and stress-pathdependent bulk modulus defined by the equations in the paper by Andersen & Schjetne (2012) and the permeability.
The results from the monotonic drained triaxial compression and extension tests are directly included in these diagrams, however, without the post-peak strain softening part. The extension from the triaxial state to the general 3D state is done by assuming coaxiallity between principal deviatoric strains and stresses.
A key assumption in PDCAM is that the response is governed by the current average effective stress state, relative density and the applied cyclic shear stress. The actual stress historythat is, the manner in which average effective stress state is reacheddoes not matter. This assumption agrees with the result shown in Fig. 9. Therefore, the effect of volumetric strain due to the dissipation of u acc has no effect on the rate of strain and pore pressure accumulation as long as D r remains close to unchanged. However, the partly drained cyclic laboratory test results in Fig. 10 show that it has an effect. The present version of the model is therefore not recommended to model conditions with small changes in p′ or to calculate the long-term accumulated displacements and tilt of OWT foundations.
Simulation of undrained and partially drained cyclic triaxial tests. To illustrate the performance of PDCAM, the same ULS load history, from the undrained and partially drained cyclic triaxial tests simulated by HCAM, was analysed with PDCAM. For the undrained condition, PDCAM perfectly fits the accumulated pore pressure against number of cycles for the tests shown in Fig. 19 since the curves are established from the contour diagram in Fig. 6.
The calculated u acc /p′ 0 plotted against number of cycles for the two drainage conditions are compared with the results obtained by HCAM in Fig. 21. For the undrained condition, the calculated u acc values with PDCAM and HCAM are somewhat different due to different procedures to account for the cyclic historythat is the equivalent number of undrained cycles in PDCAM compared to the equivalent number of drained cycles in HCAM. For the partially drained case, both models show that u acc /p′ 0 is around 0·1 at the end of the history. Based on this evaluation, PDCAM is expected to be suitable for calculating the displacements and capacity of an OWT foundation in sand during an ULS storm condition.

FINITE-ELEMENT ANALYSES
To further study the performance of the four material models, simplified two-dimensional (2D) analyses of a vertical cross-section through a monopile in dense sand Table 5. Idealised one hour ULS storm load history from the time history in Fig. 1 Amplitude group no. were carried out. The main results from these 2D analyses are then compared with the results of 3D analyses of the same monopile. The finite-element program Plaxis  was used.

Simplified finite-element analysis model of a monopile
The main purpose of the simulations was to compare the performance of the four models for stress conditions that were somewhat representative for an OWT monopile, and not to calculate accurately the response of the monopile. This would require a 3D model and accounting for installation effects, particularly for early load cycles.
The monopile has a diameter D of 6 m and a penetration depth into the sand of 30 m. The stiffness of the pile was calculated based on a steel wall thickness of 0·08 m. The ratio between the overturning moment and horizontal shear force at the seabed was 27 m, while the water depth was 25 m. The monopile was for simplicity assumed to be wished in place without any effects of the installation. The 2D model of the ground was 160 m wide and 80 m deep below the seabed and contained 1000 15-noded triangular elements with refined mesh around the monopile. Interface elements with the possibility of a tension cut-off were used between the monopile and the soil. The bottom boundary was fixed, the horizontal displacements were fixed and the vertical displacements were free along the vertical boundaries.
Effect of partial drainage during monotonic loading. A maximum reference horizontal load of 15 MN/m was applied at the monopile 27 m above the seabed. This load was then gradually scaled by an increasing load factor. Both a fully drained and a perfectly undrained case were analysed. For the partially drained case, the reference load was increased within a time period T of 1 s (Fig. 1). A loading period of 40 s was also used to illustrate the effect of loading rate. For HS-Small and Sanisand the material properties given in Tables 1 and 2 were used. For HS-Small, the triaxial compression parameter set was used at the back of the monopile and the triaxial extension set was used at the front to account for the stress-path-dependent behaviour of the sand. For Sanisand the material properties based on the monotonic tests were used because the load history prior to the peak load was neglected. A permeability of the sand of k = 5 Â 10 À6 m/s was used in the analyses. It was found that scaling k had the same effect as scaling T. Figure 22 presents the calculated horizontal displacement (u x ) of the monopile at the seabed plotted against the load factor for the 2D finite-element analyses (FEAs) with HS-Small and Sanisand. The modelled stiffness and capacity of the undrained cases were significantly larger than the drained cases. For the undrained case, there was a monotonically increasing capacity, which agrees with the undrained monotonic laboratory tests. At u x = 0·6 m (10% of D), the undrained capacity was about 3 times the drained capacity. The response for the partially drained case for T = 1 s was very close to undrained. For T = 40 s, the capacity was about 25% lower than the undrained case, but still significantly larger than the drained case. When accounting for cavitation with a cut-off in the pore water according to the actual water depth or for a tension crack with a cut-off at the back of the monopile, the undrained capacity and stiffness were reduced significantly and were closer to the drained capacity.
Sanisand gave significantly softer response than HS-Small, for both undrained and drained conditions. The main reason for this difference is that the calibrated small-strain stiffness in HS-Small, based on the initial part of the triaxial test ( Fig. 13), was stiffer than the initial stiffness in Sanisand, based on the measured shear wave velocity on the DSS specimens. It was also difficult to increase the initial stiffness in Sanisand without affecting the performance of the model at larger strains. This limitation is for instance improved in the paper by Taborda et al. (2014).
To demonstrate that the conclusions derived from the 2D analyses were reasonably representative of a 3D response of a monopile, some of the analyses were repeated with a 3D FEA Plaxis model, using HS-Small. The 3D FEA model was 160 m wide and 50 m deep below the seabed, with symmetry in the x-z plane. The mesh consisted of 38 000 ten-noded tetrahedral elements with refinements close to the monopile. The minimum element size was 0·5 m. The loading conditions were the same as for the 2D model. The average of the triaxial compression and extension parameters in Table 1 was used. The 3D results are compared in Fig. 23. One analysis under undrained conditions was done using the soil properties calibrated only with the DSS tests. The same trends as illustrated for the simplified 2D analyses (Fig. 22) were found with the 3D FEAs. However, the amount of drainage is slightly greater, as expected due to dissipation of the pore pressure normal to the loading direction. The results from the undrained analyses using the DSS properties were surprisingly very similar to the calculated drained response using the average properties from triaxial compression and extension.
Effect of cyclic loading. The idealised ULS cyclic load history (Table 5) considered in the partially drained cyclic triaxial tests was also applied to the 2D model of the monopile, using both HCAM and PDCAM. The maximum cyclic overturning moment at the seabed in the last package was 372 MNm (based on Fig. 1) over D = 6 m. The average loads were set to zero. For this load condition there will be no accumulated displacements of the monopile. The cyclic shear strains at the integration points were calculated with the non-linear cyclic model in HCAM and with the cyclic contour diagrams in PDCAM assuming fully undrained conditions. The number of cycles in the load packages was given by the time divided by the cyclic period T p = 4 s. The cyclic load history was applied with a constant permeability of the sand of k = 5 Â 10 À6 m/s. Drainage was enabled at the seabed and at the outer vertical boundary. The HCAM and PDCAM analyses indicated that u acc during the load history was very small. Fig. 24 shows the relationship of Δp′ = ( p′ 0 Àp′) with number of cycles at a point 10 m beneath the seabed close to the monopile. The results show that instead of pore pressure accumulation due to cyclic loading, the mean effective stress p′ was reduced due to stress relaxation caused by the limitation of the soil to deform according to the accumulated volumetric strain. During cyclic loading, the reduction in p′ around the monopile gradually propagated downwards. At the end of the storm, p′ was reduced significantly in the upper 6 m of the monopile over a distance of roughly 2 m outside the pile. Cyclic loading therefore has an effect on the stiffness and displacements of the monopile during an ULS storm. However, the effect depends on the monopile dimensions, soil conditions, water depth, and the tower-and turbine-dependent load conditions. PDCAM could also be used in a 3D FEA to predict this effect as done with the 2D FEA simulations. HCAM could also be used if the properties are calibrated for the actual stress level around the monopile, which for a ULS storm condition may be larger than typically used in drained cyclic laboratory tests.

CONCLUSIONS
In this paper the complex behaviour of a dense North Sea sand subjected to combined average and cyclic loading has been discussed. The results of the analyses demonstrate the requirements, parameters, advantages and limitations of four different soil models that may be used to analyse the stiffness, capacity and displacements of OWT foundations in sand. Even though FEAs are rarely used directly in the design of OWTs, they have become increasingly more popular in the process of calibrating, for instance, the p-y curves that are used in the design.
The response of a clean, dense North Sea sand with a permeability k of 5 Â 10 À6 m/s, surrounding a monopile with a diameter of 6 m, is found to be close to undrained during the application of a large cyclic load within a period of only a few seconds. The triaxial compression stiffness and shear strength of this sand are significantly higher under undrained conditions than under drained conditions. However, the high strength may be limited by the cavitation pressure in the pore water. In addition, a gap may develop on the tension side of the monopile.
The analyses further demonstrate the importance of an appropriate description of the small-strain stiffness of the sand in order to calculate the displacements of a monopile, not only at low load levels, but for the entire range of loads. The anisotropic and stress-path-dependent behaviour of the sand also needs to be taken into account, since the undrained stiffness in a DSS state of stress is significantly lower than in a triaxial compression state of stress.
The effects of several large cyclic loads prior to the maximum ULS load may, in addition to causing pore pressure accumulation, reduce the mean effective stress around the monopile, due to stress relaxation. Explicit models such as PDCAM and HCAM, which may account for accumulated strains in a proper way, can be used to estimate the effect of cyclic loading on the displacements and capacity. However, this type of model needs extensive tailor-made laboratory tests to determine the parameters required as input. Explicit models may also be used to predict the accumulated displacements and tilt of an OWT foundation during its lifetime. It is then important to account for the effect of reduced rate of strain accumulation with increasing number of cycles as observed during drained cyclic tests.
As the measured strain accumulation within a cycle is very small compared to the maximum variation of the strain within a cycle, it is generally very difficult to describe the development of accumulated strain and pore pressure with an implicit model. This may be solved by adding more input parameters in the implicit models to control the rate of strain accumulation. Another alternative to overcome this problem may be to combine a good implicit model with a proper explicit formulation of strain accumulation.
It should be noted that the mechanism of strain accumulation in dense sand is not fully understood. In particular, the effect of a varying load history under partially drained conditions on the rate of strain accumulation requires more research. In addition, more representative partially drained cyclic tests of monopiles are needed to verify advanced numerical models and to ensure that design of OWT foundations can be optimised. Finally, the required dimensions of an OWT monopile and the corresponding governing load condition are currently not well defined.