A multi-spring model for monopile analysis in soft clays

The offshore wind industry in China has seen a rapid development in recent years and is projected to account for half of global yearly installed capacity in the years to come. Monopiles are a popular foundation solution for supporting offshore wind turbines. However, due to challenging seabed soil conditions, which often feature thick normally consolidated soft clays and harsh environmental loading (i.e. frequent occurrence of typhoons), extremely long monopiles are often designed. The large monopiles are costly to fabricate and install and sometimes kills the viability of the concept in cases where the bedrock is relatively shallow and expensive piling in rock is otherwise required. However, the state-of-practice for designing these monopiles in clays is typically by using the API p-y springs, which are widely known to underestimate soil-pile interaction stiffness and capacity for large diameter monopiles. Improvement to the design method can therefore have significant economic implications to the industry. The paper presents an effort toward this direction. A multi-spring beam-column model suitable for monopile design in soil conditions in China is proposed. The model features three soil spring components, namely the lateral p-y spring, the pile tip base shear s-y spring and rotational m-θ springs along the pile shaft. The validity of the model is verified by a comprehensive suite finite element parametric analyses. Guidance for incorporating the cyclic loading effect into design is also provided. The model proposed in this paper has large potentials for application in design practice.


Introduction
The offshore wind industry in China has seen a rapid development in recent years. The accumulated installed capacity had reached 6.7 GW by the end of 2019 and this will be added with over 2-3 GW annually in the years to come. The foundations supporting the offshore wind turbines in China take a variety of forms, including monopiles, tripod or tetrapod structure resting on piles or buckets, hybrid monobuckets [1], high rise pile cap foundations supported by pile group [2]. The diversity in the foundation forms is attributed to the geotechnical conditions which often feature normally consolidated soft clays of varying thickness near the surface, followed by more competent clays, silty clays, loose to dense sands or silty sands. In some areas, shallow bedrock combined with weak seabed sediments may lead to expensive foundation construction in rock. From a geotechnical point of view, these site conditions are complex and challenging and represent a clear contrast to the conditions in the North Sea which typically consist of heavily over-consolidated clays and very dense sands. In Europe, most of the offshore wind turbines are founded on monopiles due to their proven robustness, China are considered. Lastly, the paper outlines a method to incorporate the cyclic effects into the currently proposed model.

The proposed multi-spring beam-column model for monopile analysis in soft clay
The load (moment)-displacement (rotation) response for soil-structure interaction is an integration of the stress-strain response of the soil over the volume that is mobilised by the structure. In the simplest form, the elastic stiffness coefficients are used to scale the soil stiffness to the foundation-soil interaction stiffness in ideally linearly elastic soil. For the case of nonlinear soil behaviour, it is observed that for many soil-structure interaction problems, the normalised load-displacement response is self-similar to the stress-strain response of the soil. Hence it is possible to scale the soil-structure interaction response from the stress-strain response with certain scaling coefficients. Examples can be found for shallow foundations [11,12] and a laterally loaded pile slice engaged in a plane strain flow-around mechanism [13][14][15][16][17][18].
The above-mentioned idea is exploited in the model proposed herein. This section outlines the proposed multi-spring beam-column model for monopile analyses. Formulations for lateral load-displacement (p-y) spring, base shear-displacement (s-y) spring and moment-rotation (m-θ) spring contributed by side shear are introduced in sequence. All the spring components are fundamentally linked to the soil stress-strain response. The advantage of this framework is that it allows easy implementation of site-specific soils springs that are based on soil response measured in routine laboratory tests.

