This paper describes a onedimensional (1D) computational model for the analysis and design of laterally loaded monopile foundations for offshore wind turbine applications. The model represents the monopile as an embedded beam and specially formulated functions, referred to as soil reaction curves, are employed to represent the various components of soil reaction that are assumed to act on the pile. This design model was an outcome of a recently completed joint industry research project – known as PISA – on the development of new procedures for the design of monopile foundations for offshore wind applications. The overall framework of the model, and an application to a stiff glacial clay till soil, is described in a companion paper by Byrne and coworkers; the current paper describes an alternative formulation that has been developed for soil reaction curves that are applicable to monopiles installed at offshore homogeneous sand sites, for drained loading. The 1D model is calibrated using data from a set of threedimensional finiteelement analyses, conducted over a calibration space comprising pile geometries, loading configurations and soil relative densities that span typical design values. The performance of the model is demonstrated by the analysis of example design cases. The current form of the model is applicable to homogeneous soil and monotonic loading, although extensions to soil layering and cyclic loading are possible.
Monopiles are typically the preferred foundation option for offshore wind turbine support structures in shallow coastal waters. Current design procedures for monopile foundations routinely employ a simplified analysis procedure, known as the ‘p–y’ method, in which the foundation is modelled as an embedded beam, with the lateral load–displacement interaction between the soil and pile represented by nonlinear functions known as ‘p–y curves’.
The p–y method was originally devised for the design of the long, relatively flexible, piles that are typically employed in offshore oil and gas structures. The method was initially based on data from field tests reported some decades ago (e.g. Matlock, 1970; Cox et al., 1974); early p–y curve specifications were proposed by Matlock (1970) (for clays) and Reese et al. (1974) (for sands). Although the method has evolved in the intervening years (e.g. Doherty & Gavin, 2011), current standard forms of the p–y method as specified in design guidance documents (e.g. API, 2010; DNV GL, 2016) remain broadly unchanged from this early work. Certain questions exist, however, on the extent to which standard forms of the method are applicable to offshore wind turbine monopiles, which typically employ relatively large diameters, D, and low values of L/D (where L is embedded length) and are therefore relatively stiff. Evidence highlighting the shortcomings of the conventional p–y method for monopile design applications has been observed in laboratory tests (e.g. Choo & Kim, 2015; Klinkvort et al., 2016) and at field scale (Kallehave et al., 2015; Li et al., 2017; Hu & Yang, 2018).
This paper describes a new analysis procedure, referred to as the ‘PISA design model’, for monotonic lateral and moment loading of monopiles. This design model is an outcome of a research project – known as PISA – that included field testing (Burd et al., 2019; Byrne et al., 2019b; McAdam et al., 2019; Zdravković et al., 2019a) at two onshore sites (stiff clay at Cowden, dense sand at Dunkirk) and threedimensional (3D) finiteelement modelling (Taborda et al., 2019; Zdravković et al., 2019b). The PISA design model retains the underlying simplicity of the p–y method (in which the pile is modelled as an embedded beam), but additional soil reaction components are incorporated to improve the model's performance. The model is calibrated with a set of 3D finiteelement calibration analyses; it therefore benefits from the realism that is potentially achievable with 3D finiteelement modelling, while also being rapid to compute. The PISA design model supports a wide range of practical design calculations; it is applicable to: (a) the determination of small displacement foundation stiffness (relevant to the development of dynamic models for the overall structure); (b) the analysis of serviceability limit states (i.e. relating to the displacements that occur under normal working conditions); and (c) analysis of ultimate limit states (to check for overall stability).
This paper describes the development and implementation of the PISA design model for monopiles installed in homogeneous sand for drained monotonic loading. In a companion paper, Byrne et al. (2020), the overall framework of the PISA design model is described, together with a calibration process for piles embedded in glacial clay till. The approaches employed for the clay and sand PISA design model formulations differ only in the manner in which the soil reaction components are incorporated and the way in which the model is calibrated. The particular aspects of the model that relate to the sand implementation are referred to in the current paper as the ‘sand modelling framework’.
To develop the PISA design model in a form that is applicable to sands within a practical range of densities, a set of four hypothetical representative offshore homogeneous sand sites is established, each with a specific relative density (D_{R} of 45%, 60%, 75% and 90%). The geotechnical conditions at these sites, and the modelling employed in the 3D calibration analyses, are based on the prior geotechnical characterisation of the Dunkirk site (Zdravković et al., 2019a) and the finiteelement analyses that were shown in the paper by Taborda et al. (2019) to provide a close representation of the PISA test piles at Dunkirk. The PISA design model calibration process therefore has a link, albeit an indirect one, with observations on the performance of the PISA test piles. Independent PISA design model calibrations are described for each representative site, and an optimisation process is employed to define a general model – referred to as the ‘general Dunkirk sand model’ (GDSM) – that is applicable to soils with an arbitrary value of relative density in the range 45% ≤ D_{R} ≤ 90%. The predictive capabilities of the GDSM are demonstrated by conducting analyses for monopile configurations within the calibration space, but that differ from the calibration cases.
The PISA design model provides a onedimensional (1D) representation of a monopile foundation subject to the application of a lateral load, H, applied at a distance h above seabed level (referred to in this paper as ‘ground level’) as illustrated in Fig. 1. The monopile is represented as an embedded beam with moment M_{G} and lateral force H_{G} applied to the pile at ground level, Fig. 2. Four components of soil reaction are assumed to act on the monopile. Consistent with the standard p–y method, a distributed lateral load, p (units of force/length) acts on the pile. Additionally, a distributed moment, m (units of force × length/length) is applied; this distributed moment arises as a consequence of the vertical tractions that are induced on the pile perimeter when relative vertical displacements occur at the soil–pile interface – for example, due to local rotation of the pile crosssection. A lateral force H_{B} and a moment M_{B} acting on the base of the pile are also included. The monopile is represented by Timoshenko beam theory; this allows the shear strains in the pile to be incorporated in the analysis in an approximate way. Since the influence of the shear strains on the overall pile deformation is likely to increase as L/D is reduced, the use of Timoshenko theory provides a means of maintaining the robustness of the approach as L reduces or D is increased (e.g. see Gupta & Basu (2018)). A fourcomponent model of this sort has previously been employed for the design of drilled shafts for onshore applications (e.g. Lam, 2013) and has been described in the context of the PISA research by Byrne et al. (2015), Byrne et al. (2017) and Burd et al. (2017). As discussed in the paper by Byrne et al. (2020), vertical loads are assumed to have an insignificant influence on the performance of the monopile; they are therefore excluded from the model.
The soil reactions are applied to the embedded beam using a generalised form of the Winkler assumption, in which the force and moment reactions are assumed to be related only to the local pile displacement and rotation. Functions employed in the model to relate the soil reactions and the local pile displacements (and rotations) are termed ‘soil reaction curves’. Although the Winkler approach neglects the coupling that inevitably occurs within the soil, it provides a convenient basis for design calculations, as demonstrated by the widespread use of the p–y method. A fundamental feature of the approach, however, is that soil reaction curves determined on the basis of the Winkler assumption are unlikely to be unique. Appropriate soil reaction curves may depend, for example, on the relative magnitude of the translational and rotational movements of the pile. It is considered necessary, therefore, to calibrate the soil reaction curves using pile deformation modes that are representative of those that are likely to be experienced by actual wind turbine monopile foundations. The PISA design model is therefore calibrated within a design space that is carefully selected to represent realistic loading conditions.
The PISA design model reduces to the standard p–y approach when m, H_{B} and M_{B} are set to zero (and appropriate choices are made on the relationship between the distributed lateral load, p, and the local lateral pile displacement, v). Experience has shown, however, that the m, H_{B} and M_{B} components become increasingly significant as L/D is reduced (Byrne et al., 2015, 2020). The distributed moment component, for example, depends on pile diameter, increasing as the pile diameter is increased. Similarly, the force and moment reactions H_{B} and M_{B} at the base of the pile become more significant as the pile diameter is increased. The fourcomponent model in Fig. 2(b) therefore provides a rational way of addressing a shortcoming of the p–y method, often referred to as the ‘diameter effect’, in which the standard p–y curves (e.g. API, 2010; DNV GL, 2016) are typically found to become increasingly unreliable as the pile diameter is increased, or the length is reduced (e.g. Alderlieste et al., 2011; Doherty & Gavin, 2011).
The soil reaction curves employed in the PISA design model are based on the use of dimensionless forms of the relevant soil reaction and displacement/rotation variables. This provides a convenient means of developing standard forms that, for numerical implementation, can be scaled to represent the soil reactions acting on the pile at an arbitrary depth. These dimensionless forms are specified in Table 1, where σ′_{vi} is the local value of initial vertical effective stress in the soil, G_{0} is the local value of soil smallstrain shear modulus, and v and ψ are the local pile lateral displacement and crosssection rotation, respectively.

