Application of the PISA design model to monopiles embedded in layered soils

Design procedures for monopile foundations for offshore wind turbine applications typically employ simplified models to facilitate the development of practical designs. A widely used simplified procedure, known as the p–y method, employs a beam model for the monopile and a Winkler representation of the pile–soil interaction. Non-linear functions (p–y curves) are specified to relate the pile displacement, y, to the local distributed lateral load, p, acting on the embedded pile. Functional forms of the p–y curves, and calibration parameters for sand and clay soil types, are specified in design guidance documents (e.g. API, 2010; DNV, 2016).

Although the p–y method is widely used for offshore monopile design, there is awareness that standard forms of the method may not provide realistic representations of behaviour for the relatively large-diameter monopiles that are now employed for offshore wind turbine applications. These concerns are informed by observations that the fundamental natural frequencies of wind turbine support structures – measured by way of supervisory control and data acquisition (Scada) instrumentation – are often systematically higher than those implied by the analysis models employed in the design process. A summary of issues relating to the limitations of the p–y method for monopile design is given in Doherty & Gavin (2011).

In response to perceived shortcomings of standard forms of the p–y method for monopile design applications, a new design approach, termed the ‘PISA design model’ has recently been developed (Burd et al., 2020a; Byrne et al., 2020a). This model, also referred to below as the ‘one-dimensional (1D) model’, was an outcome of a project – known as PISA – that incorporated ground characterisation (Zdravković et al., 2020a), field testing (Burd et al., 2020b; Byrne et al., 2020b; McAdam et al., 2020), three-dimensional (3D) finite-element analysis (Taborda et al., 2020; Zdravković et al., 2020b) and design model development (Burd et al., 2020a; Byrne et al., 2020a). The PISA design model employs the same basic Winkler modelling concept that forms the basis of the p–y method, but extensions and enhancements are incorporated to improve the model's performance, notably by incorporating additional soil reaction components (distributed moment, base horizontal force and base moment). The design model is intended to be calibrated using bespoke 3D finite-element calibration analyses for specific site conditions, rather than relying on general purpose calibration charts and correlations. This calibration procedure has been shown (Burd et al., 2020a; Byrne et al., 2020a) to be effective in training the design model to provide high-fidelity representations of the performance of monopiles embedded in homogeneous soils – as computed using 3D finite-element analysis – for design conditions within the parameter space employed in the calibration process.

In previous work (AWG, 2018; Burd et al., 2020a; Byrne et al., 2020a), PISA design model calibrations were developed for representative offshore sites comprising a single, homogeneous soil type. In reality, however, offshore wind farm sites typically consist of interbedded layers of soil with different geological origins. The current paper tests the application of the PISA design model to these layered soil conditions by exploring the hypothesis ‘soil reaction curves calibrated using homogeneous soil profiles can be employed, directly, to conduct 1D analyses of monopiles embedded in a layered soil’. This hypothesis follows the early work of Reese et al. (1981) for more flexible, slender piles appropriate to oil and gas applications, although there is limited supporting evidence in the literature. Subsequent studies by Georgiadis (1983) identified simplified methods to adapt standard p–y calculations to account for layering; these have since been incorporated as an option in commercial software (e.g. Isenhower et al., 2019). More recent work by Yang & Jeremic (2005) explored pile response in two- and three-layer systems using finite-element analyses, identifying that layering generally has an effect, but without providing specific design recommendations. There are, however, surprisingly few comprehensive studies into the general effect of soil layering on the lateral response of piles that could support the development of new simplified design procedures for large-diameter rigid monopiles for wind turbines.

The current paper describes a study on the response of an embedded monopile at offshore sites where the soil consists of interbedded layers with different characteristics. A range of idealised cases are considered, comprising combinations of homogeneous reference soils for which PISA design model calibrations have previously been obtained. The reference soils that have been employed (referred to here as ‘reference homogeneous soil models’) are: (a) a stiff overconsolidated clay till known as ‘Cowden till’ (Byrne et al., 2020a); (b) a marine sand known as ‘general Dunkirk sand model’ (GDSM) with relative density D_R in the range 45% ≤ D_R ≤ 90% (Burd et al., 2020a); and (c) a soft clay known as ‘Bothkennar clay’ (AWG, 2018). A set of 3D finite-element analyses of the performance of monopiles embedded in the idealised layered system cases has been conducted. Results from these analyses illustrate the various influences that soil layering has on the monopile performance. Separately, the performance of the embedded monopiles is computed using the PISA design model. Comparison of the PISA analyses with the 3D finite-element results facilitates an assessment of the extent to which the PISA design model is a satisfactory surrogate for detailed 3D finite-element analyses, for monopiles embedded in layered soils.