Lateral resistance along the pile length (p-y component)
For a laterally loaded pile slice engaged in a plane strain flow-around mechanism, Zhang and Andersen [17] proposed a framework to derive site-specific p-y curves by scaling the soil stress-strain response measured in Direct Simple Shear (DSS) tests, as illustrated in Fig. 1. A point on the stress-strain curve with mobilisation in terms of shear stress (τ/s u ) corresponds to a point on the p-y curve with an equivalent level of mobilisation in terms of lateral bearing pressure (p/p u ). The corresponding normalised lateral displacement (y/D) can be scaled from the shear strain (γ) on the DSS stress-strain curve.
Eqs. (1)-(4) express this scaling relation mathematically: where ξ e and ξ p are elastic and plastic scaling factors respectively. According to Zhang and Andersen [17]; ξ e is taken to be a constant value of 2.6 for simplicity, while ξ p is expressed as a function of pile-soil interface roughness factor (α): Note that the p-y framework proposed by Zhang and Andersen [17] is strictly only applicable for a plane-strain flow-around mechanism. Fig. 2 illustrates the possible failure mechanisms for monopiles in soft clay. For a relatively short monopile, the soil will react in a wedge mechanism in the upper part and a rotational mechanism in the lower part (Fig. 2a). Whereas for a longer monopile (Fig. 2b), there may exist an intermediate zone with flow-around mechanism. This has been illustrated by PIV centrifuge experiments [19]. Note that in Fig. 2, it has been assumed that suction is present on the rear side of the monopile therefore an active soil wedge is mobilised. This is considered a reasonable assumption in soft clay conditions concerned in this paper. Despite the difference in failure mechanisms, the current model recommends use of the proposed p-y framework along the entire pile length. The impact of this assumption will be examined later when the model is used to back-analyse the pile response calculated by finite element analyses. where N p is the lateral bearing capacity factor. The recommendation by Zhang et al. [20], derived from upper bound analyses of Yu et al. [21], for the condition with suction is adopted: where N p0 is the lateral bearing capacity factor mobilised in a passive soil wedge in ideally weightless soil. The factor of 2 in Eq. (7) reflects that when suction is assumed, an active wedge is also mobilised therefore doubling the resistance (assuming isotropy in soil strength). N pd is the limiting bearing capacity factor mobilised in a plane strain flow around mechanism, which is found to be a function of interface roughness α according to the limit analysis solution by Randolph and Houslby [22]. As suggested by Zhang et al. [20]; N p0 and N pd can be calculated by Eq. (8) and Eq. (9) respectively.

Fig. 3.
Simplification of the soil mechanism at the pile tip [7].

Fig. 2.
Soil reacting mechanism for monopiles in soft clay. Note that Eqs. (1)-(5) provide the normalised shape of the p-y curve. The ultimate strength of the p-y spring p u is calculated as: where z ¼ depth below seafloor in consideration The parameter λ is a normalised parameter which measures the strength heterogeneity with depth. s um is defined as the soil strength at the mudline or extrapolated intercept at mudline in the case of layered soil conditions(kPa), k is the strength variation gradient with depth (kPa/m).

Lateral shear resistance at the pile tip (s-y component)
Zhang and Andersen [7] examined the contribution of pile tip resistance to the overall monopile response. They found that for most practical scenarios, it is acceptable to only consider the shear component at the pile tip, i.e. simplifying the rotation scoop mechanism at the pile tip to a direct simple shear mechanism, as illustrated in Fig. 3. This is primarily because the point about which the rotation scoop occurs is typically sufficiently high that the horizontal shear component contributes the majority of the total rotational resistance from the pile tip (Fig. 3). By making this simplification, Zhang and Andersen [7] established the scaling relation from shear force (s) vs relative lateral displacement (y) to the stress-strain response of the soil at the pile tip, as described by Eq. (10). In this case, the stress-strain behaviour from the direct simple shear mode is most relevant.
where s is the currently mobilised base shear resistance and s ult is the ultimate base shear resistance. The normalised lateral displacement at pile tip y/D can be directly scaled from the stress-strain curve, as described by Eq. (2). The scaling coefficients are constant with ξ e ¼ 0.3 and ξ p ¼ 0.12. The elastic and plastic shear strain components are calculated according to Eq. (3) and Eq. (4) respectively.
The ultimate base shear resistance s ult is simply the soil undrained shear strength multiplying the cross-section area (A) of the pile, i. e. s ult ¼ As u ¼ πD 2 s u /4.