Normalised variable  Dimensionless form 

Distributed lateral load,  
Lateral displacement,  
Distributed moment,  
Pile crosssection rotation,  
Base horizontal load,  
Base moment, 
The soil reactions are implemented in the model using an appropriately calibrated algebraic function. The function selected for this purpose is, to an extent, arbitrary, provided that it is capable of providing a realistic representation of the soil reactions for behaviour ranging from small displacements (needed, e.g. to predict the natural frequencies of a wind turbine structure) to the large displacement response (required for the calculation of the ultimate limit state).
This coupling has certain implications for the numerical implementation of the model (described in the Appendix).
The PISA design model employs the 1D representation finiteelement framework illustrated in Fig. 2(b). The pile is represented by a line mesh of twonoded Timoshenko beam elements, employing the formulation in Astley (1992). The calculations described in this paper were all conducted with a shear factor κ = 0·5. In the current form of the model, consistent with the shell element formulation employed in the 3D finiteelement calibration calculations, the structural properties (area and second moment of area) are specified for the beam elements using the thinwalled approximation. Soil finite elements, with the same displacement and rotation interpolation functions that are used for the beam elements, are connected to the beam elements along the embedded length of the pile. A virtual work statement of the problem, and an outline of the development of the finiteelement equations in Galerkin form, is provided in the paper by Byrne et al. (2020). Consistent with the approach adopted for the clay framework described in the paper by Byrne et al. (2020), four Gauss points per element are adopted for both the beam and soil elements to determine the stiffness matrices and internal force vectors. Further specific implementation details for the sand modelling framework are given in the Appendix.
In connection with the PISA research, a series of pile tests (McAdam et al., 2019) was conducted at an onshore site in Dunkirk in northern France; at this site the soil consists principally of a dense Flandrian sand with a surface layer (about 3 m thick) of dense, hydraulically placed sand with the same geological origin as the deeper Flandrian deposit (Chow, 1997). This site was carefully characterised during the field testing programme (Zdravković et al., 2019a); it was convenient, therefore, to adopt the soil conditions at this site as the basis of the representative offshore sites developed to calibrate the PISA model for sand.
A detailed finiteelement study (Taborda et al., 2019) was undertaken during the PISA project to support the Dunkirk pile tests. The constitutive model employed for these analyses, described in the paper by Taborda et al. (2014), is an evolution of the bounding surface model originally proposed by Manzari & Dafalias (1997). A detailed procedure to calibrate this constitutive model for the soil conditions at the Dunkirk test site, described in the paper by Taborda et al. (2019), was conducted. The constitutive model and associated parameters (see Table 3 later) that were developed to model the Dunkirk test piles are adopted for the representative offshore site calibration calculations employed in the current work.
The ground conditions at Dunkirk have certain features not present at typical offshore sites. These are: (a) a very dense, hydraulically placed surface layer; (b) the surface soil layers are partially saturated, as a consequence of the water table being observed to be 5·4 m below the ground surface; and (c) the superficial layers are possibly lightly cemented. Adjustments to the Dunkirk soil conditions were therefore required to develop plausible representative offshore ground models; these adjustments involved (a) excluding the surface layer of hydraulic fill from the model, and (b) employing a hydrostatic pore pressure distribution. Other aspects of the representative offshore ground models were taken directly from the data in the paper by Taborda et al. (2019) on the naturally occurring Flandrian sand at Dunkirk. The constitutive model and calibration data developed to support the Dunkirk field tests (Taborda et al., 2019) were employed directly to characterise the representative ground models for the current study; the only parameter requiring adjustment is the relative density.
The 3D finiteelement calibration calculations have been conducted for a calibration space consisting of monopile dimensions and load eccentricities, h, in the range 5 m ≤ D ≤ 10 m, 2 ≤ L/D ≤ 6, 5 ≤ h/D ≤ 15 for each representative site. These dimensions were selected to span a realistic design range for current and future monopiles, on the basis of advice received from the project partners (listed in the Acknowledgements). The configurations employed in this set of pile calibration analyses, selected to provide appropriate coverage of the selected calibration space, are listed in Table 4. The analyses were conducted using the finiteelement software ICFEP (Potts & Zdravković, 1999, 2001). In total, 38 calibration analyses were conducted (see Table 8 later).