The analyses presented in this paper relate to the idealised configuration shown in Fig. 1 in which a hollow circular steel monopile of external diameter D and uniform wall thickness t is embedded to depth L at a site consisting of an arbitrary arrangement of horizontal soil layers. The monopile is loaded by a lateral load H applied at a distance h above the seabed (referred to as ‘ground level’ in the current work). The analyses are for monotonic loading only, although it is considered that small displacement secant stiffness data for monopile foundations, implied by the model, are applicable to the assessment of the dynamic performance of complete foundation–support structure systems. Design applications of the PISA design model, and other essential detail relating to calibration procedures and performance metrics, are discussed in detail by Byrne et al. (2020a) and Burd et al. (2020a). Pile installation effects are not considered, with the pile being modelled as ‘wished in place’.

Technical advice on the choice of layered soil configurations for consideration in the study was provided by the industry partners (listed in the Acknowledgements).

The PISA design model representation of an embedded monopile is illustrated in Fig. 2(a). The formulation and implementation of the model is described in Byrne et al. (2020a) and Burd et al. (2020a). The monopile is modelled as a line of two-noded Timoshenko beam finite elements. Nodes are generated at each soil layer boundary; within each layer the mesh is defined to a desired level of refinement. Four components of soil reaction act on the pile (distributed lateral load, p; distributed moment, m; base horizontal force, H_B; and base moment, M_B). Relationships between these reactions and the local lateral pile displacements, v, and cross-section rotations, ψ, are represented by functions known as ‘soil reaction curves’. Parameters are selected to calibrate the soil reaction curves for the particular soil type within each layer. The soil response is implemented in the model by way of consistent two-noded soil finite elements. Numerical solutions are obtained by employing four Gauss points within each pile and soil element and computing equilibrium solutions by Newton–Raphson iteration.

A conic function (see Appendix 1, Fig. 18) is adopted to represent each of the four soil reaction curves in the model in normalised form, $y¯=y¯(x¯)$⁠, employing the non-dimensional forms listed in Appendix 1, Table 11. Separate normalisation frameworks are adopted for clay and sand soils. The normalised soil reaction curves have four calibration parameters; k, $y¯u$⁠, $x¯u$⁠, n; these parameters in general vary with depth; functions describing these spatial variations – referred to as ‘depth variation functions’ – are specified in Appendix 1, Table 12 for the reference homogeneous soil models employed in the current study.

The PISA design model is calibrated by conducting a set of 3D finite-element analyses for monopiles embedded in a single soil type. Previous calibrations for the reference homogeneous soil models employed the parameter space 2 ≤ L/D ≤ 6; 5 ≤ h/D ≤ 15 (AWG, 2018; Burd et al., 2020a; Byrne et al., 2020a). It is demonstrated in Byrne et al. (2020a), Burd et al. (2020a) and AWG (2018) that, for monopiles with characteristics that lie within the calibration space, the relationships between the applied lateral load H and the ground-level lateral pile displacement v_G computed using 3D finite-element analysis can be represented – to a high fidelity – by the PISA design model.

For layered soil applications of the PISA design model, soil reaction curves are defined at points along the embedded pile as follows. First, the local soil type (clay or sand) is identified. This determines which of the frameworks in Table 11 is relevant. The local normalised soil reaction curve parameters are then determined, employing appropriate depth variation functions (Table 12). Dimensioned soil reaction curves are computed within the model implementation on the basis of the dimensionless forms in Table 11, employing the local values of small-strain shear modulus G₀, triaxial compression undrained shear strength s_u and initial vertical effective stress $σ′vi$ (depending on the framework) and the mathematical definition of the conic function (Byrne et al., 2020a). In other words, the original method (calibrated for homogeneous soils) is simply applied locally at each level in the pile, even for a non-homogeneous soil profile.

Three-dimensional finite-element analyses were conducted using the Imperial College finite-element program, ICFEP (Potts & Zdravkovic, 1999). Clay layers (Cowden till, Bothkennar clay) are modelled using the extended generalised non-linear elasto-plastic modified Cam Clay model that was previously employed to analyse the PISA field tests (Zdravković et al., 2020b). The use of this constitutive model to develop PISA design model calibrations for Cowden till and Bothkennar clay is described in Byrne et al. (2020a) and AWG (2018), respectively. Sand layers are represented with the bounding surface constitutive model that was employed to analyse the PISA test piles at the Dunkirk site (Taborda et al., 2020) as well as forming the basis of the calibration process for the general Dunkirk sand model (Burd et al., 2020a). Clay layers are treated as undrained; sand layers are drained.

An example finite-element mesh is shown in Fig. 2(b). The models exploit symmetry and so only one-half of the problem is discretised. The soil is represented by 20-noded hexahedral isoparametric solid elements. The soil–pile interface is modelled using 16-noded zero-thickness interface elements (Day & Potts, 1994), and the pile is modelled with eight-noded shell elements (Schroeder et al., 2007). The pile is extended above ground level to provide a convenient means of imposing the applied lateral load H (by way of prescribed displacements). Constitutive parameters and analysis procedures are identical to those described by Byrne et al. (2020a) (for Cowden till), AWG (2018) (for Bothkennar clay) and Burd et al. (2020a) (for Dunkirk sand).