Fig. 4.
Derivation of ultimate moment resistance due to side friction.

Rotational resistance along the pile shaft (m-θ component)
The moment-rotation (m-θ) spring due to the shaft friction has not been investigated in detail before. In this study, a theoretical derivation is performed. First of all, the analytical expression for the ultimate moment resistance is explored. The link between the moment-rotation (m-θ) spring to the shaft friction-axial displacement (t-z) spring is then established. Lastly, through the relation between t-z spring and the stress-strain response of the soil, the relation between m-θ spring with the soil stress-strain response is established. The following sections will detail the derivations.

Ultimate moment resistance due to side friction
Left plot of Fig. 4 presents a schematic illustration of mobilised shaft friction due to rotation of the pile cross-section. At small rotation, the largest shaft friction is mobilised furthest away from the pile centre where the relative soil-pile displacement is greatest. At large enough rotation, the shaft friction can be assumed to be fully mobilised all-around the pile shaft, which will provide the ultimate moment resistance. Right plot of Fig. 4 presents a plan view illustration for mathematical derivation. It follows that: Eq. (11) derives the theoretical ultimate moment resistance for a pile segment of unit height when the pile-soil interface strength (τ max ) is fully mobilised everywhere along the pile perimeter.

Link between t-z spring and m-θ spring
Assuming that the soil-pile interface shearing interaction under rotation is the same as that under uniform axial displacement, a link can therefore be established between the t-z spring and the m-θ spring. For a linear t-z spring that has the following relation: where z is relative axial soil-pile displacement, z max is the relative displacement at which the ultimate soil-pile interface shear stress τ max is mobilised.
Referring to Fig. 4, for a pile segment rotating about its centreline with a small rotation θ, the relative vertical soil-pile displacement (z) is: The corresponding moment resistance is Scaling relationship between t-z spring and m-θ spring.
Replacing M max as defined by Eq. (11), Eq. (14) can be re-written as: For M/M max ¼ τ/τ max , it requires that: which means that: Eq. (17) suggests that a simple scaling relationship between the t-z spring and the m-θ spring exists, which is illustrated by Fig. 5.
For the same spring mobilisation level, i.e. when M/M max ¼ τ/τ max , the rotation (θ) can be scaled from the normalised displacement (z/ D) of the t-z spring by a factor of 8/π. It should be noted that the above scaling relationship is derived based on a linear t-z spring. For a nonlinear t-z relationship, it is more involved to derive the integral of the moment resistance due to side shear analytically. However, this can be easily solved numerically in a simple Microsoft Excel spreadsheet. Fig. 6 illustrates an example where the m-θ spring is scaled from a nonlinear t-z spring which has a hyperbolic shape. It can be seen that the scaling factor derived from a linear t-z spring works well for a nonlinear t-z spring. Parametric analyses have been performed for other forms of t-z formulations, and the scaling law is found to be general.

