Offshore wind turbines in shallow coastal waters are typically supported on monopile foundations. Although threedimensional (3D) finiteelement methods are available for the design of monopiles in this context, much of the routine design work is currently conducted using simplified onedimensional (1D) models based on the p–y method. The p–y method was originally developed for the relatively large embedded lengthtodiameter 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 finiteelement 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 finiteelement calibration analysis, but at a fraction of the computational cost. Moreover, within the calibration space, the 1D model is capable of delivering highfidelity 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.
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 onedimensional (1D) analysis framework – known as the ‘p–y’ method – in which the monopile is modelled as an embedded beam and the lateral soil response is represented by ‘p–y’ curves, which are nonlinear 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 DNVGL (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. 
The development of the design model presented here was the purpose of a recently completed joint industry study – known as PISA (Byrne et al., 2015, 2017, 2019a; Zdravković et al., 2015; Burd et al., 2017) – comprising field testing and ground characterisation (Burd et al., 2019; Byrne et al., 2019b; McAdam et al., 2019; Zdravković et al., 2019a), 3D finiteelement modelling (Taborda et al., 2019; 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 (Byrne et al., 2019b; Zdravković et al., 2019b). The other was at Dunkirk in northern France, where the soil is a dense marine sand (McAdam et al., 2019; 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).
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 crosssection 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 finiteelement 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 soil–structure 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 nonlinear pile behaviour would not be straightforward.
A fourparameter function is used as the basis for each of the soil reaction curves. The curves are calibrated directly from a set of bespoke 3D finiteelement 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.
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 Galerkin form of equation (8) is obtained by discretising the pile using a line mesh of twonoded, fivedegreesoffreedom Timoshenko beam elements, based on the formulation in the book by Astley (1992). The soil is represented by a line mesh of twonoded, fivedegreesoffreedom 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 finiteelement formulation is provided below.
The integrals in equation (13) are evaluated over each element using Gauss integration with four Gauss points per element.
Finiteelement 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 nonlinear finiteelement equations are solved by Newton–Raphson.
This model formulation requires two parameters to define the behaviour of the pile (EI, GAκ). The conventional thinwalled approximation is employed for the second moment of area, I, and crosssection 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.
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 (wavedominated loading) and 15 (winddominated loading). The calibration set is specified in Table 1. 3D finiteelement analyses of the performance of each of these calibration pile configurations were conducted using the finiteelement 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).

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

C1  10  50  5  20  2  91  110 
C2  10  150  15  20  2  91  110 
C3  10  50  5  20  2  125  80 
C4  10  50  5  60  6  91  110 
C5  10  150  15  60  6  91  110 
C6  5  25  5  10  2  45  110 
C7  5  25  5  10  2  83  60 
C8  5  25  5  30  6  45  110 
C9  5  75  15  30  6  45  110 
C10  7·5  37·5  5  15  2  68  110 
C11  7·5  37·5  5  45  6  68  110 
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.
The ground conditions for the calibration analyses are based on the Cowden test site that was employed for the PISA field tests (Byrne et al., 2017; Zdravković et al., 2019a, 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 (Zdravković et al., 2019a). This site is affected by an underdrained pore water pressure profile (Zdravković et al., 2015, 2019a), 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 (Zdravković et al., 2019a, 2019b), but with the vertical effective stress determined using a bulk unit weight, γ = 21·19 kN/m^{3} (Zdravković et al., 2019a), 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 smallstrain 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.

Component  Parameters 

Strength (Van Eekelen, 1980), equation (3)  X = 0·548, Y = 0·698, Z = 0·100 
Nonlinear Hvorslev surface – shape (Tsiampousi et al., 2013), equation (1)  α = 0·25, n = 0·40 
Nonlinear Hvorslev surface – plastic potential (Tsiampousi et al., 2013), equation (2)  β = 0·20, m = 1·00 
Virgin consolidation line  ν_{1} = 2·20, λ = 0·115 
Nonlinear elasticity – swelling behaviour  κ = 0·021 
Nonlinear elasticity – smallstrain shear modulus (Taborda et al., 2016), equation (4)  G_{0}^{*} = 110 MPa, p′_{ref} = 100·0 kPa 
Nonlinear elasticity – shear stiffness degradation (Taborda et al., 2016), equation (4)  a = 9·78 × 10^{−5}, b = 0·987, R_{G,min} = 0·05 
The finiteelement 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 20noded hexahedral displacementbased isoparametric solid elements are used to model the ground. The soil–pile interface is represented by 360 16noded zerothickness interface elements (Day & Potts, 1994), and the pile itself is modelled using 600 eightnoded 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 aboveground 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 ydirection 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.
An extended generalised version of the nonlinear 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 nonlinear degradation of shear modulus with strain and its dependence on stress level (Fig. 6) were simulated with the Taborda et al. (2016) smallstrain model. The particular stiffness degradation employed in the analyses was determined from the site investigation data at the Cowden site (Zdravković et al., 2019a, 2019b). 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 elastoplastic 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 debonding 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 these assumptions, but that would require a separate suite of 3D finiteelement calibration calculations.
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 (Zdravković et al., 2019b) and Dunkirk (Taborda et al., 2019) test sites.
The 3D finiteelement calibration analyses were all taken to a groundlevel 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 finiteelement analyses were repeated using a tighter convergence tolerance.
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’.

Normalised variable  Nondimensional form 

Distributed load,  
Lateral displacement,  
Distributed moment,  
Pile crosssection rotation,  
Base horizontal force,  
Base moment, 
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.
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, 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.
Parameters to fit the soil reaction curves are determined using a twostage process. A ‘firststage calibration’ is conducted based on the numerical soil reaction curves determined from the 3D finiteelement calibration analyses. These parameters are considered to vary with depth according to functions referred to as ‘depth variation functions’. Then, a ‘secondstage 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 firststage calibration and after the secondstage optimisation.

Soil reaction component  Parameter  Depth variation functions Firststage calibration  Depth variation functions Secondstage calibration 

Distributed lateral load, p  Ultimate displacement,  200·0  241·4 
Initial stiffness, k_{p}  
Curvature, n_{p}  
Ultimate reaction,  10·21 − 7·215e^{[−0·3332(z/D)]}  10·70 − 7·101e^{[−0·3085(z/D)]}  
Distributed moment, m  Ultimate rotation,  Given by  Given by 
Initial stiffness, k_{m}  
Curvature, n_{m}  0·0  0·0  
Ultimate moment,  
Base horizontal force, H_{B}  Ultimate displacement,  300  235·7 
Initial stiffness, k_{H}  
Curvature, n_{H}  
Ultimate reaction,  
Base moment, M_{B}  Ultimate rotation,  200  173·1 
Initial stiffness, k_{M}  
Curvature, n_{M}  
Ultimate reaction, 
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.
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.
(a)  A single value of ultimate displacement,  
(b)  The initial stiffness k_{p} was selected by leastsquares fitting the linear expression  
(c)  The ultimate response, 
A bestfit 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 leastsquares 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.
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_{m} – indicated 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.
Data on the base moment soil reaction curves, together with the bestfit parameterised curves for the firststage calibration, are shown in Fig. 15(b).
Calculations employing the secondstage 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%.

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

C1  
200  0·1  12·9636  0·0000  0·7395  0·0000 
40  0·5  12·9639  0·0023  0·7396  0·0135 
20  1  12·9644  0·0062  0·7396  0·0135 
10  2  12·9658  0·0170  0·7399  0·0541 
4  5  12·9651  0·0116  0·7420  0·3381 
2  10  12·8759  −0·6765  0·7492  1·312 
C4  
120  0·5  104·0325  0·0000  0·9105  0·000 
60  1  104·0331  0·0006  0·9106  0·011 
24  2·5  104·037  0·0043  0·9111  0·066 
12  5  104·0244  −0·0078  0·9132  0·297 
6  10  104·2063  0·1671  0·9225  1·318 
4  15  103·2956  −0·7083  0·9385  3·075 
3  20  104·2368  0·1964  0·9743  7·007 
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 smallstrain 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 finiteelement 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.
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 secondstage 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 finiteelement analysis.

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

D1  7·5  37·5  5  22·5  3  68  110 
D2  8·75  87·5  10  35  4  91  97 
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 finiteelement 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 finiteelement model (Fig. 17(b)).
Figure 18 shows a comparison between the belowground 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 finiteelement results supports the use of the calibrated 1D model to conduct pile design calculations within the calibration space.
The PISA 1D design model provides a rapid means of computing the performance of a monopile foundation within a predetermined calibration space. The 1D model is calibrated for particular soil conditions and for a predefined range of design parameters based on a set of 3D finiteelement analyses. The 1D model provides a rational means of interpolating between the 3D finiteelement calibration data, employing a Winkler model that is constrained (by way of the 1D finiteelement 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 finiteelement 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 finiteelement analyses employed in the calibration process.
The PISA design model, in the form presented here, employs a fourparameter conic function to represent soil reaction curves. In the firststage calibration, the function is calibrated to provide a direct approximation of the numerical soil reaction curves, as abstracted from the 3D finiteelement analyses. The fidelity of the 1D model (i.e. the extent to which it mimics the 3D finiteelement 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 secondstage 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 entirety – that 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 firststage 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 finiteelement analysis, but the fidelity of the computed pile head behaviour is typically lower, especially for the small displacement response, than if the secondstage 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 finiteelement 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 ‘rulebased 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 finiteelement calibration process; this approach is referred to as a ‘numericalbased method’. The numericalbased approach has the capability of providing highfidelity 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. (2017, 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 sitespecific p–y relationships for more flexible (large L/D) jacket piles, or to develop sixdegreesoffreedom models for suction caisson foundations. In either case, simpler, computationally efficient calculation methods could be calibrated against advanced 3D computations that capture nonlinear 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).
The paper demonstrates an approach in which a simplified computational model of an embedded monopile foundation (based on a 1D finiteelement approach) is calibrated using data from more complex 3D finiteelement 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 finiteelement analysis.
ACKNOWLEDGEMENTS
The PISA Phase 1 project was funded by the UK Department for Energy and Climate Change (DECC) and the PISA industry partners under the umbrella of the offshore wind accelerator (OWA) programme, which was designed and is led by the Carbon Trust. The authors acknowledge the provision of financial and technical support by the following project partners: Ørsted Wind Power (formerly DONG Energy), Alstom Wind, E.ON, EDF, Equinor (formerly Statoil), innogy, SPR, Statkraft, SSE, Vattenfall and Van Oord. The authors are grateful for the contributions of Dr Christelle Abadie, who contributed to the stage 2 calibration as part of the PISA Phase 2 project.
D  pile diameter 
G_{0}  smallstrain soil shear modulus 
H  lateral load applied to pile 
H_{1D}  lateral load applied to pile, computed with the 1D model 
H_{3D}  lateral load applied to pile, computed with the 3D model 
H_{B}  horizontal force at pile base 
H_{G}  lateral load applied to pile at ground level 
H_{sd}  lateral load applied to pile at a groundlevel displacement of v_{G} = D/10 000 
H_{ult}  lateral load applied to pile at a groundlevel displacement of v_{G} = D/10 
h  load eccentricity 
k  initial stiffness of parametric soil reaction curve 
L  pile embedded length 
M_{B}  moment at pile base 
M_{G}  moment applied to pile at ground level 
m  distributed moment acting on monopile 
n  curvature parameter for parametric soil reaction curve 
p  distributed lateral load acting on pile 
s_{u}  undrained shear strength of soil 
t  pile wall thickness 
v  lateral pile displacement 
v_{B}  lateral pile displacement at pile base 
v_{G}  groundlevel lateral pile displacement 
ultimate displacement for parametric soil reaction curve  
ultimate load for parametric soil reaction curve  
z  depth coordinate along the pile 
η  accuracy metric 
η_{sd}  small displacement accuracy metric (groundlevel pile displacements up to D/10 000) 
η_{ult}  ultimate displacement accuracy metric (groundlevel pile displacements up to D/10) 
ρ  ratio metric 
ψ  rotation of the pile crosssection 
ψ_{B}  rotation of the pile crosssection at the pile base 
Discussion on this paper closes on 1 March 2021, for further details see p. ii.