The performance of the monopile is characterised by the relationship between the applied lateral load, H, and the computed ground-level pile displacement v_G (this is referred to as the ‘H–v_G response’). For the 3D analyses, v_G is determined as the average of the horizontal displacement on the front and back faces of the pile at ground level. Two reference displacements are considered. Computed values of lateral load, H_sd, at an assumed ‘small displacement’ v_G = D/10 000 are employed to quantify the secant stiffness of the monopile for low-amplitude loading. Separate data on the lateral load, H_ult, computed at an assumed ‘ultimate displacement’ v_G = D/10 are used to characterise the ultimate limit state of the pile.

Profiles of small-strain shear modulus G₀ and (for the clay soils) undrained shear strength s_u for the reference homogeneous soil models are shown in Fig. 3. For the general Dunkirk sand model, data are shown for values of relative density at the two extremes of the calibrated range, D_R = 45% and D_R = 90%. Data on submerged unit weights and the values of K₀ adopted in these models are specified in Appendix 2.

Monopile configurations employed in the study are specified in Table 1. These configurations (and the pile reference codes) are consistent with previous work (e.g. Burd et al., 2020a; Byrne et al., 2020a). Most of the analyses employ a monopile with embedment ratio L/D = 4. Two variants are considered; pile D2 (wall thickness t = 91 mm) and pile D2t (t = 150 mm). The thicker walled variant was adopted in some cases to avoid unrealistically high bending stresses developing in the pile wall. It should be noted that the original calibration of the PISA design model employed monopiles with embedment ratio of either 2 or 6; the L/D = 4 piles adopted for the majority of the current analyses are therefore regarded as independent of the calibration pile configurations. A more limited exploratory study of layered soils was conducted on piles with embedment ratio L/D = 2. This configuration was motivated by current trends, reported by the industry partners, for the specification of pile embedment ratios in the region of this value for the next generation of offshore wind-farm designs. Again, different wall thicknesses are adopted (C1 and C1t) for these exploratory analyses.

The layered soil configurations considered are illustrated in Fig. 4; they represent idealised configurations of increasing complexity and realism. To specify these configurations the term ‘matrix’ is used to describe the principal soil type and the term ‘layer’ describes any embedded or base layers that are present. Case A comprises two soil types with a boundary at the pile mid-depth. Case B employs a single embedded layer within an otherwise homogeneous matrix. In case C, a layer is placed at – or just below – the pile base; this case was used to investigate the potential benefits of extending a monopile into a relatively stiff soil layer at depth. The case D configurations comprise selected combinations of cases B and C. Case E is a set of more realistic multi-layer systems.

Cases A to E are analysed with a monopile of embedment ratio L/D = 4 (piles D2 or D2t). In a separate set of analyses – case F – selected cases are re-analysed for a monopile with L/D = 2 (piles C1 or C1t).

These analysis cases have all been previously introduced in Byrne et al. (2019). Although the analysis cases are the same, the analysis naming conventions have been modified for the purpose of the current paper.

The initial stress, strength and stiffness profiles for the layered soil models are specified in the 1D and 3D analyses as follows.

Initial pore pressures are hydrostatic. Profiles of initial vertical effective stress $σ′vi$ are determined from the submerged unit soil weights for the reference soils (Table 13). Horizontal stresses $σ′hi$ (required for the 3D analyses) are evaluated using the values of K₀ employed in the reference homogeneous soil models (Table 13 and Fig. 19).

Procedures are needed to determine the profile of small-strain shear modulus G₀ and, for clay layers, the profile of undrained shear strength, s_u. Although G₀ and s_u are specified in the representative models as functions of depth (Fig. 3), G₀ and s_u are in practice dependent on the magnitude of the local stresses. In the current layered soil models, therefore, G₀ and s_u are considered to be functions of local mean effective stress, G₀(p′), s_u(p′). To implement these functional relationships in the analyses, pairs of numerical data [G₀, p′] and [s_u, p′] are determined from the reference homogeneous soil models at selected depths, z. These data pairs are used to specify piecewise linear functions to represent G₀(p′), s_u(p′). Piecewise linear functions G₀(z), s_u(z) – for incorporation in the layered soil analysis models – are then determined in a straightforward way for layered soil cases by employing the mean effective stress profile p′ = p′(z) that corresponds to the previously determined profile of vertical and horizontal effective stresses.

To provide baseline data for comparisons with the layered soil results, 3D finite-element analyses were conducted for each of the three reference homogeneous soil models employed in the study. For Dunkirk sand, two cases are considered, D_R = 45% and D_R = 90%. Calculations were conducted employing the piles C1, D2, D2t. Pile C1t was not included in these baseline calculations.

