PISA design model for monopiles for offshore wind turbines: application to a marine sand

This paper describes a one-dimensional (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 co-workers; 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 three-dimensional finite-element 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.


INTRODUCTION
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 non-linear 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 projectknown as PISAthat included field testing 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 three-dimensional (3D) finite-element modelling 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 finite-element calibration analyses; it therefore benefits from the realism that is potentially achievable with 3D finite-element 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  and the finite-element 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 modelreferred 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.

Model overview
The PISA design model provides a one-dimensional (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 interfacefor example, due to local rotation of the pile cross-section. 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 four-component 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 consists of a circular tube with outer diameter, D, wall thickness, t, and embedded length, L. The height, h, of the load application is referred to as load eccentricity L/D is reduced (Byrne et al., , 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 four-component 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).

Soil reaction curves for the sand modelling framework
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 small-strain shear modulus, and v and ψ are the local pile lateral displacement and cross-section rotation, respectively.
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).
The current implementation of the PISA design model employs the four-parameter conic function wherex signifies a normalised displacement or rotation variable andȳ signifies the corresponding normalised soil reaction component, formulated in terms of the   Table 1. The function is illustrated in Fig. 3. The normalised soil reactions can be determined explicitly from the normalised displacements bȳ where Each of the parameters ðk;ȳ u ;x u ; nÞ has a straightforward interpretation. The parameter k specifies the initial slope;ȳ u is the ultimate value of the normalised soil reaction; andx u is the normalised displacement (or rotation) at which this ultimate value of soil reaction is reached. The parameter n (0 n 1) determines the shape of the curve; for the extreme values n ¼ 0 and n ¼ 1, the function reduces to the bilinear forms illustrated in Fig. 3(b). As discussed in the paper by Byrne et al. (2020), no particular importance is attached to the specific form that is chosen for the parametric curve, and other similar parametric functions would be possible.
For the distributed lateral load, and the base horizontal load and moment components, dimensional forms of the soil reaction curves, and their derivatives, for numerical implementation of the model, are determined in a straightforward way using local values of G 0 and σ′ vi , on the basis of the dimensionless forms in Table 1. The particular normalisation adopted for the distributed moment, however, means that a different treatment is required in this case. During the initial model development process, data from the 3D calibration calculations suggested that the distributed moment, m, appeared to scale with the current value of the local distributed lateral load, p. Since the vertical tractions induced on the pile perimeter arise as a consequence of friction at the soil-pile interface, it seems plausible that the magnitude of the distributed moment correlates with the local normal tractions; in turn, these tractions are closely related to the local distributed lateral load. It was therefore decided to adopt a dimensionless form for the distributed moment,m, by normalising the distributed moment by the local value of the distributed load, as indicated in Table 1. The use of this form form implies that the distributed moment is a function of both the local displacement, v, and the local pile cross-section rotation, ψ This coupling has certain implications for the numerical implementation of the model (described in the Appendix).
1D finite-element formulation for the sand modelling framework The PISA design model employs the 1D representation finite-element framework illustrated in Fig. 2(b). The pile is represented by a line mesh of two-noded 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 finite-element calibration calculations, the structural properties (area and second moment of area) are specified for the beam elements using the thin-walled 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 finite-element 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.

REPRESENTATIVE OFFSHORE SAND SITES
In connection with the PISA research, a series of pile tests  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 ; 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 finite-element study  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  were employed directly to characterise the representative ground models for the current study; the only parameter requiring adjustment is the relative density.
The relative density for the natural Flandrian sand at Dunkirk was estimated as D R ¼ 75%; corresponding to an initial void ratio, e 0 , of 0·629. This relative density was adopted for one of the representative offshore ground models for the current study. Three additional representative ground models were developed, with D R ¼ 45%, 60% and 90%; see Table 2. The initial stresses were determined by adopting hydrostatic pore pressure conditions with a submerged unit weight of γ′ ¼ 10·09 kN/m 3 and K 0 = 0·4 (both values correspond to data for the Dunkirk site). Variations in the submerged unit weight of the sand due to different values of relative density were not considered, as the effect of varying this parameter is regarded as minimal with respect to other aspects of sand behaviour. The small-strain shear modulus, G 0 , is obtained from the local values of mean effective stress p′ and initial void ratio e 0 (from Table 2) using the relationship proposed by Hardin & Black (1968) where B, determined from the site investigation data ) and the calibration process described in the paper by Taborda et al. (2019), is specified in Table 3 and 3D FINITE-ELEMENT CALIBRATION CALCULATIONS Specification of the finite-element calculations The 3D finite-element 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 finite-element software ICFEP (Potts & Zdravković, 1999, 2001. In total, 38 calibration analyses were conducted (see Table 8 later).
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.

Modelling procedures
The 3D finite-element 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 waythat is without the need to adjust the model parameters for soils with different relative Table 2. Initial values of void ratio, e 0 , and relative density for the representative offshore sites Relative density, D R : % Initial void ratio, e 0 45 0·741 60 0·685 75 0·629 90 0·573 Table 3. Constitutive parameters for the sand constitutive model (Taborda et al., 2014) employed in the 3D finite-element calibration analyses

Component Parameters
Critical state line These parameters are identical to those that were determined, as described in the paper by Taborda et al. (2019), to conduct 3D finite-element analysis of the PISA test piles at Dunkirk.
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 20-noded hexahedral displacement-based 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 eight-noded 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 above-ground extension is modelled with 240 shell elements. The interface between the soil and the pile exterior is modelled with 360 16-noded zero-thickness 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 finite-element mesh at the start of the analysis in a fully plugged configuration. A similar wished-in-place procedure was employed in the 3D finite-element models that were developed to analyse the Dunkirk PISA test piles .
The interface between the exterior of the embedded monopile and the soil is represented by an elasto-plastic 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.

NUMERICAL SOIL REACTION CURVES
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 finite-element 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 cross-section rotations of the pile were determined from the computed displacements of the relevant shell element nodes by averaging over the cross-section (for displacement) and by least-squares 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 ground-level 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 finite-element 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.
Quantitative comparisons between the performance of the 1D and 3D models employ the 'accuracy metric', η where the meaning of A ref and A diff are illustrated in Fig. 6. A 'ratio metric', defined by is also employed, where H 1D and H 3D are values of lateral force, computed from the 1D and 3D models, respectively, at particular values of v G , as shown on Fig. 6. The accuracy metric, η, evaluates the precision of the overall fit (and is expected to be close to 1), while the ratio metric, ρ, indicates whether the model under-predicts (,1) or over-predicts (.1) 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.

PARAMETRIC SOIL REACTION CURVES Selection and calibration of the parametric soil reaction curves
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 four-parameter 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 two-stage 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.
The conic function employed to represent the soil reactions is intended forx;ȳ ! 0 only (i.e. in the positive quadrant). Depending on the direction of the applied load and the adopted sign convention, values ofx;ȳ extracted from the calibration analyses may be negative. Also, for the distributed load and moment components, the direction ofx;ȳ may vary with position along the pile. The process of fitting the conic function to the numerical data is conducted by first mapping all of the numerical data into the positive quadrant. In the subsequent implementation in the 1D finite-element model, the soil reaction curves for the full range ofx (positive and negative) are specified on the basis that the response in the third quadrant (x;ȳ , 0) is identical to that in the first quadrant, but with appropriate sign changes.
First-stage calibration Distributed lateral load soil reaction curves. 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 post-peak 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,p u , is taken as the value of the numerical soil reaction curve at large displacement (i.e. the final increment of the analysis). For softening behaviour, this value ofp u is initially reached earlier in the analysis; in this case the ultimate normalised displacementv pu is selected as the value at the first increment of the numerical soil reaction curve at whichp u is exceeded. Otherwisev pu is taken as the value at the final analysis increment.   (b) The initial stiffness k p is determined by proportional least-squares fitting the linear expressionp ¼ k pv to the numerical soil reaction curve for 0 ,p , 0Á1. (c) The curvature parameter, n p , is determined by minimising the proportional least-square error between the numerical data and the conic function, for the full range of the data.
Distributed lateral load parameters determined for all of the piles in the calibration set, for D R ¼ 75%, are plotted in Fig. 9. To develop functions (referred to as 'depth variation functions') to represent the dependency of the parameters on depth, z, it is convenient to employ normalised depth parameters that collapse the data (as far as possible) onto a single variable. Adopting a normalised depth z/D for k p and n p , and an alternative normalised depth z/L forv pu andp u , appeared to provide the best approach; the parameters are plotted with respect to these normalised depth variables in Fig. 9 for all of the piles in the calibration set for D R ¼ 75%.
The data in Fig. 9 exhibit a certain amount of variability and scatter along the pile. Some of these patterns can be related directly to physical aspects of the problem. For example, the cluster of points in Fig. 9(a) with relatively high values of k p close to z/D ¼ 2 all relate to the short monopiles L/D ¼ 2 employed in the calibration set. It appears that these short, relatively stiff, monopiles attract a larger lateral soil stiffness near their base than the more flexible L/D ¼ 6 piles. The apparent discontinuity in the k p data close to z/D ¼ 2·6 is associated with the behaviour of the L/D ¼ 6 piles near the pivot point (where the direction of the lateral displacements changes sign with increasing depth). A similar influence of the pivot point (the location of the pivot tends to increase in depth as displacements increase) is seen in the data in Fig. 9(c). Other features of Fig. 9 relate to the calibration process. For example, the soil near the base of the L/D ¼ 6 piles was not taken to failure in the calibration analyses (since the lateral displacements induced near the base of the piles were relatively small). As a consequence, thep u data in Fig. 9(c) for relatively large depths seem unrealistically low. This is actually of little consequence for the PISA design model, since the model provided is only used within the calibration space, soil failure will not be approached near the pile base in any design calculations. A further aspect of the data relates to the actual physics of the problem being represented by an imperfect (Winkler) model. It is assumed in the current model, for example, that the lateral distributed load depends only on the lateral displacement, but there is likely also to be a dependency on local rotation. Additionally, the data are normalised with respect to the local soil stiffness and strength; the actual lateral distributed load as determined from the finite-element analysis doubtless depends on non-local spatial stiffness/strength variations. Moreover, the spatial coupling within the soil is ignored. The influences of these various approximations will be likely to vary with the dimensions of the pile and the loading eccentricity. These factors combine to generate the significant scatter observed in Fig. 9 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.
Distributed moment curves. 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. Base horizontal load curves. 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 least-squares fitting the expressionH B ¼ k HvH to the numerical data for 0 ,H B , 0Á01. (b) A displacementv HðHÀmaxÞ is established at which the peak value ofH B is first reached. The normalised ultimate response parameter,H Bu , is calculated as the average of the normalised base horizontal force values forv H .v HðHÀmaxÞ . (c) The ultimate displacementv Hu is selected as the first normalised displacement at which the normalised numerical soil reaction is equal toH Bu . (d ) The curvature parameter, n H , is determined by minimising the proportional least-square error between the numerical data and the conic function.
Base moment curves. The base moment reaction curves extracted from the calibration analyses for D R ¼ 75% are shown in Fig. 10(b). Soil reaction curve parameters are determined as follows.
(a) The initial stiffness, k M , is calculated using proportional least-squares regression for 0 ,ψ M , 0Á05. (b) A value of ultimate rotation parameter is selected, arbitrarily, atψ Mu ¼ 50. This value exceeds the computed normalised rotations and allows reasonable values of the curvature parameter to be selected. (c) The curvature parameter, n H , and the ultimate response parameter,M Bu , are selected by minimising the proportional least-square error between the numerical data and the conic function.
It is seen from the above that threshold values for the distributed lateral load (p , 0Á1), base horizontal force (H B , 0Á01) and base moment (M B , 0Á05) were adopted to determine the relevant initial stiffness parameters. These threshold values are essentially arbitrary and were selected for the current work, on the basis of experimentation, to ensure a satisfactory match between the finite-element calibrations and the calibrated 1D model, for small displacements.  C1  C2  C3  C4  C5  C6  C7  C8  C9  C10  C11 First stage GDSM Fig. 9. Depth variations of the normalised distributed load soil reaction curve parameters for D R = 75%. The markers indicate data determined from stage 1 for all of the calibration piles; also shown are regression lines that are fitted to these data. For comparison purposes, depth variation functions corresponding to the GDSM are also indicated: (a) initial stiffness parameter, k p ; (b) ultimate displacement parameter,ˉv u ; (c) ultimate response parameter,ˉp u ; (d) curvature parameter, n p

Second-stage optimisation
To improve the fit between the 1D model and the 3D finite-element calibration data, adjustments are made to the depth variation function parameters to minimise the cost function, k p2 ¼ À0·9178 where η ult,i and η sd,i are the ultimate and small displacement accuracy metrics, respectively, for pile C i , i ¼ 1:N and the summation is taken over the piles in the calibration set (N ¼ 11 for D R ¼ 75% and N ¼ 9 for the other relative density cases). This process was conducted, separately, for each relative density using optimisation routines implemented in Matlab.
The parameters from the first-stage 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 non-negative 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.

Third stage optimisation; relative density functions
The GDSM employs simple functionslinear and constantto 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 loadthe dominant reaction component in terms of overall pile responsewere 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 parametersfor 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 casesfor 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 least-squares 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 re-determining 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) 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).