Relative density, D_{R}: %  Initial void ratio, e_{0} 

45  0·741 
60  0·685 
75  0·629 
90  0·573 

Component  Parameters 

Critical state line  p′_{ref} = 101·3 kPa; e_{CS,ref} = 0·910; λ = 0·135; ξ = 0·179 
Strength  
Model surfaces  
Hardening modulus  h_{0} = 0·4; α = 1·0; γ = 0·0; β = 0·0; μ = 1·0 
Nonlinear elasticity – smallstrain stiffness  B = 875·0; ν = 0·17 
Nonlinear elasticity – shear stiffness degradation  a_{1} = 0·40; γ_{1} = 1·031 × 10^{−3}; κ = 2·0 
Fabric tensor  H_{0} = 0·0; ζ = 0·0 
These parameters are identical to those that were determined, as described in the paper by Taborda et al. (2019), to conduct 3D finiteelement analysis of the PISA test piles at Dunkirk.

Pile reference  D: m  h: m  h/D  L: m  L/D  t: mm  D/t 

C1  10  50  5  20  2  91  110 
C2  10  150  15  20  2  91  110 
C3  10  50  5  20  2  125  80 
C4  10  50  5  60  6  91  110 
C5  10  150  15  60  6  91  110 
C6  5  25  5  10  2  45  110 
C7  5  25  5  10  2  83  60 
C8  5  25  5  30  6  45  110 
C9  5  75  15  30  6  45  110 
C10  7·5  37·5  5  15  2  68  110 
C11  7·5  37·5  5  45  6  68  110 
Procedures to calibrate the model were initially developed on the basis of the D_{R} = 75% representative site. This initial calibration exercise was conducted for all of the calibration piles in Table 4. It is noted that piles C3 and C7 (incorporated in the pile calibration set to check whether the inferred soil reaction curves are influenced by pile wall thickness) are similar to C1 and C6, respectively, except for differences in wall thickness. Results from the D_{R} = 75% calibration indicated that the influence of wall thickness on the soil reaction curves is negligible; piles C3 and C7 were therefore excluded from the calibration sets employed for the other representative sites.
Results from the D_{R} = 75% calibration are employed later in the paper to illustrate the various stages in the calibration process.
The 3D finiteelement calculations employed a critical state constitutive model, based on the state parameter framework for sands (Taborda et al., 2014), to represent the soil at the representative sites. The state parameter framework employed in the model ensures that the influence of soil void ratio, and mean effective stress, on the mechanical behaviour of soil is accounted for in a consistent way – that is without the need to adjust the model parameters for soils with different relative densities. The constitutive parameters employed in the analyses (determined as described in the paper by Taborda et al. (2019)) are listed in Table 3. It is noted that these constitutive parameters were developed for monotonic loading and have not been calibrated for cyclic loads.
A typical mesh employed for the calibration analyses (for pile C4) is shown in Fig. 4. By exploiting symmetry in the geometry and in the applied load, only half of the problem is discretised. The soil domain is represented with 10 530 20noded hexahedral displacementbased solid elements. A refined mesh is employed for the soil in the region below the base of the pile to ensure that the computed base reactions are reliable. The embedded pile is discretised with 360 eightnoded shell elements (Schroeder et al., 2007) arranged in 30 rings of elements; distributed load and moment soil reactions curves could therefore be extracted from the calibration analyses at 30 discrete depths. The aboveground extension is modelled with 240 shell elements. The interface between the soil and the pile exterior is modelled with 360 16noded zerothickness interface elements (Day & Potts, 1994). Fully rough boundary conditions are prescribed to the base of the mesh and a zero normal displacement boundary condition was prescribed to the vertical cylindrical boundary (at a radial distance of 100 m from the pile central axis).
No attempt is made to model the stress and state changes that occur in the soil due to the pile installation process. Instead, the monopile is modelled as ‘wished in place’; it is incorporated in the finiteelement mesh at the start of the analysis in a fully plugged configuration. A similar wishedinplace procedure was employed in the 3D finiteelement models that were developed to analyse the Dunkirk PISA test piles (Taborda et al., 2019).
The interface between the exterior of the embedded monopile and the soil is represented by an elastoplastic Mohr–Coulomb model. The elastic part of the interface model is defined by a shear and a normal stiffness, both set to 1·0 × 10^{5} kN/m^{3}, and the plastic part by zero cohesion (c′ = 0) and an angle of shearing resistance (32^{°}) that is equal to the triaxial compression critical state friction angle. The monopile is modelled as an elastic material with properties representative of steel; Young's modulus, E = 200 GPa and Poisson's ratio, ν = 0·3. The pile wall thickness is specified as an additional model parameter.
Loading is applied to the top of the pile in a displacementcontrolled manner, by prescribing increments of uniform horizontal displacements in the ydirection around the halfperimeter of the pile. The resulting load is obtained as a reaction to these prescribed displacements; its magnitude is one half of the total lateral load, H, acting on the pile.
Numerical representations of the soil reaction curves (referred to as ‘numerical soil reaction curves’) for the distributed lateral load and moment components were determined from the 3D finiteelement analyses by extracting the nodal forces acting at the soil–pile interface, and the stresses in the interface elements between the pile exterior and the soil. The force and moment reactions at the pile base were determined by integrating the stresses in the layer of soil elements immediately below the pile base, and incorporating the reactions computed from the nodal forces at the base of the shell elements representing the pile. Local lateral displacements and crosssection rotations of the pile were determined from the computed displacements of the relevant shell element nodes by averaging over the crosssection (for displacement) and by leastsquares fitting on the vertical displacements (for rotation).
Checks were conducted to confirm that the computed nodal forces acting on the monopile boundary were in equilibrium with the externally applied lateral load (within an acceptable tolerance). If boundary checks were satisfactory, then the data were further processed to develop the soil reaction curves, as described below. Alternatively, if this boundary equilibrium check indicated the presence of unacceptable equilibrium errors, then the 3D analyses were repeated using a tighter calculation tolerance. For the calibration calculations (listed in Table 8 later) a maximum equilibrium error of 1·81% was achieved; this level of equilibrium error is considered to be well within the bounds of acceptability.
Check calculations were conducted using a form of the 1D model, referred to as ‘1D (numerical)’, that is based directly on the numerical soil reaction curves. In this approach, dimensionless forms of the numerical soil reaction curves are determined using the normalisations in Table 1. Normalised numerical soil reaction curves at the depth location of each Gauss point in the 1D model are computed by interpolation; the corresponding dimensional forms are then determined on the basis of the local values of σ′_{vi} and G, and the dimensionless form definitions in Table 1. The H–v_{G} performance (where v_{G} is groundlevel pile displacement) computed using the 1D (numerical) model for piles C1 (L/D = 2) and C4 (L/D = 6) for D_{R} = 75% is shown in Fig. 5. The 1D (numerical) model is seen to provide a close fit to the 3D finiteelement calibration data. A similarly close match is obtained for other calibration piles (data not presented here). These checks confirm that the procedures used to determine the numerical soil reaction curves are robust. They also indicate a likely upper bound on the accuracy of the PISA modelling approach.
Separate 1D (numerical) calculations have been conducted to investigate the significance of individual soil reaction components. Example results for piles C1 and C4 for D_{R} = 75%, for cases where soil reaction components are selectively excluded from the model, are also shown in Fig. 5. In case P, only the distributed lateral load terms are included; in case PM, only the distributed lateral load and distributed moment terms are included. It is clear from Fig. 5(b) that the lateral distributed load is the dominant soil–pile interaction mechanism for the relatively long pile, C4 (i.e. the case P data match closely the 1D (numerical) results). For the shorter pile (C1) Fig. 5(a), however, the case P data differ significantly from the 1D (numerical) model, indicating that, in this case, neglecting the three other soil reaction components causes a significant loss of fidelity. The case PM data provide an improved fit for pile C1, indicating the importance of the distributed moment in this case. These results confirm the pattern observed in the paper by Byrne et al. (2020) for a stiff glacial clay till, that for relatively long piles a p–y type method (distributed lateral load only) is capable of providing a robust model of the load–displacement behaviour, but that additional soil reaction components need to be included for piles with relatively low values of L/D.
Accuracy metric values have been computed for the 1D (numerical) model for ‘ultimate displacements’ η_{ult} determined for 0 < v_{G} < D/10 and ‘small displacements’, η_{sd}, determined for 0 < v_{G} < D/10 000 for all of the piles in the calibration set for all relative densities. For D_{R} = 75%, the accuracy metrics are in the range from 0·92 to 0·98 for η_{ult} and from 0·89 to 0·98 for η_{sd}. Values of the ratio metric evaluated at v_{G} = D/10 and v_{G} = D/10 000 for D_{R} = 75%, are in the range from 0·93 to 1·07 and from 0·88 to 1·08, respectively. These results indicate a close match between the 1D (numerical) model and the 3D calibration data. Similarly close agreement was obtained for the other relative density cases.
For a practical design tool, general forms of the soil reaction curves are required that are applicable to pile configurations not included in the calibration set. The current form of the PISA design model employs the fourparameter form in equation (1) to represent the soil reaction curves. Soil reaction curves based on this function are referred to as ‘parametric soil reaction curves’.
Values of the parameters required to fit the parametric soil reaction curves to the numerical data for each particular relative density are determined by way of a twostage process, conducted over the full set of piles in each calibration set. A final, third, stage is employed to determine the calibration parameters for the GDSM. These calibration procedures are described below and summarised in Fig. 7.
Example data, for pile C4; D_{R} = 75%, on the normalised distributed lateral load numerical soil reaction curves at selected depths, z, are shown in Fig. 8(a). At shallow depths, (i.e. z/D = 0·23 and 1·08) where the displacements are relatively large, a peak, followed by postpeak softening, is apparent in the numerical curves. This behaviour is likely to be associated with the dilation characteristics of the soil as represented in the calibration analyses, and was typically observed in the distributed lateral load numerical data. Since softening cannot be represented with the selected conic function, a simplified representation is adopted. At greater depths (e.g. z/D = 2·33 and 5·97 in Fig. 8(a)) the soil reaction curve does not reach a peak. The following calibration process is adopted.
(a)  The value of the ultimate normalised lateral load,  
(b)  The initial stiffness k_{p} is determined by proportional leastsquares fitting the linear expression  
(c)  The curvature parameter, n_{p}, is determined by minimising the proportional leastsquare error between the numerical data and the conic function, for the full range of the data. 
Linear depth variation functions determined by leastsquares fitting to these data are also indicated in Fig. 9. Although more complex depth variation functions could be employed, the overall pile behaviour can be captured remarkably well using just a simple linear fit, as discussed in the paper by Burd et al. (2017) and further demonstrated later in this paper. This is in spite of the significant variability in the individual soil reaction curve parameters. Also shown in Fig. 9, for comparison purposes, are the depth variation functions determined using the final GDSM calibration.
An example set of numerical distributed moment soil reaction curves, for pile C4; D_{R} = 75%, is shown in Fig. 8(b). The response typically tends to a limiting value after a sharp initial rise. At shallow depths a peak is observed in the response. A bilinear form of the parametric curve (n_{m} = 0) was selected in this case; only two parameters therefore require calibration, as follows.
(a)  A high value of initial stiffness is chosen, arbitrarily, as k_{m} = 20.  
(b)  The ultimate normalised moment, 
The base horizontal load numerical soil reaction curves extracted from all of the calibration analysis for D_{R} = 75% are shown in Fig. 10(a). Soil reaction curve parameters are determined as follows.
(a)  The initial stiffness k_{H} is selected by proportional leastsquares fitting the expression  
(b)  A displacement  
(c)  The ultimate displacement  
(d)  The curvature parameter, n_{H}, is determined by minimising the proportional leastsquare error between the numerical data and the conic function. 
(a)  The initial stiffness, k_{M}, is calculated using proportional leastsquares regression for  
(b)  A value of ultimate rotation parameter is selected, arbitrarily, at  
(c)  The curvature parameter, n_{H}, and the ultimate response parameter, 
The parameters from the firststage calibration were used as initial values for this optimisation process. All parameters were allowed to vary by up to ± 50% of their initial value, subject to an upper limit of 1·0 on the curvature parameters, and the need for the soil reaction curve parameters to be nonnegative at all pile locations.
The form of the depth variation functions developed during this process is indicated in Table 5. These functions require the specification of a total of 22 parameters; parameter values determined for D_{R} = 75% at the end of this second stage (stage 2) are listed in Table 5. Values calculated at the ground surface and at the base of a pile of length L/D = 2 and L/D = 6, for D_{R} = 75%, are also tabulated.