Computed data for piles C1 and D2 are plotted in Fig. 5 for both the small and ultimate displacement ranges. These piles exhibit a considerable range of initial stiffness and ultimate capacity; piles embedded in Cowden till have the greatest initial secant stiffness, whereas piles embedded in Dunkirk sand, D_R = 90%, have the greatest ultimate capacity. Piles embedded in Bothkennar clay lie at the lower end of the stiffness and capacity ranges in all cases. Numerical data on reference lateral loads for these baseline analyses are listed in Table 2.

The influence of the embedded layer and/or base layer on the H–v_G response computed using the 3D analyses is quantified at the reference ground-level pile displacements (v_G = D/10 000 and D/10) in terms of a ‘mobilisation index’ χ and an ‘influence factor’ IF defined by

where H is the computed lateral load for the layered analysis at a particular reference ground-level pile displacement; H_matrix and H_layer are the computed lateral loads for the pile when embedded in reference homogeneous soils corresponding to the matrix and layer, respectively, at the same reference pile displacement (data on H_matrix and H_layer correspond to the baseline results in Table 2). Values of these metrics are determined separately at small (χ_sd and IF_sd at v_G = D/10 000) and ultimate (χ_ult and IF_ult at v_G = D/10) reference pile displacements. A mobilisation index χ = 0 signifies that the layer has no effect on the monopile performance; conversely χ = 1 signifies that the monopile behaviour is determined entirely by the layer. In actual cases the mobilisation index is expected to lie in the range 0 ≤ χ ≤ 1 (although exceptions can occur). The influence factor, IF, provides a separate measure on the proportional change in computed lateral load associated with the presence of the layer.

Values of the mobilisation indices and influence factors for all of the case A, B, C and D analyses are listed in Table 3. For case A the lower soil is arbitrarily considered to be ‘layer’ and the upper soil ‘matrix’. It can be seen that the mobilisation factor varies from less than 5% (nine out of 24 cases for χ_sd and four out of 24 cases for χ_ult), indicating that the layer has only a small influence, to over 30% (two out of 24 cases for χ_sd and four out of 24 cases for χ_ult), in which the strength and stiffness of the layer is clearly significant. In two cases, both involving a clay layer at or near the tip of a pile in sand, the computed effect of the layer at the ultimate displacement appears to be ‘negative’ (i.e. the overall capacity is enhanced despite H_layer being lower than H_matrix); these cases are discussed in further detail later.

To quantify the extent to which the PISA model is able to mimic the H–v_G response computed using 3D finite-analysis, an ‘accuracy metric’, η – consistent with previous work by Byrne et al. (2020a) – is employed. The accuracy metric is defined (Fig. 6) as

An accuracy metric of η = 1·0 indicates a perfect match between the 1D and 3D models. Two sets of η have been determined. A small displacement accuracy metric, η_sd, is evaluated for 0 ≤ v_G ≤ D/10 000; an ultimate displacement accuracy metric, η_ult, is determined for 0 ≤ v_G ≤ D/10.

A separate measure – referred to as the ‘ratio metric’ – is determined from the computed values of lateral load, H, at the reference ground-level displacements as

where H^(1D) and H^(3D) are values of lateral load computed with the 1D and 3D models, respectively (Fig. 6). Values of ratio metric are determined at the small and ultimate displacement reference pile displacements (ρ_sd and ρ_ult, respectively). Data on ratio and accuracy metrics for all of the analyses in the current study are listed in Table 4. For comparison purposes, averages of the ratio and accuracy metrics for all homogeneous soil calculations previously completed in connection with the development of the PISA design model (Byrne et al., 2019) are also provided in the table.

The current study had the twin aims of (a) exploring the performance of the PISA design model for a wide range of realistic soil layering conditions, and (b) identifying cases where the approach adopted in the current form of the model – in which homogeneous soil reaction curves are assigned to each of the layers – appears to have limitations.

Data on the full set of 35 analysis configurations considered in the study (Table 4) indicate that in all cases the small displacement accuracy metric, η_sd, is 0·9 or better, demonstrating that the 1D model matches closely the equivalent 3D finite-element calculation for v_G ≤ D/10 000. There is a greater variation in the ultimate displacement metrics; for the majority of the analyses (26 cases) η_ult ≥ 0·9; for eight cases 0·85 ≤ η_ult ≤ 0·89 and there is a single case (F5, discussed later) for which η_ult = 0·66.

Summary discussions are provided below on each of the layered soil configurations that have been considered. These discussions deliberately focus on cases where the PISA design model performed less well.

