PISA design model for monopiles for offshore wind turbines: application to a stiff glacial clay till

Offshore wind turbines in shallow coastal waters are typically supported on monopile foundations. Although three-dimensional (3D) finite-element methods are available for the design of monopiles in this context, much of the routine design work is currently conducted using simplified one-dimensional (1D) models based on the p – y method. The p – y method was originally developed for the relatively large embedded length-to-diameter ratio ( L / D ) piles that are typically employed in offshore oil and gas structures. Concerns exist, however, that this analysis approach may not be appropriate for monopiles with the relatively low values of L / D that are typically adopted for offshore wind turbine structures. This paper describes a new 1D design model for monopile foundations; the model is specifically formulated for offshore wind turbine applications, although the general approach could be adopted for other applications. The model draws on the conventional p – y approach, but extends it to include additional components of soil reaction that act on the pile. The 1D model is calibrated using a set of bespoke 3D finite-element analyses of monopile performance, for pile characteristics and loading conditions that span a predefined design space. The calibrated 1D model provides results that match those obtained from the 3D finite-element calibration analysis, but at a fraction of the computational cost. Moreover, within the calibration space, the 1D model is capable of delivering high-fidelity computations of monopile performance that can be used directly for design purposes. This 1D modelling approach is demonstrated for monopiles installed in a stiff, overconsolidated glacial clay till with a typical North Sea strength and stiffness profile. Although the current form of the model has been developed for homogeneous soil and monotonic loading, it forms a basis from which extensions for soil layering and cyclic loading can be developed. The general approach can be applied to other foundation and soil – structure interaction problems, in which bespoke calibration of a simplified model can lead to more efficient design. 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 accuracyof the parameters. 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.


INTRODUCTION
Monopile foundations are currently the preferred option for offshore wind turbine structures in shallow coastal waters. The design of monopile foundations in this application requires analysis tools that can deliver reliable predictions of performance for lateral and moment loading conditions due to the action of wind, waves and current. Accurate predictions of the stiffness of the foundation are also required to allow estimates to be made of the natural frequency of the wind turbine support structure.
Current design procedures typically employ a simplified one-dimensional (1D) analysis frameworkknown as the 'p-y' methodin which the monopile is modelled as an embedded beam and the lateral soil response is represented by 'p-y' curves, which are non-linear relationships between the distributed lateral load acting on the pile and the local pile lateral displacement. The p-y method was originally developed for the relatively slender piles used in the offshore oil and gas industry (e.g. Reese & Matlock, 1956;Matlock, 1970;Reese et al., 1975). Standard forms of the p-y method (e.g. as specified in API (2010) and DNV-GL (2016)), however, are widely regarded as having significant limitations when applied to piles with the relatively low L/D ratios of about 6 or less (where L is embedded length and D is pile diameter) that are typically employed for offshore wind turbine applications (e.g. Doherty & Gavin, 2011;Byrne et al., 2017). Particular concerns include (a) evidence of inconsistencies between design predictions based on the p-y approach and measured performance of actual offshore wind turbine structures (e.g. Kallehave et al., 2012) (b) lack of clarity in relating the parameters needed to define the p-y curves to soil data obtained during routine site investigation (e.g. Jeanjean et al., 2017) (c) lack of robust extensions of the p-y method to address cyclic loading effects.
This paper describes a new modelling approachreferred to as the 'PISA design model'for rapid, high-fidelity design calculations for offshore monopile foundations. The model draws on the traditional p-y method, but extends it to include additional soil reaction components that have been identified as significant for piles with relatively low values of L/D (e.g. Byrne et al., 2015Byrne et al., , 2017Byrne et al., , 2019aSchroeder et al., 2015). Importantly, the model is calibrated against a set of bespoke three-dimensional (3D) finite-element analyses.
Once the model has been calibrated it can be used to conduct rapid computations, with results that closely match those obtained using more detailed, but computationally costly, 3D finite-element methods. The PISA design model addresses some of the limitations of the traditional p-y method when applied to monopile foundations (although the current model does not incorporate cyclic effects). In particular, it facilitates accurate predictions for piles with relatively low values of L/D. Furthermore, the 3D finite-element procedures, against which the method is calibrated, have themselves been validated through detailed comparison with field tests Zdravković et al., 2019b). Thus there is an audit trail from field data through 3D analysis and to the new design method. The current application of the PISA design model is concerned with piles embedded in a single soil type, but the approach is capable of extension to layered soils. Future extensions to cyclic loading are feasible based on the computed monotonic backbone curve, but these extensions are beyond the scope of this paper. The development of the design model presented here was the purpose of a recently completed joint industry studyknown as PISA Zdravković et al., 2015;Burd et al., 2017) comprising field testing and ground characterisation Byrne et al., 2019b;McAdam et al., 2019;Zdravković et al., 2019a), 3D finite-element modelling Zdravković et al., 2019b) and 1D model development (described here and in the paper by Burd et al. (2020)). The various aspects of the PISA research were all based on the conventional design assumption that the lateral behaviour of monopiles for wind turbine applications is unaffected by vertical loads (caused by the weight of the turbine and the support structure). This is on the basis that vertical loads applied to monopiles in this application are invariably a small fraction of the vertical load capacity of the foundation.
Field tests and associated numerical modelling were conducted at two test sites; one at Cowden (on the north east coast of England) where the soil consists principally of an overconsolidated glacial till Zdravković et al., 2019b). The other was at Dunkirk in northern France, where the soil is a dense marine sand Taborda et al., 2019). The current paper describes the new modelling framework and provides an example application for a 'representative offshore glacial clay till site', with ground conditions based closely on those at the Cowden test site. A companion paper describes a separate application to 'representative marine sand' sites (Burd et al., 2020).