Soil reaction component  Soil reaction parameter  Depth variation functions  Depth variation parameters for D_{R} = 75%  Value at surface  Value at base of L/D = 2 pile  Value at base of L/D = 6 pile 

Distributed lateral load, p  Ultimate displacement,  64·78  64·78  64·78  
Initial stiffness, k_{p}  8·64  7·02  3·78  
Curvature, n_{p}  n_{p}  n_{p} = 0·966  0·966  0·966  0·966  
Ultimate reaction,  20·86  15·03  15·03  
Distributed moment, m  Ultimate rotation,  n/a  n/a  n/a  n/a  
Initial stiffness, k_{m}  k_{m}  k_{m} = 18·1  18·1  18·1  18·1  
Curvature, n_{m}  n_{m}  n_{m} = 0·0  0·0  0·0  0·0  
Ultimate moment,  0·23  0·18  0·18  
Base horizontal force, H_{B}  Ultimate displacement,  n/a  1·51  0·27  
Initial stiffness, k_{H}  n/a  2·54  1·06  
Curvature, n_{H}  n/a  0·714  0·482  
Ultimate reaction,  n/a  0·49  0·21  
Base moment, M_{B}  Ultimate rotation,  n/a  49·4  49·4  
Initial stiffness, k_{M}  k_{M}  k_{M} = 0·30  n/a  0·30  0·30  
Curvature, n_{M}  n_{M}  n_{M} = 0·86  n/a  0·86  0·86  
Ultimate reaction,  n/a  0·29  0·09 
Depth variation parameters for D_{R} = 75% determined from the stage 2 optimisation, with values for selected cases also listed.
The GDSM employs simple functions – linear and constant – to represent the dependency of each depth variation parameter on relative density; these are referred to as ‘relative density functions’. If linear functions were to be adopted for all of the (22) model parameters, then a total of 44 relative density parameters would require calibration. It is desirable, therefore, to reduce the calibration space by assigning at least some of the relative density parameters to be constant.
The relative density function forms were chosen in two stages. Initially (stage 3a) the m, H_{B} and M_{B} components were considered. (The relative density functions for the distributed lateral load – the dominant reaction component in terms of overall pile response – were determined in a subsequent process.) Depth variation parameters, from stage 2, for m, H_{B} and M_{B} were inspected. Some of the depth variation parameters – for example, the parameter n_{M} plotted in Fig. 11(a) – indicated a dependency on relative density; linear relative density functions were assigned to these parameters. In other cases – for example, the parameter k_{M} plotted in Fig. 11(b) – where no obvious trend was apparent, constant relative density functions were assigned. An initial set of calibrated relative density functions for the m, H_{B} and M_{B} components, based on these chosen relative density function forms, was then determined by leastsquares fitting to the stage 2 data.
In a subsequent stage (stage 3b) choices were made on the relative density forms for the distributed lateral load. This was done by redetermining the individual depth variation parameters for the distributed lateral load only, for each reference relative density, by minimising the cost function in equation (10). These computations employed the relative density functions from stage 3a to define the model parameters for m, H_{B} and M_{B}. It was discovered that this process reduced the scatter in the distributed lateral load depth variation parameters and therefore facilitated the selection of appropriate relative density function forms for this component. Distributed lateral load parameters which, at this stage, exhibited a consistent dependency on relative density (e.g. the parameter n_{p}), Fig. 12(a), had a linear function of relative density assigned to them. The one parameter that did not exhibit an obvious trend (k_{p2} shown in Fig. 12(b)) was assigned to be a constant.
The system of relative density functions developed in this way is specified by a set of 39 parameters. A final optimisation (stage 3c) was conducted over all of these parameters to minimise the cost function in equation (10). The relative density functions employing this final set of GDSM parameters are specified in Table 6. Note that the fitting process across the relative densities leads to marginal differences between the evaluation of the functions in Table 6 (stage 3c) and the results shown in Table 5 (stage 2).