Four analyses employing the configuration in Fig. 4(a) were conducted, as specified in Table 5. In cases A1, A3 and A4, combinations of soils in which lower layers are stronger than upper ones (e.g. as quantified by the relative magnitudes of H_ult in Table 2) were employed, reflecting the likely conditions at most sites. Pile D2 was employed for cases A1 and A3; the thicker walled variant D2t was selected for cases A2 and A4 on the basis of initial assessments using the 1D model.

Values of H^(3D) computed at small and ultimate reference displacements are listed in Table 4. Corresponding data on mobilisation index, χ, and influence factor, IF, are listed in Table 3. The data on mobilisation index all lie in the range 0 < χ < 1, indicating that the computed response is intermediate between the relevant baseline data.

The PISA design model accuracy metrics, Table 4, indicate that cases A2, A3 and A4 all have both η_sd ≥ 0·95 and η_ult ≥ 0·95. The computed H–v_G response for case A2 is plotted in the bottom row in Fig. 7, indicating the closeness of the match between the 1D and 3D models. For case A1, although a close match with the 3D model is obtained at small displacements (η_sd = 0·98), the ultimate response accuracy metric (η_ult = 0·89) is not quite so good. Case A1 data are plotted in the top row of Fig. 7. Case A1 involves the combination Bothkennar clay and Dunkirk sand, D_R = 90%. These soils lie at opposite extremes of the strength range, as indicated in Fig. 5 and by the baseline data on H^(3D)_ult in Table 2; combinations of these two soils were consistently found in the current work to present a challenging test of the 1D model.

Figure 7 also provides comparisons of the bending moments in the embedded piles, computed from the 1D and 3D analyses, for cases A1 and A2. These bending moment data relate to the load (H^(3D)_ult) computed from the 3D finite-element calculation at v_G = D/10, and also H^(3D)_ult/2. The bending moment data computed from the 1D and 3D models are in close agreement. An exact match is naturally achieved at the ground surface, since in this comparison the 1D and 3D data both correspond to the same pile head moment.

The magnitudes of the bending moments induced in the pile are strongly conditioned by the pile head moment – which is determined solely by the applied loading. Consequently, subsurface bending moments computed with the PISA design model were found to provide a well-conditioned match with equivalent 3D finite-element data for all of the cases considered in the study. Bending moment data for other cases are therefore not further discussed in this paper.

Two configurations have been considered: a soft layer embedded in a relatively stiff matrix (case B1) and a stiff layer in a relatively soft matrix (cases B2 and B3). These cases are illustrated in Fig. 8 and defined in Table 6. The layer is located either in the upper half of the pile (mean depth L/3) or in the lower half (mean depth 2L/3). Two layer thicknesses are considered t_layer = L/15 and t_layer = L/4. Cases B1 and B3 involve combinations of the two least similar reference soils considered in the study (Bothkennar clay and Dunkirk sand, D_R = 90%). Case B2 comprises dissimilar clays: stiff layer (Cowden till) and soft matrix (Bothkennar clay). Initial analyses with the 1D model indicated that the influence of layers placed in the lower half of the pile is relatively small, and this is of course explained by the fact that a stiff pile tends to rotate under moment loading about a pivot point at a depth of about 2/3 of its length. Thin layer cases (t_layer = L/15) in the lower half of the pile were therefore excluded from the schedule of 3D finite-element runs.

Mobilisation indices and influence factors (Table 3) indicate, as expected, that incorporating a weak layer within a relatively strong matrix (case B1) causes a reduction in stiffness and strength of an embedded monopile (negative values of IF). In cases B2 and B3 (stiff layer in a soft matrix), incorporating a layer increases the stiffness and strength (positive values of IF). In all cases, when a thick layer is present in the upper half of the pile it has a greater influence on monopile performance than when it is in the lower half.

The PISA model provides a close match (Table 4) to the 3D finite-element data for combinations of Bothkennar clay matrix and Cowden till layer (case B2), and for Bothkennar clay and Dunkirk sand. In every case the metrics η_sd and η_ult are 0·90 or greater, except for B1BK only (η_ult = 0·85). Example data, for cases B2UN and B1BK, are shown in Fig. 9.

The configurations that have been considered are illustrated in Fig. 10 and specified in Table 7. The configuration in Fig. 10(a) (e = L/8) represents a monopile with a significant embedment in the base layer. For e = 0 (Fig. 10(b)), the pile base coincides with the upper boundary of the base layer. In this configuration, the base force and moment components are determined in the PISA design model by the base layer material; the distributed load and moment components are determined by the matrix. The configuration e = −L/8 (Fig. 10(c)) is a special case, incorporated in the study to investigate whether a deep-lying layer is ‘felt’ by the pile in the 3D finite-element calculations. The PISA design model does not include any representation of deep layers with a clearance below the pile tip; the e = −L/8 cases are therefore equivalent to matrix-only analyses in the model.