Model framework
The PISA design model is concerned with computing the performance of a monopile foundation for the loading conditions shown in Fig. 1. A monotonic lateral load H is applied at a height h above the seabed level (referred to in this paper as 'ground level'). The PISA model provides a means of determining the resulting displacements of the pile as well as the distribution of bending moment and shear force in the pile. The distance h, referred to as 'load eccentricity', is chosen to provide an appropriate representation of the environmental loads. Vertical loading applied to the foundation (e.g. due to the weight of the turbine, rotor and tower) is assumed to be small compared to the axial pile capacity, and is neglected in the current model.
The loading configuration in Fig. 1 applies a lateral load H G ¼ H and a moment M G ¼ Hh to the monopile at ground level, as indicated in Fig. 2. Four separate soil reaction components are assumed to act on the pile at the soil-pile interface, as illustrated in Fig. 2(a). These are: (a) distributed lateral loads; (b) vertical shear tractions; (c) a horizontal force at the pile base; and (d ) a moment at the pile base. The vertical shear tractions arise partly as a consequence of vertical displacements of the pile perimeter caused by local rotation of the pile cross-section and partly due to vertical relative movements of active and passive soil wedges that develop near the ground surface as soil failure conditions are approached (see, e.g. Jeanjean et al. (2017)). The shear tractions combine to form a distributed moment along the pile; any net vertical load applied to the pile by the tractions is neglected. In the PISA design model implementation illustrated in Fig. 2(b), the monopile is represented as an embedded beam; a distributed lateral load p and a distributed moment m are assumed to act on the pile along its length, with the distributed moment providing a means of incorporating the moment associated with the vertical shear tractions induced at the soil-pile interface. Additionally, a horizontal force H B and a moment M B act on the pile base.
The PISA design model is formulated in a 1D finite-element framework in which the pile is represented as a line mesh of beam finite elements. The soil is modelled with a separate set of finite elements, with displacement interpolation functions that are identical to those employed in the elements representing the pile. Each embedded pile element has a soil element associated with it, attached to its two nodes. These procedures ensure compatibility of displacement and rotation of the pile and the soil along the length of each embedded pile element. Each of the soil reaction components is related in the model to the local lateral displacement or rotation (i.e. adopting a 'Winkler' approach) by a calibrated parametric function referred to as a 'soil reaction curve'. The soil reaction curves have a similar role to the p-y curves in a conventional p-y analysis, but the PISA model is extended to include a representation of additional soil reaction components. These extensions follow previous work by Davidson (1982), Lam & Martin (1986) and Lam (2013) for the design of drilled shafts, principally for onshore applications. It is acknowledged that the Winkler assumption (soil reactions depend solely on the local displacement and rotation) is an oversimplification of the actual response of the soil body, as it does not represent the spatial coupling that occurs within the soil. The Winkler assumption has, however, a long track record of satisfactory application to soilstructure interaction problems, and has the advantage that it allows highly efficient computational techniques to be employed. The efficacy of adopting the Winkler approach is demonstrated later in this paper by the fidelity that it achieves in reproducing the more complex 3D analyses. Versteijlen et al. (2018) employ a simplified monopile model for elastic soil behaviour, in which spatial coupling within the soil is included. Extension of these ideas to incorporate non-linear pile behaviour would not be straightforward.
A four-parameter function is used as the basis for each of the soil reaction curves. The curves are calibrated directly from a set of bespoke 3D finite-element calibration analyses, tailored to a particular offshore site and range of monopile dimensions and loading eccentricities of interest. This calibration process ensures that the 1D model is able to provide a close representation of the performance of the monopile, as computed with the more detailed 3D analyses, for arbitrary pile dimensions and loading eccentricities within the calibration space.