Theoretical t-z curves
Based on the work of Randolph and Wroth [23]; Kraft et al. [24] presented a model for theoretical t-z response of axially loaded piles by assuming: 1) The displacement pattern of the soil can be modelled as concentric cylinders in shear. 2) Radial soil displacements due to pile loads are assumed negligible when compared to vertical soil deformations. Therefore, simple shear conditions prevails in the soil.
3) Shear stress decreases with distance such that τ r ¼ τ 0 r 0 ; in which τ ¼ the shear stress at distance r; τ 0 ¼ the shear stress at pile soil interface; and r 0 ¼ the pile radius. 4) The shear stresses are negligible beyond a radial distance r m (or zone of influence) and soil does not deform beyond that point. Fig. 7 illustrates the concentric deformation mechanism assumed in the theoretical derivation. The vertical displacement z at the pile-soil interface is therefore an integration of shear strain γ over the radial distance from r 0 to r m : where G r is the shear modulus of the soil at a radial distance (r) from the centreline of the pile. For a linear elastic material, i.e. constant G, Eq. (18) gives the elastic component of vertical displacement z e : Manipulating Eq. (19) leads to: where γ e 0 is the shear strain at the pile-soil interface under a shear stress of τ 0 .
Eq. (20) simply suggests that for a linear elastic material, the axial displacement for the t-z curve can be scaled by an elastic scaling factor (ξ e ) of 0.5ln (r m /r 0 ). Apparently ξ e increases (i.e. elastic t-z stiffness decreases) with increase of r m /r 0 . For practical purpose, r m / r 0 ¼ 10 is assumed in this study, which leads to ξ e ¼ 1.15.
Kraft et al. [24] further considered the case where the shear modulus reduces with shear strain according a hyperbolic relationship. However, the considered relationship is found not flexible enough and therefore an alternative relation adopted by the NGI-ADP model [25] is considered in this study. The expression adopted by the NGI-ADP model is demonstrated to be capable of describing the stress-strain response of different clays under undrained conditions using two parameters: G max /s u and plastic shear strain at failure γ p f [17]. The G max can be determined through resonant column test or shear wave velocity measured through bender elements integrated in a triaxial test or a DSS test. The γ p f , despite its physical meaning of plastic shear strain at failure, is a curve fitting parameter that provides the best fit to the stress-strain response measured in the soil element shearing test.
The NGI-ADP model expresses the shear strain as summation of an elastic component (γ e ) and a plastic component (γ p ), as described by Eqs. (21) and (22): The axial displacement for the t-z response is therefore also a summation of an elastic component and a plastic component. The solution for the elastic component is already given by Eq. (20). The plastic component is studied below.
At a distance r from the pile centreline with a shear stress of τ 0 (r 0 /r), the corresponding plastic shear strain γ p can be solved from the following equation given by: It can be demonstrated that the solution for Eq. (23) is: Fig. 7. Deformation pattern and shear stress, strain distribution assumed for developing theoretical t-z response of the axially loaded pile.
where b is τ 0 /s u and x is r/r 0 .
The solution provided above is approximate as it removes the (r m /r 0 ) term. Effectively, this corresponds to the solution for an influence zone that extends infinitely (r m /r 0 ¼ ∞).
Eq. (25) suggest a plastic scaling relation factor: ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi where γ p 0 is the plastic shear strain at the pile-soil interface due to the shear stress τ 0 . It can be seen that the plastic scaling factor (ξ p ) is not a constant, but changes with the mobilisation level (b, which is τ 0 /s u ). This is illustrated in Fig. 8. The value of ξ p falls in the range from 1/6 to 1/2. For practical purpose, a constant value of 0.45 is suggested, which leads to slight over-estimation of the displacement at high mobilisation levels (i.e. softer t-z response).
The total normalised vertical displacement z/D, which consists of an elastic and a plastic component, is: where γ e and γ p are calculated according to Eq. (3) and Eq. (4) respectively. Based on the above discussions, the scaling coefficients are assigned the following constant values: ξ e ¼ 1.15. ξ p ¼ 0.45 Note that the above solution is established for an idealised situation where soil remains intact around the pile and a fully rough interface (i.e. the maximum interface shear stress is equal to the shear strength of adjacent soil). However, in reality, there will be a zone of 1-1.5 times the wall thickness of the pile that is highly remoulded immediately after pile installation [26][27][28]. After set-up, the strength in the remoulded zone partially recovers and so does the stiffness. There are two aspects of the installation effects: i) the interface roughness α which limits the maximum attainable interface shear stress compared to the strength of intact soil; ii) the impact of changed soil stiffness of in the remoulded zone on the t-z stiffness. Incorporation of the interface roughness into the t-z formulation is straightforward. With regard to the latter effect, it is however considered relatively insignificant in that the axial displacement (z) is an integration of shear strain over a large radial distance. The thin remoulded zone, which may have a different (typically degraded) stiffness compared to the soil further away, only contributes to a small proportion to the total axial displacement. Base on the assessment, it is most important to account for the interface roughness on the maximum strength of the t-z spring.