The ultimate mobilisation indices, χ_ult, for cases C2T and C2B (Table 3) have small negative values, signifying in these cases that the ultimate capacity of the monopile exceeds the capacity for the separate cases in which the monopile is embedded in homogeneous soils corresponding to the layer and matrix materials. This apparently paradoxical result has been investigated as follows. Pile D2t has a lower ultimate capacity when embedded in Cowden till (the layer material) than when embedded in D_R = 45% sand (the matrix material). However, the horizontal force and moment developed at the base of the pile at v_G = D/10 are actually greater when the pile is embedded in Cowden till than when embedded in D_R = 45% sand. In configurations C2T and C2B, therefore, the strength of the pile is slightly enhanced by the presence of the Cowden till base layer. Additionally C2 is something of a special case since the values of H_layer and H_matrix are numerically similar, with the consequence that the denominator in the definition of the index (equation (1)) is relatively small, thereby amplifying the value of the index.

Influence factor data in Table 3 indicate a substantial uplift (IF_ult = 42·2%) in the ultimate capacity for C1E (matrix = Bothkennar clay, base layer = Dunkirk sand, D_R = 90%) due to embedment into the stiff base layer. Case C3E (matrix = Bothkennar clay, base layer = Cowden till) also indicates a significant ultimate capacity uplift (IF_ult = 11·6%). Monopile embedment in a stiff base layer is clearly beneficial for the ultimate capacity of the monopile (and to a lesser extent the initial stiffness) in these cases. For case C1, the stiff base layer has a significant influence on the ultimate capacity even when the top of the layer lies below the pile base (case C1B, IF_ult = 11%).

Ultimate displacement accuracy metric data (Table 4) indicate that cases C2 and C3 all have η_ult ≥ 0·95. However, the C1 cases (layer = Dunkirk sand, D_R = 90%; matrix = Bothkennar clay) all have 0·87 ≥ η_ult ≥ 0·89; in these cases the 1D model systematically underestimates the ultimate displacement response as computed with the 3D model (ρ_ult < 1). This tendency presumably reflects the fact that the 1D model does not incorporate the full beneficial effect of the stiff material below the pile base. Example data, for C1E and C1T, are plotted in Fig. 11.

The case D analyses consider the combined effect of an embedded layer and pile embedment in a stiff material (Dunkirk sand, D_R = 90%), with a relatively soft matrix (Bothkennar clay), Fig. 12 and Table 8. The two cases considered comprise a combination of C1E with B3UN (case DN) and C1E with B3UK (case DK).

The influence factors for cases DN and DK are comparable in all cases to the sum of the relevant influence factors from the corresponding B and C cases (Table 3).

Figure 13 shows results for both D cases; for DK η_ult = 0·86, but the other accuracy metrics are all 0·95 or better (Table 4). The D cases were devised specifically to further test the capabilities of the 1D model for highly dissimilar soils and it is noteworthy that the model performs well in these configurations, especially for case DN.

Table 9 and Fig. 14 specify the analyses that were conducted.

Cases E1A and E1B were intended to test the capabilities of the general Dunkirk sand model, which provides a general sand calibration in the relative density range 45% ≤ D_R ≤ 90%. Consistent with likely density variations at offshore sites, the relative density in both E1 cases increases with depth.

Case E2 is a three-layer system employing all of the reference soil models. Cases E3 and E4 are more complex systems, proposed by the industry partners, to resemble cases encountered in typical design applications.

For all of these multi-layer systems, the 1D model performed well, with ultimate displacement metrics all in the range η_ult ≥ 0·95 (Table 4). Example data, for E2 and E4, are plotted in Fig. 15.

The analyses that have been conducted are specified in Table 10.

The small displacement accuracy metrics η_sd obtained for these cases are comparable to the equivalent cases with the L/D = 4 piles (Table 4). In all cases, however (except for F4), the ultimate response accuracy metrics for the F cases are lower than those from the equivalent cases A to E. This is especially the case for F5, for which the ultimate displacement metric is η_ult = 0·66. This is the lowest value of accuracy metric obtained in the entire study, and indeed the only value below 0·85. Data for F3 and F5 are plotted in Fig. 16.

The paper describes a set of 35 3D finite-element analyses on monopiles embedded in a variety of layered soil configurations and 12 reference calculations for homogeneous soils. The layered soil analysis results are all broadly consistent with expectations. A soft layer in a stiff matrix reduces the strength and stiffness of an embedded monopile. Conversely, incorporating a stiff layer in a soft matrix causes the strength and stiffness of an embedded monopile to increase.

Cases are examined in which a stiff layer is incorporated below the level of the pile base (cases C1B, C2B and C3B). In case C1B (matrix = Bothkennar clay; layer = Dunkirk sand, D_R = 90%), the stiff deep-lying layer has a significant influence on the computed ultimate strength of the monopile (uplift of 11% computed with the 3D model). These two soil types (Bothkennar clay and Dunkirk sand, D_R = 90%) have highly contrasting strength and stiffness characteristics. Significant strength uplift was not observed for cases C2B and C3B (involving combinations of Dunkirk sand, D_R = 45%/Cowden till and Bothkennar clay/Cowden till). It appears that an uplift in monopile capacity due to a deep-lying layer is only available if the layer material is very substantially stronger/stiffer than the matrix. The potential benefit of a deep-lying stiff layer is not incorporated in the current implementation of the PISA design model.