Soil reaction component  Soil reaction parameter  Relative density functions 

Distributed lateral load, p  Ultimate displacement,  
Initial stiffness, k_{p}  k_{p1} = 8·731 − 0·6982D_{R}  
k_{p2} = −0·9178  
Curvature, n_{p}  n_{p} = 0·917 + 0·06193D_{R}  
Ultimate reaction,  
Distributed moment, m  Ultimate rotation,  Given by 
Initial stiffness, k_{m}  k_{m} = 17·00  
Curvature, n_{m}  n_{m} = 0·0  
Ultimate moment,  
Base horizontal force, H_{B}  Ultimate displacement,  
Initial stiffness, k_{H}  k_{H1} = 6·505 − 2·985D_{R}  
k_{H2} = −0·007969 − 0·4299D_{R}  
Curvature, n_{H}  n_{H1} = 0·09978 + 0·7974D_{R}  
n_{H2} = 0·004994 − 0·07005D_{R}  
Ultimate reaction,  
Base moment, M_{B}  Ultimate rotation,  
Initial stiffness, k_{M}  k_{M} = 0·3515  
Curvature, n_{M}  n_{M} = 0·300 + 0·4986D_{R}  
Ultimate reaction,  
Note: In these relative density functions, the value of D_{R} is expressed as a decimal (i.e. D_{R} = 0·75 for sand with 75% relative density). The relative density functions relate to the depth variation function forms specified in Table 5. The relative density functions are specified in the table to a precision of four significant figures; parameters with this precision were adopted in the 1D model computations described in the current paper. This relatively precise form of the data, selected to be suitable for numerical computations, should not be interpreted as being indicative of the perceived accuracy of these expressions. For a general consideration of the trends and characteristics of the soil reaction curves, employing the data at a lower level of precision (e.g. two significant figures) might be more appropriate.
An indicative convergence study has been conducted in which the 1D (GDSM) model (i.e. a form of the 1D model in which the soil reaction curves are determined by the GDSM) is employed with D_{R} = 75%, for piles C1 and C4, to investigate the sensitivity of the results to the size of the embedded pile and soil elements employed in the model. Calculations were conducted for embedded element lengths of between 0·1 m and 10 m for C1 (L = 20 m) and between 0·5 m and 20 m for C4 D_{R} = 75%. Computed values of the lateral loads η_{sd} and H_{sd} at v_{G} = D/10 and v_{G} = D/10 000, respectively, are listed in Table 7; this table also lists errors in the computed lateral load relative to the finest mesh used in each case.