Link from stress-strain curve to moment-rotation spring
In Sections 3.3.2 and 3.3.3, the links between the m-θ spring, the t-z spring and eventually the stress-strain behaviour of the soil are explored. The scaling relation from the soil stress-strain behaviour to the moment-rotation spring is established through t-z spring and illustrated in Fig. 9. For an interface roughness factor α less than 1 (i.e. not fully rough), a cut-off is applied to the scaling relationship.
An important point to note is that the moment resistance spring due to shaft friction proposed herein is formulated only in terms of pile cross-section rotation. In reality, since the clay is loaded undrained (i.e. incompressible), relative soil-pile vertical displacement will also be induced due to lateral movement of the pile alone. This is however difficult to quantify in a mathematical form and is neglected herein. The current model will therefore tend to under-estimate the rotational moment due to shaft friction for a given crosssection rotation. Furthermore, the link between the moment-rotation spring and the soil stress-strain response is established assuming that the interface t-z response under rotation is similar to that under uniform axial displacement. This is an important assumption and the performance of the model is to be verified against the finite element simulations.

Summary of the multi-spring beam-column model
The above sections have detailed how to derive the p-y, s-y and m-θ springs from stress-strain response of the soil. Table 1 provides a summary of the multi-spring beam-column model.

Finite element analyses
In order to check the performance of the proposed model, a comprehensive set of finite element analyse (FEA) was carried out. The results of the FEA are back-analysed with the proposed model to check its performance. The FEA investigate several soil strength profiles and monopile geometries. The details of the FEA model, geometry, soil strength profile and soil constitutive model will be described below. Fig. 10 illustrates an example FEA model developed in ABAQUS [29] for a monopile with L/D ¼ 5. Due to symmetry, only half of the pile cross-section is simulated. The monopile is modelled as a solid pile instead of a pipe pile in reality. Equivalent elastic modulus is specified to the pile material so that the same bending stiffness (EI) as the actual pile is modelled. In the current analyses, monopiles with a diameter (D) of 6 m and a uniform wall thickness of 60 mm are considered. The penetrated length (L) varies from 18 m to 60 m, covering an L/D range from 3 to 10. The pile stick-up above the soil surface is kept 30 m for all L/D ratios, where a lateral force is applied. First order, reduced integration brick elements (C3D8R) were used to discretise both the soil and the pile. For numerical accuracy, a cylindrical zone with a radius of 2D from pile centre is discretised with refined mesh. Thin layers of elements are used to accurately calculate the soil resistance at the pile tip. The finite element mesh is consistent with the model adopted in Zhang and Andersen [7].

Finite element model and geometry
On the pile-soil interfaces (both shaft and tip), fully rough interface is modelled by defining surface to surface contact pairs. No separation is allowed on the interfaces, therefore simulating a condition without tension gapping. This is considered a reasonable assumption for the soil profiles considered in this study, as will be detailed below.

Soil model
In all FEA analyses, the Mohr-Coulomb model was used. By specifying a 0 � friction angle, it simplifies to a Tresca model. A Poisson's ratio υ of 0.48 was used, simulating undrained conditions. An elastic modulus E was specified which remains constant and was used to calculate the elastic deformation over the entire stress range. The equivalent elastic shear modulus G can therefore be calculated by: For undrained conditions, with υ ¼ 0.5, thus G ¼ E/3. In ABAQUS, it allows for specifying a table of plastic shear strain γ p versus mobilised shear stress τ to consider the nonlinear stressstrain material behaviour. In this study, the stress-strain response (plastic component) was calculated using the hardening relation adopted by the NGI-ADP model [25], as expressed by Eq. (21) and then used as input in ABAQUS. The complete stress-strain curve is captured using two parameters: G max /s u and γ f p and the total strain under a mobilised shear stress τ is a summation of an elastic component γ e and a plastic component γ p , as expressed by Eq. (22). Fig. 9. Scaling of moment-rotation spring from soil stress-strain response.
In the FEA parametric analyses, the following combination of stress-strain parameters were considered: � G max /s u of 333, 1000 � γ f p of 0.02, 0.10.