Assessments of the performance of the PISA design model are based on the numerical values of the ratio and accuracy metrics, η and ρ. For all small displacement cases, these metrics indicate a close match between the PISA design model and corresponding 3D finite-element analysis (η_sd ≥ 0·9 and 0·91 ≤ ρ_sd ≤ 1·08).

In 26 of the cases considered a similarly close match between the PISA design model and the corresponding 3D finite-element analysis for ultimate displacements is obtained (η_ult ≥ 0·9 and 0·89 ≤ ρ_ult ≤ 1·07). In eight other cases, however, the ultimate displacement match is less close (0·9 > η_ult ≥ 0·85 and 0·88 ≤ ρ_sd ≤ 1·22). A single case (F5, η_ult = 0·66, ρ_ult = 1·48) is an outlier. The cases with η_ult < 0·9 all involve combinations of Bothkennar clay and Dunkirk sand, D_R = 90%. Combinations of soft clay and very dense sand – at opposite ends of the strength spectrum – present a challenge to the PISA modelling approach in its current form.

Certain patterns are apparent in the cases with η_ult < 0·9. The computed lateral load at ultimate displacements is underestimated by the 1D model in five cases (A1, C1E, C1T, C1B, F4). These cases have the common characteristic that they incorporate very dense Dunkirk sand base layers. In these cases the sand layer appears to constrain the softer matrix soil – to increase its effective strength – in a way that is not captured by the PISA model.

In the four other cases with η_ult < 0·9 (including the worst case), the 1D model overestimates the ultimate capacity of the monopile (B1BK, DK, F3, F5). Three of these cases (DK, F3, F5) involve a thick layer of very dense Dunkirk sand in the upper half of the pile. It appears that in these cases the presence of the softer soil above and below the stiff layer has the effect of degrading the strength of the layer. The current 1D model does not include a means of incorporating this degraded strength effect. In case B1BK the configuration involves a soft embedded layer (Bothkennar clay) in a very dense Dunkirk sand matrix. In this case, it is considered that the presence of the soft embedded layer degrades the performance of the sand matrix material in a way that is captured by the 3D model but not the 1D analysis.

For some combinations of soft clay and very dense sand, however, the ultimate response performance of the 1D model does give η_ult ≥ 0·9 (i.e. all of the B3 cases and DN). In case DN the tendency for the 1D model to underestimate the influence of the stiff base layer may have compensated for the tendency to overestimate the influence of the embedded stiff layer. The B3 cases are perhaps anomalous; the performance of the 1D model is better than might be expected on the basis of the other soft clay/very dense sand combinations.

The case F calculations provide a means of comparing the performance of the 1D model for low embedment ratio monopiles (L/D = 2) with the larger embedment ratio monopiles (L/D = 4) that form the basis of the main study. The small displacement metrics for the F cases are comparable to those from the main study. The ultimate displacement case F metrics are, however, lower than the equivalent L/D = 4 analyses (with the exception of F4). Piles with L/D = 2 lie at the lower end of the parameter space adopted to calibrate the reference soils. At this relatively low embedment ratio, the lateral reactions reduce in relative importance and the other soil reaction components become increasingly significant (Burd et al., 2020a; Byrne et al., 2020a). These factors combine (as addressed in the discussion in Burd et al. (2020a) and Byrne et al. (2020a) for piles embedded in homogeneous soils) to cause a tendency for the accuracy of the 1D model to be lower for low embedment ratio piles. Improved performance may be obtained by restricting the homogeneous soil calibration process to a tighter geometrical space. Losses in accuracy of the 1D model due to the separate influences of low embedment ratio and layers with contrasting strength tend to reinforce each other, as is evident in the relatively low ultimate displacement accuracy metric for case F5.

Figure 17 indicates the accuracy metrics for the full set of layered soil analyses presented in this paper, together with the accuracy metrics for previous analyses on homogeneous soils (Byrne et al., 2019). It is apparent that the PISA design model small displacement response provides a close match with the corresponding 3D analyses in all layered and homogeneous cases. The ultimate displacement performance of the PISA model falls below the general small displacement trend for only a small number of cases, mostly for the piles with L/D = 2.