Number of embedded elements  Embedded element size: m  H_{ult}: MN  Error_{ult}: %  H_{sd}: MN  Error_{sd}: % 

C1  
200  0·1  25·5510  0·0000  0·5384  0·0000 
40  0·5  25·5507  −0·0012  0·5385  0·0054 
20  1  25·5491  −0·0077  0·5386  0·0209 
10  2  25·5529  0·0074  0·5388  0·0739 
4  5  25·6273  0·2986  0·5410  0·4812 
2  10  25·6202  0·2707  0·5487  1·8962 
C4  
120  0·5  174·3406  0·0000  0·7556  0·0000 
60  1  174·3434  0·0016  0·7557  0·0119 
24  2·5  174·3623  0·0124  0·7562  0·0902 
12  5  174·4155  0·0429  0·7583  0·3591 
6  10  174·7061  0·2096  0·7666  1·4576 
4  15  174·8272  0·2791  0·7873  4·1995 
3  20  175·2942  0·5470  0·8153  7·9028 
The results indicate that H_{ult} is remarkably tolerant of employing a relatively coarse mesh for both piles. In all cases, even for the coarsest meshes, the error is less than 1%. The small displacement response appears more sensitive to element size, however. In this case, for both piles, embedded element lengths of 5 m or less are required to achieve an error of less than 1%.
The process conducted to calibrate the GDSM employed a standard embedded element length of 2·5 m. This convergence study suggests that modelling errors associated with mesh discretisation effects in the model calibration process are likely to be negligible.
The H–v_{G} performances of piles C1 (D = 10 m, L = 20 m) and C4 (D = 10 m, L = 60 m) computed using the 1D (GDSM) model for D_{R} = 75% are shown in Fig. 13. A close fit is obtained between the 1D model and the calibration data. The numerical soil reaction curves, together with the parametric curves determined using the GDSM, for the distributed lateral load, are plotted in Fig. 14(a) (for the full range of displacements) and in Fig. 14(b) (for small displacements). It is clear that differences exist between the two sets of data. Although the GDSM soil reaction curves are tailored to provide a representation of the 3D finiteelement data across the complete set of calibration analyses, they can exhibit a tendency, apparent in Fig. 14, to depart from the 3D calibration data for individual piles at a local level. Experience from the use of the 1D model indicates, however, that it is able to reproduce the overall behaviour of the calibration piles to a high accuracy, although at a local level, significant differences can exist between the calibration data and the parametric soil reaction curves.
Table 8 provides the performance metrics for the application of the GDSM to the full range of calibration piles, showing an excellent fit of the model to the data.