Soil strength profiles
The FEA modelling examines four different undrained strength profiles, as illustrated by Fig. 11. Soil profile 1 corresponds to an idealised normally consolidated soft clay with strength increasing linearly with depth. Soil profile 2 is lightly over-consolidated with a small strength intercept at the seafloor. Soil profile 3 and profile 4 are layered profiles, where the seabed consists of a normally consolidated soft clay layer followed by a more competent layer with constant strength with depth.
Due to the presence of soft clay near the surface of the seabed, the collapse of which will block any seepage path that is required for a tension gap to develop, suction (i.e. no tension gap) is assumed in all FEA performed in this study. It should be noted that the soil model adopted in this study is isotropic. The strength and stiffness anisotropy is an important aspect of natural soil behaviour. However, it is currently not considered to limit the complexity. Zhang et al. [20] evaluated the effect of strength anisotropy on the ultimate soil resistance to laterally loaded piles. They found that the DSS strength governs the ultimate resistance for the flow-around mechanism. For the wedge mechanism, if both active and passive wedges are mobilised (as currently modelled), the mechanism involves a mix of triaxial compression, triaxial extension and DSS shear modes and the DSS strength can be used as an average strength to computer the ultimate soil resistance. In the analyses performed herein, the isotropic strength modelled can be taken as the strength measured in DSS tests. The impact of strength anisotropy is considered to be minor, while the impact of stiffness anisotropy (different normalised stress-strain responses in different shear modes) remains an aspect to be looked further into in future studies.

Table 1
Summary of the multi-spring beam-column model. It is should be noted that the self-weight of the soil does not have an impact on its mechanical response when the undrained behaviour is modelled with a total stress based soil model. Furthermore, the self-weight does not contribute to the soil resistance when double wedges are mobilised as the self-weight induced soil pressure on both sides of the pile cancels out. For these reasons, the FEA assumed weightless soil for simplicity.

Notes on the "wished-in-place" assumption
It should be noted that in all analyses performed in this study, the pile is assumed to be "wished-in-place". During pile installation, the soil immediately adjacent to the pile wall in a zone that is 1-1.5 times of the wall thickness will be highly remoulded, as already mentioned before. Partial regain of the strength will occur after installation through a combination of immediate pore pressure change as the pile tip passes by, thixotropy, and pore pressure dissipation (which is often referred as "set-up" altogether). In terms of lateral p-y curves, the stiffness of the curves is expected to be mainly dominated by the bulk of soil outside the influence zone of pile installation; while for ultimate capacity of the p-y curves, it is considered sufficient to use the interface roughness factor to capture the impact of pile installation. This argument is also applicable to the moment-rotation spring, as already discussed in section 3.3.3. In terms of the base shear s-y curve, it is mainly contributed by the internal soil plug, which is believed to be largely intact for the large diameter monopiles installed by impact driving.
In this study, the pile installation process is not explicitly modelled. Instead, intact soil condition is assumed around the pile as the purpose of the numerical analyses is to verify the proposed framework. In application of the framework for actual design, the installation effect should be accounted for by evaluating an interface roughness and incorporating the effect in the strength and stiffness of the soil springs.     Fig. 12 presents the finite element analysis (FEA) results and back-analyses (B-A) using the proposed beam-column model in the ideally normally consolidated clay (Profile 1). Results for four L/D ratios are presented, which are for a stress-strain response described by the NGI-ADP model with a G max /s u ¼ 333, and a γ f p ¼ 0.1. To illustrates the effect of different soil springs, back-analyses are performed with:

FEA results in single layer clay and back-analyses
� Only lateral p-y springs � Lateral p-y springs and base shear s-y spring � Lateral p-y springs, base shear s-y spring and side shear m-θ springs.
The following observations can be made: 1) For a pile with L/D ¼ 10, the pile exhibits large deformation due to bending and the pile toe has little displacement/rotation. The proposed beam-column model captures the pile response very well. The effect of the base-shear spring and the moment-rotation springs is minimal and can be neglected for practical purposes. 2) For a pile with L/D ¼ 7, large deformation due to bending is observed at relatively small loads. However, the pile toe kicks when the load is high enough. The proposed beam-column model with all spring components captures the pile response well. The effect of the base shear and moment-rotation springs is noticeable. However, considering the p-y curves alone still provides a reasonably good prediction (albeit slightly on the softer side). 3) For a pile with L/D ¼ 5, most of the pile deflection is due to the rigid body displacement and there is only a small contribution from the bending of the pile. The rotation is highest at the loading point and decrease only moderately along the pile length. The proposed beam-column model provides a fair prediction of the pile response, although over-predicts the deflection at the pile head by around 30%. The base shear and moment-rotation springs are very important and neglecting of these springs will lead to large over-estimation of the pile defection and rotations. 4) For a pile with L/D ¼ 3, the pile cross-section rotation is almost constant along the pile length. Despite the extremely small aspect ratio, the proposed beam-column model with all spring components provides a fair prediction of the response that lies on the softer side. Neglecting the base shear and moment-rotation springs leads to almost twice amount of lateral deflection and rotation.
Very similar results are obtained in soil profile 2, which has a small strength intercept at the mudline. Therefore, the full set of comparison between the FEA results and the back-analyses is not presented herein. To put the results into perspective, the FEA results for L/D ¼ 5 are compared to predictions made by the API p-y springs and by the currently proposed model, as illustrated in Fig. 13.
While the currently proposed model provides a good match, the API p-y springs (constructed using ε 50 ¼ 0.01 and J ¼ 0.5) excessively over-predicts the pile deflection and rotation. Predictions are also made by API model after including the currently proposed base shear (s-y) and moment-rotation (m-θ) side shear springs. Fig. 13 illustrates that the lateral deflection and rotation of the pile are still grossly over-predicted. This comparison highlights that the poor performance of the API model for monopile design is not only due to ignoring the base shear and side shear resistances, but also due to the p-y formulation itself. Amongst others, Jeanjean et al. [18] advocated a shift from the API p-y spring formulation in clay even for conventional slender pile design as it poorly captures the lateral soil resistance and stiffness along the pile. Zhu et al. [30] presented field test results of two offshore driven piles in soft clay, which echoed similar conclusions.

FEA results in layered clay and back-analyses
The performance of the proposed model is also checked in the layered soil profiles. Fig. 14 presents the FEA results together with back-analyses in soil profile 4. Despite the contrast in strength between the soil layers, the proposed model performs very well and capture the pile response even for small L/D ratios. The comparison in soil profile 3 is similar and is therefore not presented herein.