Convergence study
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. Table 5. General forms of depth variation functions employed in the sand modelling framework, calibrated within the parameter space set out in Table 4 Soil reaction component Ultimate rotation,ψ muψmu n/a n/a n/a n/a Depth variation parameters for D R ¼ 75% determined from the stage 2 optimisation, with values for selected cases also listed. 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.

ANALYSIS OF THE CALIBRATION CASES USING THE GDSM
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 Ultimate rotation,ψ mu Given bym u =k m Initial stiffness, k m 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. curves are tailored to provide a representation of the 3D finite-element 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.

DESIGN EXAMPLES
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.
The load-displacement responses computed for pile D2t using the 1D (GDSM) model and, separately, with corresponding 3D finite-element 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 finite-element 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 finite-element 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.

DISCUSSION AND CONCLUDING REMARKS
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 finite-element 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 well-conditioned 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 high-fidelity 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 high-fidelity predictions of behaviour within a predefined calibration space. It is also necessary to recognise that any modelling errors inherent in the 3D finite-element 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 feasiblefor example by the development of cycle-by-cycle 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  The dependence of distributed moment on the distributed lateral soil reaction in the sand modelling framework requires some special consideration for the numerical implementation of the PISA design model. The finite-element formulation adopted for the model employs soil finite elements that have compatible displacements with the finite elements used to model the pile (see Byrne et al. (2020) for further details). A soil finite element of length L e connected to nodes 1 and 2 with depth coordinates Z 1 and Z 2 , respectively, is illustrated in Fig. 16.
The lateral displacement within this soil element is defined by the interpolation where V 1 , V 2 are the values of lateral displacement and Ψ 1 , Ψ 2 are the pile cross-section rotation at nodes 1 and 2, respectively. The functions N i are the conventional set of Hermite cubic interpolation polynomials, given by N 2 ¼ αL e ð1 À 2α þ α 2 Þ N 3 ¼ 3α 2 À 2α 3 N 4 ¼ αL e ðÀα þ α 2 Þ where α ¼ (z À Z 1 )/L e and z is a general depth coordinate of a point within the element. The shear strain, γ 0 , is assumed constant within each beam element; since the soil and beam elements share the same interpolation, the shear strain γ 0 in the beam elements also appears as a degree of freedom for the soil elements in equation (11). The local displacement v and rotation ψ within each soil element is given by where v ¼ ð v ψ Þ T is the local displacement/rotation vector, U ¼ ð V 1 Ψ 1 γ 0 V 2 Ψ 2 Þ T is a vector containing the element displacement degrees of freedom and where N′ i denotes a shape function derivative with respect to z. The stiffness matrix k for the soil elements is determined from where D is an appropriate constitutive matrix. The tangent constitutive matrix, D t , for the soil element is η ult ultimate displacement accuracy metric (ground-level pile displacements up to D/10) ρ ratio metric ρ sd small displacement ratio metric (ground-level pile displacements up to D/10 000) ρ ult ultimate displacement ratio metric (ground-level pile displacements up to D/10) σ vi initial vertical effective stress in the soil Ψ 1 , Ψ 2 pile cross-section rotation at nodes 1 and 2 ψ rotation of the pile cross-section ψ B rotation of the pile cross-section at the pile base