Relative density  Pile reference  ρ_{sd}  ρ_{ult}  η_{sd}  η_{ult} 

D_{R} = 45%  C1  0·96  1·09  0·96  0·93 
C2  0·93  1·09  0·93  0·93  
C4  1·03  1·04  0·97  0·95  
C5  0·99  1·03  0·99  0·96  
C6  1·02  0·98  0·98  0·98  
C8  1·09  1·00  0·91  0·98  
C9  1·05  0·99  0·94  0·99  
C10  0·98  1·05  0·98  0·96  
C11  1·05  1·03  0·95  0·96  
D_{R} = 60%  C1  0·94  1·08  0·95  0·96 
C2  0·92  1·07  0·92  0·96  
C4  1·01  1·03  0·99  0·97  
C5  0·98  1·02  0·98  0·98  
C6  1·00  0·93  0·99  0·91  
C8  1·07  0·98  0·92  0·99  
C9  1·04  0·96  0·95  0·98  
C10  0·96  1·02  0·97  0·98  
C11  1·04  1·01  0·96  0·98  
D_{R} = 75%  C1  0·93  1·02  0·94  0·97 
C2  0·91  1·02  0·91  0·96  
C3  0·94  1·02  0·95  0·96  
C4  1·00  1·00  1·00  0·99  
C5  0·97  0·99  0·97  0·99  
C6  0·99  0·92  0·99  0·90  
C7  1·00  0·89  0·98  0·88  
C8  1·07  0·95  0·92  0·98  
C9  1·03  0·94  0·96  0·97  
C10  0·95  0·98  0·96  0·94  
C11  1·03  0·98  0·97  0·99  
D_{R} = 90%  C1  0·91  1·04  0·92  0·98 
C2  0·89  1·04  0·89  0·99  
C4  0·99  1·02  0·99  0·95  
C5  0·95  1·01  0·95  0·97  
C6  0·97  0·94  0·98  0·95  
C8  1·05  1·00  0·94  0·97  
C9  1·02  0·99  0·97  0·98  
C10  0·93  1·00  0·95  0·99  
C11  1·01  1·01  0·98  0·96  
Average  0·99  1·00  0·96  0·96  
CoV  5·14%  4·52% 
To demonstrate the predictive capability of the 1D (GDSM) model, various design examples have been considered. The geometries of these example cases, specified in Table 9, are selected to fall within the calibration space but to differ from the geometric conditions employed for the calibration piles. Values of relative density have been chosen that fall within the calibration space but not at the original calibration densities.

Pile reference  D: m  h: m  h/D  L: m  h/L  L/D  t: mm  D/t 

D1  7·5  37·5  5  22·5  1·67  3  68  110 
D2  8·75  87·5  10  35  2·5  4  91  96 
D2t  8·75  87·5  10  35  2·5  4  150  58 
The load–displacement responses computed for pile D2t using the 1D (GDSM) model and, separately, with corresponding 3D finiteelement models, are shown in Fig. 15, for relative densities 55% and 85%. A close match is observed between the two data sets. Fig. 15 also shows excellent agreement of the bending moments induced in the embedded portions of the piles, determined for D_{R} = 85% (where H_{ult} is the lateral load determined from the 3D finiteelement analysis at v_{G} = 0·1D) and also for H = 0·5H_{ult}.
Values of accuracy and ratio metrics for a set of 13 design example cases are listed in Table 10. These data, which indicate a close match between the 1D (GDSM) model and corresponding 3D finiteelement results, support the assumption implicit in the PISA methodology, that the 1D model provides an efficient means of interpolating the overall pile response computed using the 3D calibration calculations to other pile geometries and relative densities within the calibration space.

Relative density  Pile reference  ρ_{sd}  ρ_{ult}  η_{sd}  η_{ult} 