Effect of cyclic loading
A framework based on a conventional beam-column approach for analysing the monopile response in soft clays and layered profiles is described and validated above. The framework fundamentally links the soil-structure interaction to the stress-strain responses of the soil. Monopile foundations supporting offshore wind turbines are subjected to cyclic loading, which alters the soil's strength and stiffness from those measured under monotonic loading. The cyclic effects must be properly accounted for in design. In this section, an approach to incorporate the cyclic effects is briefly outlined.
To account for the effect of cyclic loading in offshore foundation design, the concept of equivalent number of cycles (N eq ) which relates to the largest load cycle, is often used to describe the equivalent effect of a load history comprising cycles of different amplitudes. Cyclic shear strain or accumulated pore pressure is typically used as a state parameter for calculating N eq of a cyclic load history. Andersen [31] provides a comprehensive description of the cyclic accumulation procedure. Once N eq has been determined for the given design load history, the strength and stress-strain response of the soil can then be extracted from cyclic strain contour diagrams, which are established from systematic testing of soil elements under cyclic shearing. An example is shown in Fig. 15a, which illustrates a cross-section of a strain contour diagram for 10 cycles established from cyclic DSS tests. A point on the contour diagram represents the average and the cyclic shear strains that are developed after 10 cycles of shearing under the corresponding combination of average and cyclic shear stresses. By considering the relevant load path (i.e. the ratio between the cyclic amplitude, τ cy , to the average component, τ a , which is typically assumed to be the corresponding ratio for the pile head load), the points at which the loading path intersects the contour lines form the stress-strain response of the soil under cyclic loading, as illustrated in Fig. 15b. It should be clarified that the stress-strain response illustrated does not describe the soil behaviour within an individual cycle, but rather describes the total strain (sum of the average component and the cyclic amplitude) after 10 cycles under the combination of average and cyclic shear stresses. Using the total stress-strain curve, the soil reaction curves (p-y etc.) can then be constructed using the framework described above and the same scaling coefficients as listed in Table 1. The pile response under the peak loading (the sum of the average component and the cyclic amplitude) can then be evaluated.  However, if the soil-pile interaction stiffness within an individual load cycle is needed, for example, to be used in the dynamic analysis of the wind turbine system, the cyclic component of the stress-strain curve (τ cy versus γ cy ), which describes the peak-trough secant stiffness within an individual cycle, should be used to construct the soil springs (using the model and scaling coefficients  Table 1). Accordingly, the cyclic amplitude of the pile head load should be applied to calculate the cyclic foundation stiffness.
It should be noted that the method outlined above assumes the same N eq along the depth of a soil layer, i.e. the cyclic effect is independent of depth. This is clearly a rough assumption as the soil along the upper part of the pile is more heavily loaded and therefore greater cyclic effect is expected. This will lead to a redistribution of stresses along the pile. Despite that, the approach outlined above is considered reasonable as a preliminary method to take cyclic loading into account in design. A more refined approach which captures the different levels of cyclic effect along the pile can be found in Zhang et al. [32,33].

Limitations
In the current study, the soil has been simplified as an isotropic material. However, natural soil shows anisotropic strength and stress-strain response. This aspect is to be further investigated in future studies. As an interim recommendation, the shear strength and stress-strain response measured in the DSS shear model is suggested to be used when applying the proposed model. The DSS strength is found to be typically a good average of triaxial compression, DSS and triaxial extension strengths [34,35]. A further extension of the study is to investigate the model performance in soil profiles with more complex layering, for example inclusion of sand layers.

Concluding remarks
This paper has presented a multi-spring beam-column model for monopile analyses in soft clays. The model takes into account the lateral p-y springs, the base shear s-y spring and moment-rotation (m-θ) spring contributed by shaft friction. All the spring components are fundamentally linked to the stress-strain response of the soil that is measured in conventional laboratory testing. The model provides an analytical framework to derive different soil-pile interaction springs which are conventionally derived empirically from field or laboratory model testing. However, pile load testing is time consuming and expensive, and the applicability of the empirically derived soil springs is often limited to the soil condition tested. The proposed model allows for easy implementation of site specific curves in design, without the need of time consuming and often complicated finite element analyses. Examined against a comprehensive suite of 3D finite element parametric analyses with consideration of several typical soil conditions encountered in Chinese offshore wind farms, the proposed framework demonstrates good predictions to the pile responses, with the important contribution from the pile tip and side shear friction for monopiles (especially for lower L/D cases) well captured. Compared with the current design practice (i.e. the API p-y curves), the model represents a significant improvement, which can have important implications for the monopile design.
The monopiles supporting offshore wind turbines are subjected to constant cyclic loading. An approach to incorporate the cyclic loading effect into the proposed beam-column framework is also briefly outlined in this paper.