1D finite-element formulation
Timoshenko beam theory is used to model the behaviour of the embedded pile; providing an approximate means of including the effect of shear deformations in the pile within the analysis. Gupta & Basu (2018) demonstrate that, although for many cases the effect of using the more accurate Timoshenko theory rather than Euler-Bernoulli theory is small, for piles of low L/D ratio and a low ratio of pile to soil stiffness, the effects of using the more accurate theory can result in a significant increase in pile head displacement, and especially of pile head rotation (see e.g. their Fig. 6). The use of Timoshenko beam theory in the current model, therefore, provides assurance that shear effects are captured in any cases where they might be significant.
The assumed displacement field in the pile is indicated in Fig. 3(a). The axial and lateral displacements are where y is distance from the neutral axis (assumed to coincide with the centroid of the pile cross-section); ψ is the rotation of the pile cross-section; and v 0 (z) is the lateral displacement of the neutral axis. The rotation ψ is defined to be clockwise positive, Fig. 3(b), consistent with a positive rotation about a right-handed x-axis implied by the coordinate directions in Fig. 3(a). The corresponding strains are (illustrated in Fig. 3 where E is the Young's modulus of the pile material; I is second moment of area of the pile cross-section; G is the shear modulus; A is the pile cross-section area; and κ is a shear factor. The shear factor employed in the current model is determined as the ratio of the average shear stress acting on the cross-section and the shear stress at the neutral axis for a thin-walled tube, determined using conventional beam theory; this approach gives κ ¼ 0·5. An alternative analysis, Cowper (1966), gives the shear factor as where ν is the Poisson ratio of the beam material; ν ¼ 0·3 in the current work. The shear factor determined from equation (5) is 0·53, similar to the value κ ¼ 0·5 employed in the current model. For equilibrium configurations of the model in Fig. 2(a), for arbitrary virtual lateral displacements and rotations δv and δψ, the internal virtual work, δW I , is equal to the external virtual work, δW E The external virtual work is where δv G is the virtual displacement of the pile at ground level and δψ G is the virtual cross-section rotation of the pile at ground level. The internal virtual work is where δv B and δψ B are virtual displacement and cross-section rotation, respectively, at the base of the pile. In this equation, the soil reactions ( p, m, H B , M B ) are considered as internal force resultants, determined in the model by specified functions of the local displacement and rotation. The Galerkin form of equation (8) is obtained by discretising the pile using a line mesh of two-noded, five-degreesof-freedom Timoshenko beam elements, based on the formulation in the book by Astley (1992). The soil is represented by a line mesh of two-noded, five-degrees-of-freedom elements using displacement and rotation interpolation procedures that are consistent with those employed for the beam elements. The pile and soil meshes share the same nodal degrees of freedom. An overview of the finite-element formulation is provided below.
The lateral displacement within each pile element is determined by the interpolation where N i are the conventional set of Hermite cubic interpolation functions, and V j , Θ j are the nodal values of lateral displacement and pile neutral axis rotation, PISA DESIGN MODEL FOR MONOPILES FOR OFFSHORE WIND TURBINES respectively. The shear strain is assumed constant, γ yz ¼ γ 0 , within each element. The neutral axis rotation at each node is Θ j ¼ γ 0 À Ψ j , where Ψ j are the nodal values of the cross-section rotation (considered to be continuous at the element nodes). The lateral displacement is therefore where γ 0 is treated as an additional element degree of freedom. The detailed formulation of the finite-element equations for the beam element proceeds on the basis of standard approaches (Astley, 1992). The formulation employed for the soil elements adopts the displacement interpolation in equation (10). The local displacement and rotation in each soil element is where is a vector containing the element displacement degrees of freedom and where N′ i denotes a shape function derivative with respect to z. The internal force vector f and the tangent stiffness matrix k for the element are where p is a vector of force resultants p ¼ ð p mÞ T and the constitutive matrix, D, is The integrals in equation (13) are evaluated over each element using Gauss integration with four Gauss points per element.
Finite-element equations for the pile, the soil and separate lumped models at the pile base for H B and M B are assembled in the conventional manner; the resulting non-linear finite-element equations are solved by Newton-Raphson.
This model formulation requires two parameters to define the behaviour of the pile (EI, GAκ). The conventional thin-walled approximation is employed for the second moment of area, I, and cross-section area A. Additional parameters are required to define the soil reaction curves, as described below. In any practical application of the model, choices need to be made on the number of embedded pile elements (and associated soil elements) to employ in the analysis. A discussion of this issue is provided later in the paper in connection with the observed convergence characteristics of the model.

3D FINITE-ELEMENT CALIBRATION ANALYSES Specification of the calibration analyses
Calibration analyses have been conducted for the idealised problem geometry shown in Fig. 1. The piles (referred to as the 'calibration set') employed in the calibration include pile diameters, D, of 5 m, 7·5 m and 10 m, and values of L/D of 2 and 6. The load eccentricity, h, for an offshore monopile depends on whether the loading is dominated by wind or wave action; in the current calibration set the normalised eccentricity h/D is assumed to be between 5 (wave-dominated loading) and 15 (wind-dominated loading). The calibration set is specified in Table 1. 3D finite-element analyses of the performance of each of these calibration pile configurations were conducted using the finite-element software ICFEP (Potts & Zdravković, 1999, 2001 using the same constitutive model and procedures that had been validated against field test data; see Byrne et al. (2019b) and Zdravković et al. (2019b).
Most of the piles in the calibration set have a pile wall thickness, t ¼ D/110; this ratio provides wall thickness values that are regarded as being typical for realistic monopiles. Two additional calibration calculations with thicker pile walls (pile C3 with t ¼ D/80 and pile C7 with t ¼ D/60) were included in the calibration set; these cases were intended to explore the influence of pile wall thickness on the computed soil-pile interface reactions, and therefore to assess whether variations in D/t need to be included in the calibration process.

Initial ground conditions
The ground conditions for the calibration analyses are based on the Cowden test site that was employed for the PISA field tests Zdravković et al., 2019aZdravković et al., , 2019b. The ground conditions at this site consist of a deep layer (depth about 40 m) of overconsolidated glacial till. A range of historical data for the site (e.g. Powell & Butcher, 2003) were supplemented by additional triaxial and in situ data collected during the PISA project to characterise the soil conditions at the site . This site is affected by an under-drained pore water pressure profile , which is unrepresentative of offshore conditions. A separate ground model for a representative offshore glacial clay till site was therefore developed for the calibration analyses. This was based on the soil parameters determined for the Cowden test site , but with the vertical effective stress determined using a bulk unit weight, γ = 21·19 kN/m 3 , and hydrostatic pore water pressure from the ground surface. The (relatively small) adjustments to the effective stresses, as a consequence of adopting a hydrostatic pore pressure variation, imply small changes to the profiles of in situ undrained triaxial compression shear strength, s u , and small-strain shear modulus, G 0 . Profiles of s u , G 0 (G 0 ¼ 1100p′ is employed, based on the calibration in the paper by Zdravković et al. (2019b), where p′ is the mean effective stress), the coefficient of earth pressure at rest, K 0 , and the overconsolidation ratio employed in the ground model for the representative offshore glacial till site, which are consistent with the selected soil constitutive model (see Table 2 and Zdravković et al. (2019b)), are plotted in Fig. 4. Below 50 m the profile of s u continues linearly at the same gradient.

Finite-element meshes and boundary conditions
The finite-element meshes exploit symmetry; so that only one half of the problem needs to be discretised. An example mesh, for pile C4, is shown in Fig. 5. In this case 10 530 20-noded hexahedral displacement-based isoparametric solid elements are used to model the ground. The soil-pile interface is represented by 360 16-noded zero-thickness interface elements (Day & Potts, 1994), and the pile itself is modelled using 600 eight-noded shell elements (Schroeder et al., 2007). In the axial direction, each pile is discretised  with 30 rows of elements below the ground level and 20 rows of elements above. (The extension of the pile above ground level provides a convenient means of applying the desired moment M G to the monopile. The deformations induced in the above-ground extension, however, do not have any significance for the calibration study.) Similar meshes were employed for the other calibration piles. The nodes on the base of the mesh are fixed. Displacements normal to the vertical cylindrical boundary are prescribed to be zero, and appropriate symmetry conditions are applied to the nodes on the plane of symmetry. The lateral load at the pile top (i.e. at z ¼ Àh) is applied in a displacement controlled manner, such that the increments of displacement in the y-direction are applied uniformly around the pile perimeter. The lateral load, H, on the pile is determined from the computed reactions. The numerical calculations were conducted incrementally with an appropriate number of increments to ensure sufficient resolution for calibration purposes.

Constitutive models
An extended generalised version of the non-linear elastoplastic modified Cam Clay model was adopted for the soil (see Zdravković et al. (2019b)). The model utilises a Hvorslev surface on the dry side (Potts & Zdravković, 1999;Tsiampousi et al., 2013). The general Van Eekelen (1980) expression for strength variation in the deviatoric plane was employed. The non-linear degradation of shear modulus with strain and its dependence on stress level (Fig. 6) were simulated with the Taborda et al. (2016) small-strain model. The particular stiffness degradation employed in the analyses was determined from the site investigation data at the Cowden site . Undrained conditions were enforced in the analyses by prescribing a large value for the bulk modulus of the pore fluid (Naylor, 1974;Potts & Zdravković, 1999); the developed 1D model therefore applies only to undrained loading conditions.
The soil-pile interface was represented by an elasto-plastic Tresca model with a tensile capacity equal to the local hydrostatic pressure; this ensures that the appropriate hydrostatic pressure acts on the pile when gap formation occurs. As the pile is loaded, the horizontal total stress reduces on the active face; if it reaches the tensile capacity then a gap will start to form (from the surface downwards in the current analyses). Consequently, at the start of the analysis there is full bonding between the pile and the soil. However, as the lateral load increases, and a gap starts to open, then de-bonding progresses down the pile. In the current calculations the phreatic surface was located at ground level. The shear and normal stiffness employed in the interface model were both set to 1·0 Â 10 5 kN/m 3 . The shear strength in compression of the interface was set to the local value of undrained soil shear strength, approximating a fully rough interface. The monopile was modelled as an elastic material (representing steel), with Young's modulus E = 200 GPa and Poisson ratio ν = 0·3. The pile wall thickness (Table 1) is specified as an additional constitutive parameter. The 1D model presented below inherently accounts for this assumed constitutive behaviour. Variations of the 1D model could be developed that explore the sensitivity of the design model to . Shear modulus degradation curve employed for the representative offshore glacial clay till site. The generalised deviatoric strain, E d , is defined in the paper by Zdravković et al. (2019b) these assumptions, but that would require a separate suite of 3D finite-element calibration calculations.

Analysis procedures
The piles were 'wished in place' in the 3D calibration analyses, to avoid introducing uncertainties into the analyses around how installation effects can be realistically quantified. This analysis approach reflects that laterally loaded piles mobilise a relatively large volume of soil around the pile, which extends well beyond any interface zone disturbed during installation. This contrasts with axially loaded piles for which the performance of the pile depends strongly on the soil-pile interface conditions and installation effects may be significant. The appropriateness of the wished in place assumption is further reinforced by the excellent agreement between the numerical predictions and field measurements of the PISA test piles at both the Cowden  and Dunkirk  test sites.
The 3D finite-element calibration analyses were all taken to a ground-level pile displacement of at least v G ¼ D/10; this value of displacement was considered to correspond to the ultimate limit state of the pile. To ensure that the analyses had converged sufficiently tightly at a local level for model calibration purposes, an equilibrium check was performed on the numerical results. To conduct this check, the nodal forces acting on the pile and soil plug at the boundary with the surrounding soil were used to compute the overall force acting on the pile. These reactions were compared with the force implied by the applied loading. Numerical solutions were only accepted if these boundary checks were satisfactory (maximum overall force discrepancies of less than 5% were considered acceptable); in other cases the finite-element analyses were repeated using a tighter convergence tolerance.

NUMERICAL SOIL REACTION CURVES
The soil reaction curves were initially determined by extracting numerical data on the soil reactions from the 3D finite-element results at each loading increment. 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 of the vertical displacements (for rotation). The local distributed lateral load was determined by integrating the y-component of the computed normal and shear tractions at the external soil-pile interface. The local distributed moments were computed from where t z are the computed vertical tractions (defined positive upwards acting on the pile) and ϕ is defined in Fig. 7. The base horizontal force and moment were determined from the computed nodal forces at the base of the shell elements used to represent the pile and the computed nodal forces acting on the nodes at the base of the soil plug. The soil reaction data determined in this way were normalised using the dimensionless groups listed in Table 3. The local soil strength and stiffness data used in this normalisation process were calculated at the average depth of each element. This normalisation process is a key part of the procedure, as it allows the soil reaction curves to be computed within the 1D model for local values of s u and G 0 . The soil reaction data normalised in this way are referred to as 'numerical soil reaction curves'.
To check the robustness of the procedures used to compute the numerical soil reaction curves, and to assess the potential performance of a 1D model based on the PISA modelling approach, check calculations were conducted for all of the calibration piles using a form of the 1D model, referred to as 1D (numerical), in which the numerical soil reaction curves are incorporated directly. In these calculations, the numerical soil reaction curves at the depth location of each Gauss point in the 1D model were determined by interpolation; local values of soil strength and stiffness were then used to compute the local soil reactions.
Comparisons of the relationships between the lateral load, H, and the ground-level pile displacement, v G , computed using the 3D and 1D models are facilitated by an accuracy measure referred to as the 'accuracy metric', η, which is defined below, based on the nomenclature in Fig. 8.
An accuracy metric of η ¼ 1·0 indicates a perfect match between the 1D and 3D models. Two separate sets of η have been determined. To assess the accuracy of the 1D model for pile displacements up to an assumed limit state of v G ¼ D/10, values of the accuracy metric, η ult , were determined for 0 , v G , D/10. To obtain a separate indication of the small displacement performance of the 1D model (e.g. relevant for computing the natural frequencies of a wind turbine support structure) a small displacement accuracy metric, η sd , was determined for 0 , v G , D/10 000. The values of η determined from this process, for all of the calibration piles,  Fig. 7. Vertical shear tractions developed on an elemental length of pile at the soil-pile interface. The angle ϕ is defined relative to the loading direction as shown are 0·90 η ult 0·99 (mean = 0·95) and 0·95 η sd 0·99 (mean = 0·97); these results indicate a close match with the calibration data, providing confidence in the robustness of the modelling approach. Example comparisons between the H-v G performance determined from the 3D calibration data and the 1D (numerical) model are indicated in Fig. 9 (for pile C4) and Fig. 10 (for pile C1).
A further exercise has been conducted to investigate the significance of the individual soil reaction components on the performance of the 1D model. Two forms of the model have been considered; 'case P' (which incorporates only the distributed lateral load, p) and 'case P + H and M' (in which the base reactions, H B and M B , are included together with the distributed lateral load; the distributed moment is omitted). Results for case P for piles C4 and C1 are plotted in Figs 9 and 10. These data indicate that for pile C4 (which is relatively long, L/D ¼ 6) the 3D calibration data are matched closely by the case P model (Fig. 9). The loss of performance of the 1D model caused by omitting the distributed moment and the base force and moment terms, for this particular case, appears small. For pile C1 (which is relatively short, L/D ¼ 2), however, considerable loss of accuracy is apparent in the case P model (Fig. 10). Also plotted on Figs 9 and 10 are the results for case P + H and M. For the long pile (C4) there is no significant loss of fidelity when the distributed moment term is omitted from the model (case P + H and M). For the short pile (C1), however, the omission of the distributed moment does cause a discernible loss of accuracy. These observations support the underlying assumption of the PISA model in Fig. 2, that a complete set of four soil reaction components is required for a realistic 1D model for the behaviour of relatively short, large diameter monopiles, such as pile C1.
Although the 1D (numerical) model is successful in reproducing the 3D calibration data, it does not have any predictive capability. However, by developing and calibrating general analytical expressions for the soil reaction curves, it is possible to formulate a form of the 1D model, referred to as '1D (parametric)' that can be applied to any pile and loading geometry within the calibration space. The development of the 1D (parametric) model is described below.

Functional form
The 1D (parametric) model is based on the use of the following conic function to represent the soil reaction curves À nȳ y u Àx wherex is a normalised displacement (or rotation) variable andȳ is the corresponding normalised load (or moment) variable as listed in Table 3. This function requires the specification of four separate parameters: ultimate displacement,x u , ultimate load,ȳ u , initial stiffness, k, and the curvature parameter, n, where 0 n 1 andx u .ȳ u =k.
The general form of the function is plotted in Fig. 11(a).
For the extreme cases of n = 0 and n ¼ 1, the function reduces to a bilinear form, as illustrated in Fig. 11(b). The choice of this particular function is not essential for the methods described here. The function was chosen because it allows fitting of (a) initial stiffness, (b) an ultimate load, (c) a displacement at which ultimate load is reached and (d) a curvature parameter which determines the sharpness of the transition from high initial stiffness to softer non-linear response. A number of other functions have been proposed in the literature with similar characteristics, and in principle the fitting exercise described below could be carried out using an alternative kernel function to fit each curve. The normalised soil reactions can be determined explicitly from the normalised displacements bȳ where, Parameters to fit the soil reaction curves are determined using a two-stage process. A 'first-stage calibration' is conducted based on the numerical soil reaction curves determined from the 3D finite-element calibration analyses. These parameters are considered to vary with depth according to functions referred to as 'depth variation functions'. Then, a 'second-stage optimisation' process is conducted in which small changes are made to the depth variation parameters to improve the fit between the 1D model and the calibration data for the pile head performance. These processes are described below, with Table 4 reporting the depth variation parameters for both the first-stage calibration and after the second-stage optimisation.
First-stage calibration Distributed lateral load curves. As the lateral load reaction on the pile contributes the largest term to the overall pile response (see case P in Fig. 10), this is the most important load component to model accurately. Parameters for the normalised distributed lateral load response were determined from the numerical soil reactions as follows. The initial stiffness k p was selected by least-squares fitting the linear expressionp ¼ k pv to the numerical data in the small displacement region of 0 ,v , 10. (c) The ultimate response,p u , and curvature, n p , were determined by minimising the proportional least-square error between the numerical data and the parametric expressions.
Comparisons between the numerical and parametric soil reaction curves for pile C4 determined in this way are plotted in Fig. 12(a). The parameters k p , n p andp u are plotted as a function of normalised depth z/D in Fig. 13 for all of the piles in the calibration set. These data do not indicate any systematic influence of pile wall thickness (i.e. the results for piles C3 and C7 are consistent with the other data). A similar insensitivity to pile wall thickness is observed in the other soil reaction components. This has the desirable implication that variations in pile wall thickness may not be required for future calibration exercises.
A best-fit line, representing the depth variation function (shown as 'First stage' in Fig. 13) is determined for the initial stiffness, k p ( Fig. 13(a)), and curvature parameter, n p ( Fig. 13(b)), using weighted least-squares regression (where the weighting factors are proportional to the number of data points at each discrete depth). The normalisation adopted for the initial stiffness parameter k p typically leads to unrealistically large values when the shear modulus tends to zero near the ground surface; the fitting of the linear depth variation function for k p was therefore only conducted for z/D . 0·2.
Data onp u for all of the calibration piles are plotted in Fig. 13(c). The regions of poor fit around normalised depths of z/D = 1·4 and z/D = 4·4 occur near the pile rotation point (for values of L/D of 2 and 6, respectively), where the lateral displacement developed in the calibration analysis appears to be too small to mobilise the ultimate lateral capacity. Apart from these local anomalies, the data indicate a general tendency forp u to increase with depth, from a value of about 3 at the soil surface to about 10 at z/D ¼ 6. This tendency is consistent with the numerical solutions, based on a surface wedge mechanism, in the paper by Murff & Hamilton (1993). For long piles and relatively large depths, the soil is expected to fail in a flow-around mechanism, implyingp u ¼ 11Á94 for a rough pile (Randolph & Houlsby, 1984;Martin & Randolph, 2006), which is consistent with the current modelling approach. The data in Fig. 13(c) suggest that, in all the cases considered here, the piles are too short to mobilise this flow mechanism.
Following Murff & Hamilton (1993) the depth variation function for the ultimate resistancep u is selected as Distributed moment, m Ultimate rotation,ψ mu Given bym u =k m Given bym u =k m The soil reaction curve parameters 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 the parameters. 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.
where N 1 , N 2 and ξ are parameters determined from the numerical data by non-linear regression. The PISA derived variation, specified in Table 4, broadly follows that given by Murff & Hamilton (1993), despite different calculation methods being applied (see Byrne et al. (2017)). Further recent work exploring the variation of the parameters in the equation (22) expression for different clays has been completed by Zhang et al. (2016), with further recommendations given in the paper by Jeanjean et al. (2017). Distributed moment curves. Figure 12(b) shows an example set of numerical distributed moment soil reaction curves, for pile C4. These data indicate an initial peak followed by softening. A bilinear form of the parametric curve is adopted in this case by selecting n m = 0 and ensuring that ψ mu .m u =k m ; the resulting parametric curves are shown in Fig. 12(b). This approach implicitly removes the peak response. However, the peak response occurs over a relatively small range of normalised rotation, and at values of rotation that vary along the pile; as a consequence, this simple bilinear model has been found to be broadly satisfactory for modelling the overall behaviour of the pile. The ultimate momentm u is determined at each depth by taking the mean of the values that satisfym . 0Á9m final , wherem final is the value ofm at large displacements. A linear depth variation function is employed form u . The initial stiffness parameter k m is determined using a least-squares fit to the values in the regionm ,m u =10. Figure 14 shows variations of the initial stiffness parameter, k m , and the ultimate moment parameter, m u . In a similar way to the parameter k p , the initial stiffness parameter k m can lead to unrealistically large values near the ground surface. A linear depth variation function is adopted for k mindicated as 'First stage' in Fig. 14. This function is only fitted to the numerical data z/D . 0·2. A linear variation is fitted to m u , although there is more scatter in the abstracted reaction curve data for this parameter.
Base force and moment curves. The base horizontal force soil reaction curves extracted from all of the piles in the calibration set are shown in Fig. 15(a). The parameterv Hu is assumed to be independent of pile L/D. The parameters k H ; n H ;H Bu (determined at L/D ¼ 2 and L/D ¼ 6) are fit to a depth variation function that varies linearly with L/D.
Data on the base moment soil reaction curves, together with the best-fit parameterised curves for the first-stage calibration, are shown in Fig. 15(b).

Second-stage optimisation
The effectiveness of the first-stage calibration process is assessed using the accuracy metric specified in equation (11). The values of η ult determined using this process were in the range 0·90-0·98 and η sd values were in the range 0·77-0·92. Based on the sign conventions in Fig. 2(b), the numerical values of displacement in (a) and the corresponding values of force are negative; for convenience these data are plotted here in the first quadrant While this indicates a broadly satisfactory performance, an improved calibration can be obtained by a 'second-stage optimisation' in which the 1D model parameters are adjusted by small amounts to minimise the cost function where η ult,i and η sd,i are the accuracy metrics for pile C i (i ¼ 1:11). This second-stage optimisation process is conducted using standard procedures in Matlab. The final set of optimised parameters for the 1D model are specified in Table 4. The depth variation functions, following the second-stage optimisation, are plotted in Fig. 13 (indicated 'Second stage') for three of the distributed lateral load parameters. The value of the fourth parameter,v pu , was adjusted from its initial value of 200 to 241·4 by the second-stage optimisation process. The second-stage results are also plotted on Fig. 14 for the distributed moment parameters and Fig. 15 for the base shear and moment; it is apparent that a significant adjustment has been made through the optimisation process. Values of accuracy metrics for the individual piles, based on these parameters, are shown in Fig. 16. These data indicate that the small displacement metric for all calibration piles is in the range 0·95 η sd 0·99; the ultimate displacement metric is in the range 0·9 η ult 0·96. These accuracy ranges are similar to those achieved with the 1D (numerical) model. An alternative 'ratio metric', defined below, has been found to be useful for quantifying the accuracy of the 1D model where for a particular value of ground level displacement, H 1D is the value of lateral load computed using the 1D model and H 3D is the corresponding value from the 3D model (as defined in Fig. 8). The ratio metric quantifies the extent to which the 1D model overestimatesor underestimatesthe 3D calibration data at particular values of v G . For completeness, values of the ratio metric are also provided in Fig. 16. These data indicate no significant overall bias (for the ultimate load case: mean = 1·03, CoV = 5·0%; for the small displacement case: mean = 1·00, CoV = 3·3%), although it is noted that ρ 1 for the piles of larger diameter and ρ ! 1 for the piles of smaller diameter.
Convergence and parameter study for the 1D model Calculations employing the second-stage optimised parameters have been completed with different numbers of embedded elements for piles C1 and C4, to compute the applied force at v G ¼ D/10 and v G ¼ D/10 000. Computed force values (H ult at v G ¼ D/10 and H sd at v G ¼ D/10 000) are listed in Table 5, along with the percentage errors relative to the finest mesh used in each case. The results indicate that H ult is remarkably tolerant to the employment of a coarse mesh for both piles. For pile C1, for example, even with only two embedded elements, the error is less than 1%.
The pattern is slightly different for H sd . For C4, the error rises to 7% when the number of embedded elements is reduced to three. A mesh of 12 embedded elements in this case is needed to achieve an error of less than 1%. Pile C1, however, is less sensitive to mesh coarseness: four embedded elements in this case are sufficient to achieve an error of less than 1% on H sd . This behaviour is perhaps due to the fact that, for C4 (and to a lesser extent C1), the small-strain behaviour is strongly conditioned by bending deformations in the pile; this may imply the need for an increased number of elements to achieve an accurate result. All 1D model calculations presented in this paper, for comparison with the 3D finite-element calculations, adopt 20 embedded pile elements.
Additional calculations have been conducted for C1 and C4, employing the most refined meshes in each case, in which the value of the shear factor κ is artificially increased to a large number (1000), to suppress shear deformations in the pile. This process indicates the extent to which the shear deformations, incorporated in the 1D model, have an influence on performance. For pile C1, suppressing the shear deformations caused an increase of 0·02% and 0·6%, respectively, to H ult and H sd . For pile C4, suppressing the shear deformations caused an increase of 1·5% and 2·9%, respectively, to H ult and H sd . The greater sensitivity of C4 to shear deformations is considered to be due to the greater compliance of the pile. Overall, these results suggest that the influence of shear effects in the pile, for the soil conditions explored here, is relatively small compared with the uncertainties in other aspects of the problem.

PREDICTIVE CAPABILITY OF THE 1D (PARAMETRIC) MODEL
Once the 1D model has been calibrated, it can be used to determine the performance of a monopile for arbitrary geometric and loading conditions within the calibration space. To demonstrate the predictive capability of the 1D model, two design cases within the calibration space (D1 and D2) were analysed (using the 1D (parametric) model employing the second-stage parameters in Table 4). The geometric configurations for these test cases are specified in Table 6. These same pile configurations were separately analysed using 3D finite-element analysis.
The H-v G responses computed using the 1D and 3D models are shown in Fig. 17. The results indicate a close fit between the 1D and 3D finite-element data over the full range of applied loading (Fig. 17(a)). In particular, the small displacement response computed using the 1D model is seen to be almost identical to the response computed using the 3D finite-element model ( Fig. 17(b)). Figure 18 shows a comparison between the below-ground bending moments for the two design cases, computed using the 1D and 3D models. These data also indicate a close match.
The close fit that is obtained between the 1D and 3D finite-element results supports the use of the calibrated 1D model to conduct pile design calculations within the calibration space.

DISCUSSION
The PISA 1D design model provides a rapid means of computing the performance of a monopile foundation within a pre-determined calibration space. The 1D model is calibrated for particular soil conditions and for a pre-defined range of design parameters based on a set of 3D finite-element analyses. The 1D model provides a rational means of interpolating between the 3D finite-element calibration data, employing a Winkler model that is constrained (by way of the 1D finite-element formulation) to satisfy the requirements of equilibrium and compatibility, and the assumptions of Timoshenko beam theory for the pile. Whereas the analysis of a monopile using 3D   finite-element analysis invariably requires substantial computational resources, the 1D model can be computed using a standard desktop computer with run times of the order of seconds. Such reduced order, or 'surrogate', modelling is well established in other areas of engineering and holds promise for application in geotechnical engineering, as demonstrated here for monopile foundations. Clearly the fidelity that can be achieved with the PISA design model is limited by the quality of the 3D finite-element analyses employed in the calibration process.
The PISA design model, in the form presented here, employs a four-parameter conic function to represent soil reaction curves. In the first-stage calibration, the function is calibrated to provide a direct approximation of the numerical soil reaction curves, as abstracted from the 3D finite-element analyses. The fidelity of the 1D model (i.e. the extent to which it mimics the 3D finite-element model) is limited by the following factors (a) imperfect fitting of the numerical soil reaction curves by the conic function (b) the use of unique, simple functions in the 1D model to represent the depth variation of the conic function parameters, whereas the values of these parameters determined from the 3D calibration analyses exhibit considerable scatter and often complex variations with depth (c) the Winkler approach employed in the 1D model necessarily neglects various soil-pile interaction mechanisms (e.g. spatial coupling within the soil) that develop in the 3D calibration analyses.
The 1D model employing the first-stage calibration data was found in many cases to provide a close representation of the pile head response as determined by the 3D finite-element calibration analyses. The fidelity of the model, however, is fundamentally limited by the three factors outlined above. The second-stage optimisation process therefore improves the fidelity of the pile head performance computed using the 1D model, by adjustment of the depth variation parameters. The parameters determined from the second-stage optimisation do not necessarily provide an improved representation of the local soil-pile behaviour, and in some pile locations this may be less good. Instead, the second-stage optimisation is regarded as an expedient to compensate for the limitations of the Winkler framework by limited adjustment of the model parameters.
The data in Fig. 16 indicate that the 1D model employing the second-stage optimised parameters provides high-fidelity representations of the 3D finite-element results. The second-stage optimised form of the PISA design model is considered to be appropriate for practical design predictions of the pile head performance and the bending moments in the embedded monopile. The model should be employed in its entiretythat is, forms of the model in which selected soil reaction curve components are excluded from the analysis are not recommended. Calibration data obtained from the first-stage calibration process provide an alternative option for the design model. These data will tend to provide a closer representation of the local soil-pile behaviour as computed using 3D finite-element analysis, but the fidelity of the computed pile head behaviour is typically lower, especially for the small displacement response, than if the second-stage parameters are employed. For all applications of the PISA design model it is recommended that final designs, or at least a selection of them, are confirmed by independent 3D finite-element analysis.
It is important to note that the 1D model calibration parameters presented in Table 4 relate to the particular profiles of undrained shear strength and shear modulus employed in the calibration process. The limitations of the Winkler modelling approach mean that these parameters do not provide a general representation for arbitrary strength and stiffness profiles. Although the model is unlikely to be sensitive to small variations in these profiles (although this has not been systematically investigated), depth profiles of strength and/or stiffness that differ significantly from those considered here will require separate calibration. Similarly, soils formed from other types of clay geological units, or cases where the soil-pile interface behaviour is considered to differ from the assumptions employed in the current analysis, may also require a separate calibration process.
The PISA model can be used in two separate modes. For initial design studies, previously determined calibration data (such as those presented in Table 4) could be used; this approach is referred to as the 'rule-based method'. Although these parameters relate to a single soil type (glacial till) and for a single profile of soil strength and stiffness, it is expected that a database of appropriate parameters for different soil types will be assembled in due course. Initial work in this direction is described in the paper by Byrne et al. (2019a). For detailed design calculations, it is anticipated that sitespecific model parameters will be determined by way of a bespoke 3D finite-element calibration process; this approach is referred to as a 'numerical-based method'. The numericalbased approach has the capability of providing high-fidelity representations of the underlying 3D calibration models, as demonstrated in the current work. The reliability of the model is principally limited by the quality of the 3D model and the in situ and/or laboratory test data that are used to calibrate it. Further details on how these analysis procedures might be embedded within a general design framework are given in the papers by Byrne et al. (2017Byrne et al. ( , 2019a. Although developed for offshore wind turbine monopile foundations, the PISA modelling concept is capable of being generalised to other foundation design and soil-structure interaction problems. For example, it could be employed in developing site-specific p-y relationships for more flexible (large L/D) jacket piles, or to develop six-degrees-of-freedom models for suction caisson foundations. In either case, simpler, computationally efficient calculation methods could be calibrated against advanced 3D computations that capture non-linear effects through suitable constitutive models and soil characterisation. The model calibration processes could focus specifically on aspects of behaviour that are most relevant for design, such as initial stiffness, ultimate capacity, or the foundation response for both small and large displacements (as considered in this paper).

CONCLUSIONS
The paper demonstrates an approach in which a simplified computational model of an embedded monopile foundation (based on a 1D finite-element approach) is calibrated using data from more complex 3D finite-element analysis. For design applications the calibration analyses are required to span the likely range of pile embedded lengths, diameters and load eccentricities. For a representative offshore glacial till site, the 1D model is shown in the current paper to provide a close representation of the 3D calibration data; this confirms the validity of the modelling assumptions and calibration procedures.
The predictive capabilities of the model are demonstrated, by way of a set of design calculations based on pile geometries within the calibration space but not included in the calibration data. These example calculations indicate a close match between the 1D model and an equivalent 3D finite-element analysis.