D_{R} = 45%  D1  0·97  1·04  0·97  0·96 
D2  1·03  1·06  0·97  0·94  
D2t  1·07  1·06  0·93  0·94  
D_{R} = 60%  D1  0·95  0·99  0·96  0·99 
D2  1·02  1·01  0·98  0·98  
D_{R} = 75%  D1  1·03  0·93  0·96  0·95 
D2  1·01  0·96  0·99  0·98  
D_{R} = 90%  D1  0·92  0·96  0·93  0·97 
D2  0·99  1·00  0·99  0·96  
D2t  1·02  0·99  0·98  0·96  
D_{R} = 55%  D2t  1·05  1·04  0·94  0·96 
D_{R} = 70%  D2t  1·04  0·97  0·95  0·98 
D_{R} = 85%  D2t  1·02  0·97  0·97  0·98 
Average  1·01  1·00  0·96  0·97  
CoV  4·13%  4·15% 
The PISA design model provides a rapid means of conducting design calculations for monopile foundations for offshore wind turbines. This paper demonstrates an application of the model to homogeneous marine sand sites, complementing the modelling approach described in the paper by Byrne et al. (2020) for glacial clay till soils. The model is capable of delivering predictions of performance that closely match the results obtained from equivalent 3D finiteelement models.
The paper describes a calibration process based on the soil conditions at the PISA sand site in Dunkirk. This calibration is considered to provide a realistic model for monopiles installed at offshore sand sites where the characteristics of the sand are similar to the Flandrian sand encountered at Dunkirk and where the monopile dimensions fall within the calibration space. In other cases, application of the model may require a separate calibration exercise. The model has been demonstrated for monopiles with uniform wall thickness. However, the model can be applied straightforwardly, to piles with variations in wall thickness along their embedded length, by the specification of appropriate structural properties for the beam elements in the 1D model. The normalisations employed in the model do not explicitly include the load eccentricity, h, although the optimised calibration parameters are likely to depend on the range of h/D employed in the calibration process. It therefore follows that the model should not be used for values of h/D (or indeed any other pile parameters) that fall outside the calibration space.
The PISA design model is shown to reproduce the overall behaviour of the calibration piles, even though at a local level significant differences can exist between the numerical soil reaction curves and the calibrated model. This apparently wellconditioned aspect of the model is considered to be due to the overall pile performance being obtained by integrating the soil reaction curves along the entire length of the foundation. Provided that significant systematic errors are absent, this averaging process appears to have the consequence that the model is remarkably tolerant of imperfect fitting of the data at a local level.
The approximate nature of the Winkler modelling approach adopted in the PISA design model has a number of implications. First, as is indicated by the considerable scatter in the data in Fig. 9, the model is unable to represent the pile–soil interaction at all points along the pile in a highfidelity manner. The stage 1 calibration process is modestly successful at representing the overall monopile performance, but the performance of the model was found to be enhanced by the use of a further optimisation process stage 2. Although the stage 2 process (and to an extent the stage 3 process) improves the overall performance of the model, it does not necessarily lead to an improved representation of the actual physics of the local soil–pile interaction. Instead, the stage 2 and stage 3 optimisation should be understood as a pragmatic expedient to calibrate an imperfect model (Winkler) to provide highfidelity predictions of behaviour within a predefined calibration space. It is also necessary to recognise that any modelling errors inherent in the 3D finiteelement calibration analyses will be inherited by the design model.
The current form of the PISA design model is restricted to monotonic loading. Extensions to cyclic loading are feasible – for example by the development of cyclebycycle soil reaction curves, or the implementation of approaches in which the (monotonic) soil reaction curves are modified to reflect the influence of previous load cycling. The model is demonstrated for homogeneous soil deposits only, whereas offshore sites usually consist of layered profiles, often involving interbedded clays and sands. This can be addressed using the PISA design model by assigning clay soil reaction curves (Byrne et al., 2020) to the clay layers and employing the current model for the sand layers; Byrne et al. (2019a) describe an initial evaluation of this approach.
ACKNOWLEDGEMENTS
The PISA Phase 1 project was funded by the UK Department for Energy and Climate Change (DECC) and the PISA industry partners under the umbrella of the offshore wind accelerator (OWA) programme, which was designed and is led by the Carbon Trust. The authors acknowledge the provision of financial and technical support by the following project partners: Ørsted Wind Power (formerly: DONG Energy), Alstom Wind, E.ON, EDF, Equinor (formerly: Statoil), innogy, SPR, Statkraft, SSE, Vattenfall and Van Oord. Development of the specific GDSM aspects of the PISA design model were conducted as part of the PISA Phase 2 project, which was supported by the Scottish Government and all partners from the PISA project except Statkraft.
D  constitutive matrix 
D_{t}  tangent constitutive matrix 
D  pile diameter 
D_{R}  soil relative density 
e_{0}  initial soil void ratio 
G_{0}  smallstrain soil shear modulus 
H  lateral load applied to pile 
H_{1D}  lateral load applied to pile, computed with the 1D model 
H_{3D}  lateral load applied to pile, computed with the 3D model 
H_{B}  horizontal force at pile base 
H_{G}  lateral load applied to pile at ground level 
H_{sd}  lateral load applied to pile at a groundlevel displacement of v_{G} = D/10 000 
H_{ult}  lateral load applied to pile at a groundlevel displacement of v_{G} = D/10 
h  load eccentricity 
k  initial stiffness of parametric soil reaction curve 
k  soil finiteelement stiffness matrix 
L  pile embedded length 
L_{e}  length of soil finite element 
M_{B}  moment at pile base 
M_{G}  moment applied to pile at ground level 
m  distributed moment acting on monopile 
N_{i}  conventional set of Hermite cubic interpolation polynomials 
n  curvature parameter for parametric soil reaction curve 
p  distributed lateral load applied to pile 
p  local distributed load/moment vector 
s_{u}  undrained shear strength of soil 
t  pile wall thickness 
U  vector containing element displacement degrees of freedom 
V_{1}, V_{2}  lateral displacement values 
v  lateral pile displacement 
v  local displacement/rotation vector 
v_{B}  lateral pile displacement at pile base 
v_{G}  groundlevel lateral pile displacement 
ultimate displacement for parametric soil reaction curve  
ultimate load for parametric soil reaction curve  
Z_{1}, Z_{2}  depth coordinates of nodes 1 and 2 
z  depth coordinate along the pile 
α  normalised depth coordinate 
η  accuracy metric 
η_{sd}  small displacement accuracy metric (groundlevel pile displacements up to D/10 000) 
η_{ult}  ultimate displacement accuracy metric (groundlevel pile displacements up to D/10) 
ρ  ratio metric 
ρ_{sd}  small displacement ratio metric (groundlevel pile displacements up to D/10 000) 
ρ_{ult}  ultimate displacement ratio metric (groundlevel pile displacements up to D/10) 
σ_{vi}  initial vertical effective stress in the soil 
Ψ_{1}, Ψ_{2}  pile crosssection rotation at nodes 1 and 2 
ψ  rotation of the pile crosssection 
ψ_{B}  rotation of the pile crosssection at the pile base 
Discussion on this paper closes on 1 March 2021, for further details see p. ii.