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

The PISA design model is a procedure for the analysis of monopile foundations for offshore wind turbine applications. This design model has been previously calibrated for homogeneous soils; this paper extends the modelling approach to the analysis of monopiles installed at sites where the soil profile is layered. The paper describes a computational study on monopiles embedded in layered soil configurations comprising selected combinations of soft and stiff clay and sand at a range of relative densities. The study comprises ( a ) analyses of monopile behaviour using detailed three-dimensional (3D) finite-element analysis, and ( b ) calculations employing the PISA design model. Results from the 3D analyses are used to explore the various influences that soil layering has on the performance of the monopile. The fidelity of the PISA design model is assessed by comparisons with data obtained from equivalent 3D finite-element analyses, demonstrating a good agreement in most cases. This comparative study demonstrates that the PISA design model can be applied successfully to layered soil configurations, except in certain cases involving combinations of very soft clay and very dense sand.


INTRODUCTION
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 structuresmeasured by way of supervisory control and data acquisition (Scada) instrumentationare 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 Byrne et al., 2020a). This model, also referred to below as the 'onedimensional (1D) model', was an outcome of a projectknown as PISAthat incorporated ground characterisation , field testing Byrne et al., 2020b;McAdam et al., 2020), threedimensional (3D) finite-element analysis Zdravković et al., 2020b) and design model development 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 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 soilsas computed using 3D finite-element analysisfor 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 largediameter 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' ; (b) a marine sand known as 'general Dunkirk sand model' (GDSM) with relative density D R in the range 45% D R 90% ; 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).

MODELLING PROCEDURES PISA design model
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 crosssection 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,ȳ ¼ȳð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,ȳ u ,x u , n; these parameters in general vary with depth; functions describing these spatial variationsreferred 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 representedto a high fidelityby 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 0 , 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 . 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 model
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 . 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  as well as forming the basis of the calibration process for the general Dunkirk sand model . 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

Reference displacements for monopile performance metrics
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.

Reference homogeneous soil models
Profiles of small-strain shear modulus G 0 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 0 adopted in these models are specified in Appendix 2.

Monopile configurations
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.

Layered soil configurations
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 ator just belowthe 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 analysescase Fselected cases are re-analysed for a monopile with L/D ¼ 2 (piles C1 or C1t). 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.

Layered soil model implementation
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 0 employed in the reference homogeneous soil models (Table 13 and Fig. 19).
Procedures are needed to determine the profile of smallstrain shear modulus G 0 and, for clay layers, the profile of undrained shear strength, s u . Although G 0 and s u are specified in the representative models as functions of depth (Fig. 3), G 0 and s u are in practice dependent on the magnitude of the local stresses. In the current layered soil models, therefore, G 0 and s u are considered to be functions of local mean effective stress, G 0 ( p′), s u ( p′). To implement these functional relationships in the analyses, pairs of numerical data [G 0 , 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 0 ( p′), s u ( p′). Piecewise linear functions G 0 (z), s u (z)for incorporation in the layered soil analysis modelsare 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.

THREE-DIMENSIONAL FINITE-ELEMENT BASELINE ANALYSES
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.

Influence of soil layers
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.

Comparisons between the PISA design model and equivalent 3D finite-element analyses
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 measurereferred 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 For case A the upper soil is designated as 'matrix' and the lower soil as 'layer'. Subscript 'sd' signifies small displacement and 'ult'

LAYERED SOIL ANALYSES Overview
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 layersappears 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.

Case A: two-layer systems
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. Table 4. Values of horizontal force H (3D) determined from the 3D finite-element layered soil analyses at the two reference displacements (small (v G = D/10 000) and ultimate (v G = D/10)). Values of ratio metric ρ and accuracy metric η are listed at the two reference displacements

Case
Code Pile 3D finite element 1D/3D model comparisons 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 ult (3D) 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 ult (3D) ) computed from the 3D finite-element calculation at v G ¼ D/10, and also H ult (3D) /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 momentwhich 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.

Case B: single embedded layer
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,   (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.

Case C: base layer
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. Forcase 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.

Case D: conceptual multi-layer systems
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).
Note: the dimension e represents the base embedment (see Fig. 10).  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.
Case E: exploratory multi-layer systems 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.
Case F: selected repeated cases with a low embedment ratio monopile (L/D = 2) 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.

DISCUSSION
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 sandat opposite ends of the strength spectrumpresent 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 soilto increase its effective strengthin 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 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

CONCLUSIONS
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. 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. APPENDIX 1. FORMULATION DETAILS FOR PISA DESIGN MODEL 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. Table 11. Non-dimensional forms employed in PISA design model; s u is local value of triaxial compression undrained shear strength; G 0 is local value of small-strain shear modulus; σ′ vi is local value of initial vertical effective stress For the GDSM, relative density is expressed as a decimal in the calibration range 0·45 D R 0·9.