The mean value of the ratio metrics for all the layered soil analyses is close to 1 for both the small and ultimate displacement cases (Table 4); this is similar to the previously reported data (also listed in Table 4) for homogeneous soil analyses conducted using the PISA design model (Byrne et al., 2019). The coefficient of variation (CoV) of the small displacement ratio metric, ρ_sd, for the current layered study (3·9%, Table 4) is consistent with previous homogeneous soil calculations, indicating that the PISA design model is able consistently to deliver small displacement analyses of layered soil cases at the same level of fidelity as homogeneous cases. For ultimate displacement cases, the slightly larger CoV of the ratio metric, ρ_ult, compared to the homogeneous cases (Table 4) is driven largely by a few challenging configurations (e.g. F5, F3, B1BK) involving combinations of dense sand and soft clay.

The PISA design model, suitably calibrated for homogeneous soils following the methodology reported in Byrne et al. (2020a) and Burd et al. (2020a), captures the fidelity of 3D finite-element analyses of monopile behaviour, but at a fraction of the computational cost. This enables a wide range of scenarios to be explored during the design phase for an offshore wind farm. The current paper describes an extension of the model to layered soils, relevant to many offshore wind farm sites. This extended form of the model is shown to provide close representations of the behaviour of an embedded monopile, as computed with 3D finite-element analysis, for a wide range of layered soil cases.

The study suggests that for many practical layered soil configurations the PISA modelling approach can be used with high confidence to conduct design calculations. Challenging cases, where the model performs less well, involve specific combinations of very soft clay and very dense sand. In all design applications employing the PISA design model it is recommended that the final designs are checked using bespoke 3D finite-element analysis. This is especially the case when combinations of soft clay and very dense sand (e.g. the conditions in case F5) are present at the site.

Various options are available to develop the PISA approach further for cases where substantial differences in strength and stiffness exist in adjacent soil layers. One option would be to assume that the pile is only embedded in the weaker material (i.e. ignore the stronger material completely), although this is likely to be excessively conservative. Another possibility is to develop rules on how the strength of stronger layers is degraded by nearby softer soils. A consideration of these possibilities could be developed in a straightforward way in the future, most probably motivated by a site-specific application of the modelling approach rather than for generic application.

The PISA project was funded through the joint industry Offshore Wind Accelerator (OWA) programme, designed and led by the Carbon Trust, with PISA phase 1 supported by the UK Department for Energy and Climate Change (DECC) and PISA phase 2 supported by the Scottish government. The authors acknowledge the provision of financial and technical support by the following project partners: Ørsted (formerly DONG Energy, lead partner), E.ON, EDF, GE Renewable Energy (formerly through Alstom Wind), Iberdrola, innogy, SSE, Statkraft (phase 1 only), Equinor (formerly Statoil), Van Oord and Vattenfall.

Full details of the PISA design model formulation and implementation are given in Byrne et al. (2020a) and Burd et al. (2020a). Brief details are provided here.

The soil reaction curves employed in the model are represented in normalised form by the four-parameter conic function illustrated in Fig. 18. Separate frameworks are employed for clay and sand soil types (Table 11). Variations of the calibration parameters with depth (known as ‘depth variation functions’) for the three reference homogeneous soils employed in the current work are listed in Table 12.

The monopile is represented in the model by Timoshenko beam theory employing a shear factor κ = 0·5. Young's modulus and Poisson's ratio are 200 GPa and 0·3, respectively. Second moment of area and cross-sectional area are determined using thin-walled theory.

The data provided here are employed to formulate the layered ground models for incorporation in the PISA design model and the 3D finite-element analysis (see Table 13 and Fig. 19).

D

pile diameter

D_R

relative density

G₀

small-strain 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^(3D)_sd

lateral load applied to pile, computed with the 3D model at ground-level pile displacement v_G = D/10 000

H^(3D)_ult

lateral load applied to pile, computed with the 3D model at ground-level pile displacement v_G = D/10

h

load eccentricity

IF_sd

influence factor, ground-level pile displacement v_G = D/10 000

IF_ult

influence factor, ground-level pile displacement v_G = D/10

K₀

earth pressure coefficient at rest

k

initial stiffness of parametric soil reaction curve

L

pile embedded length

M

bending moment

M_B

moment at pile base

m

distributed moment acting on pile

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_G

ground-level lateral pile displacement

$x¯u$

ultimate displacement for parametric soil reaction curve

$y¯u$

ultimate load for parametric soil reaction curve

z

depth coordinate along the pile

η

accuracy metric

η_sd

small displacement accuracy metric, ground-level pile displacement up to D/10 000

η_ult

ultimate displacement accuracy metric, ground-level pile displacement up to D/10

ρ

ratio metric

ρ_sd

small displacement ratio metric, ground-level pile displacement v_G = D/10 000

ρ_ult

ultimate displacement ratio metric, ground-level pile displacement v_G = D/10

$σ′vi$

initial vertical effective stress

χ

mobilisation index

χ_sd

mobilisation index, ground-level pile displacement v_G = D/10 000

χ_ult

mobilisation index, ground-level pile displacement v_G = D/10

ψ

rotation of the pile cross-section