Finite-element modelling of laterally loaded piles in a stiff glacial clay till at Cowden

The PISA project was a combined field testing/numerical modelling study with the aim of developing improved design procedures for large-diameter piles subjected to lateral loading. This paper describes the development of a three-dimensional finite-element model for the medium-scale pile tests that were conducted in Cowden till as part of the PISA work. The paper places particular emphasis on the consistent interpretation of the soil data determined from the available field and laboratory information. An enhanced version of the modified Cam clay model was employed in the numerical analyses, featuring a non-linear Hvorslev surface, a generalised shape for the yield and plastic potential surfaces in the deviatoric plane and a non-linear formulation for the elastic shear modulus. Three-dimensional finite-element analyses were performed prior to the field tests. Excellent agreement between the measured and simulated behaviour for a range of pile geometries was observed, demonstrating the accuracy of the numerical model and the adequacy of the calibration process for the constitutive model. The developed numerical model confirmed the premise of the PISA design method that site-specific ground characterisation and advanced numerical modelling can directly facilitate the development of additional soil reaction curves for use in new design models for laterally loaded piles in a stiff clay till.


INTRODUCTION
The design of laterally loaded piles in clay is a complex problem, involving aspects of soil-structure interaction which are difficult to consider in most methods of analysis, such as the possibility of a gap opening on the active soil side. However, despite these challenges, there are various published studies investigating the behaviour of laterally loaded piles, with solution strategies including elastic methods (e.g. Poulos & Davis, 1980;Basu et al., 2009) and numerical approaches, such as those based on the finite-element (FE) method (e.g. Brown & Shie, 1990;Yang & Jeremić, 2002;Jung et al., 2015;Haiderali & Madabhushi, 2016). In these studies, the response of clay is typically simulated using a total stress approach with a Tresca (e.g. Jung et al., 2015) or a von Mises failure criterion (e.g. Brown & Shie, 1990;Yang & Jeremić, 2002). The undrained shear strength, S u , is often assumed to be either uniform (e.g. Yang & Jeremić, 2002) or increasing with depth (e.g. Jung et al., 2015), while the associated Young's modulus is specified as being uniform (e.g. Brown & Shie, 1990), a function of S u (e.g. Jung et al., 2015), or to depend directly on the confining pressure (e.g. Yang & Jeremić, 2002). Although the abovementioned studies provide significant insight into the response of laterally loaded piles in clay, particularly in terms of ultimate limit states, the accurate prediction of displacements requires the use of more comprehensive modelling approaches, capable of simulating accurately the non-linear pre-yield response of the soil, including the dependency of stiffness on the stress and strain levels (see e.g. Jardine et al. (1986) for a discussion on the impact of this aspect of soil response on the behaviour of geotechnical structures). In this paper, the procedure adopted in the advanced three-dimensional (3D) FE analyses carried out in the context of the PISA project, which include the abovementioned features, is presented and discussed.
As outlined in the paper by Zdravković et al. (2019), the PISA (Pile-Soil Analysis) project proposed a new design model for laterally loaded piles. Referred to as the PISA design model, it is consistent with the existing onedimensional (1D) Winkler-type p-y approach, but extended to include soil reactions in addition to that representing just a distributed lateral load. The principal premise of this development was that the new 1D PISA design model could be derived from the results of site-specific 3D FE modelling. Consequently, the 3D FE analyses discussed in the current paper were not back-analyses of pile tests. Instead, they were performed before the medium-scale piles, installed in a glacial clay till at the Cowden test site, were tested under lateral loading, thus providing Class A predictions according to Lambe (1973). The results of these analyses were used to: (a) initially aid the design of the field testing programme ; (b) subsequently validate the established numerical model by comparison with field testing results; and (c) provide insight into the main load-transfer mechanisms of laterally loaded monopiles in glacial clay, which would form the basis for formulating soil-reaction curves for new Winkler-type design models. Similar field testing  and numerical  studies were performed on piles installed in a dense marine sand at the Dunkirk test site.
The central focus of the PISA project was to develop a design model for monotonic loading only (see Zdravković et al., 2019), as a key step for future extensions of the model to cyclic and other loading conditions. Consequently, the 3D FE model and most of the field testing programme were developed with this objective. Some field testing explored the effects of different loading rates and of load cycles, to maximise the output from the field tests, but these aspects were not considered at the current stage of the new design model development.
The paper starts with elaborating a consistent 3D numerical model for laterally loaded piles installed at the Cowden test site, based on available ground information interpreted and discussed in the paper by Zdravković et al. (2019). This is followed by demonstrating the capabilities of the developed numerical model through comparison of the predicted response with the field measurements presented in the paper by Byrne et al. (2019).

CONSTITUTIVE MODEL
The interpretation of the Cowden ground model discussed in the paper by Zdravković et al. (2019) demonstrated the applicabilityof the critical state framework to describe the mechanical behaviour of the low-plasticity, overconsolidated clay till. The constitutive model employed for the numerical analyses presented here is selected to (a) satisfy design requirements, as outlined in the paper by Zdravković et al. (2019), for realistic representation of initial and ultimate loading behaviour of the test piles, and (b) simulate accurately aspects of the mechanical behaviour of Cowden till, observed from the available historic and new field and laboratory investigations. The model adopts elements of the original modified Cam Clay model formulation (Roscoe & Burland, 1968), but is enhanced with: (a) a nonlinear Hvorslev-type surface (Tsiampousi et al., 2013) to capture accurately the undrained strength of Cowden till dry of critical; (b) a generalised shape for the yield and plastic potential surfaces in the deviatoric plane (Van Eekelen, 1980) to account for the effect of the intermediate principal stress, as evidenced by the different strength of Cowden till in triaxial compression and extension; and (c) a modified hyperbolic equation for simulating the non-linear variation of the elastic shear modulus with mean effective stress and deformation level (Taborda & Zdravkovic, 2012;Taborda et al., 2016). The latter component of the numerical model can be coupledwith a reversal-detection algorithm and non-linear scaling laws (see Taborda et al., 2016) in order to replicate the hysteretic response exhibited by geomaterials when subjected to cyclic loading. However, given that the objective of the PISA project was to establish the response of piles subjected to monotonic lateral load, this aspect of the constitutive model formulation was deemed unimportant to the analyses presented herein and was deactivated.
It is noted that the selected constitutive model does not reproduce cyclic and strain-rate effects, as this was not the objective of the project, nor does it reproduce strength anisotropy of the soil, as characterisation of these aspects of Cowden till behaviour was not possible from the available experimental information.

Formulation
The Hvorslev surface, presented in detail by Tsiampousi et al. (2013), introduces alterations to the shape of the yield and plastic potential surfaces in the p′-J plane (where p′ is mean effective stress and J is generalised deviatoric stress, see Appendix 1) on the dry (or supercritical) side of the critical state, thus defining the elastic region for overconsolidated clays and controlling their plastic response. The adopted yield surface, F ðσ′; p′ 0 Þ, is defined by where θ is the Lode's angle (see Appendix 1); gðθÞ defines the inclination of the critical state line (CSL) to the p′ axis for the current value of θ; α and n are parameters defining the non-linearity of the yield surface on the dry side (i.e. p′ , p′ 0 =2); and p′ 0 is the pre-consolidation mean effective stress and the state/hardening parameter for this model. Fig. 1 illustrates the shape of the yield surface defined by equation (1) in the p′-J plane. The plastic potential, Pðσ′; p′ 0 Þ, is modified accordingly where β and m are parameters defining the shape of the plastic potential and therefore control the rate at which the material reaches the critical state, while p′ c and θ c denote the current mean effective stress and Lode's angle, respectively (see Tsiampousi et al., 2013). As shown in the expressions above, on the wet side of critical (i.e. p′ ! p′ 0 =2) the yield surface and plastic potential are similar to the original proposal by Roscoe & Burland (1968), where associated plasticity was assumed.
As shown in Fig. 1, the CSL inclination is controlled by the function gðθÞ, which defines the shape of the yield surface in the deviatoric plane (see Potts & Zdravković (1999) for a detailed discussion on the importance of this function). In the present case, the generalised shape proposed by Van Eekelen (1980) was employed  Lade & Duncan (1975) surfaces.
The final modification to the formulation is the introduction of a non-linear variation for the tangent shear modulus, which is based on the modified hyperbolic expression used in the ICG3S model described in the papers by Taborda & Zdravkovic (2012) and Taborda et al. (2016). According to this, the tangent non-linear elastic shear and bulk moduli are given by where G Ã 0 and K Ã 0 are the basic stiffness properties of the material (i.e. without taking into account the effects of void ratio, mean effective stress and strain level); p′ ref is a reference pressure; f G ðeÞ and f K ðeÞ are user-defined functions expressing the influence of void ratio on the values of the shear and bulk moduli, respectively; m G and m K are parameters controlling the non-linearity of the variation of stiffness with mean effective stress; E d is the generalised deviatoric strain (i.e. second invariant of the strain tensor, see Appendix 1); ε vol is the volumetric strain (i.e. first invariant of the strain tensor, see Appendix 1); and a; b; R G;min , r, s and R K;min are parameters defining the reduction in stiffness with deformation level. Clearly, there is an incompatibility between the ICG3S model and the original formulation of the modified Cam Clay model, since the assumption that the swelling lines have a constant slope (κ) in the plane plotting the specific volume (v) against ln p′ implies that the value of the tangent elastic bulk modulus is Therefore, the use of this small-strain stiffness overlay is only possible if the adopted function of void ratio is defined as f K e ð Þ ¼ 1 þ e ¼ v and m K and R K;min are set to 1·0 (the latter parameter renders the values of r and s irrelevant). In this situation, the ICG3S model will reproduce swelling lines of constant slope, and hence be compatible with the original formulation of the modified Cam Clay, providing the value of K Ã 0 is chosen to be As a result, the integrity of the assumptions central to the modified Cam Clay model is guaranteed, meaning that, for example, the closed-form solutions for the undrained shear strength (Potts & Zdravković, 1999) remain unchanged. Concurrently, the use of the ICG3S model enables a more accurate representation of the considerably non-linear variation of the shear modulus exhibited by Cowden till , which would not be captured by the alternative formulations for the elastic shear modulus typically employed when using the modified Cam Clay model (constant Poisson ratio, ν, constant shear modulus, G, or constant ratio G=p′ 0 ).

Calibration
The interpretation of the historic and newly performed laboratory tests on Cowden till in terms of strength at critical state under triaxial compression and triaxial extension loading conditions is discussed in detail in the paper by Zdravković et al. (2019). These data indicate angles of shearing resistance of ϕ′ TXC ¼ 27°and ϕ′ TXE ¼ 32°, in triaxial compression and extension, respectively. Since no information on the strength of the material is available for intermediate values of the Lode's angle, a number of combinations are possible for the parameters X, Y and Z in equation (3), ensuring the convexity of the yield surface. To achieve the experimentally observed difference of up to 6°in the angle of shearing resistance between triaxial and plane-strain loading conditions for most soils (e.g. Bishop, 1966;Gens, 1982), an iterative procedure suggested that it was appropriate to assume a value of 0·1 for Z, thus allowing the calculation of the remaining quantities using In the expressions above, M J;TXC and M J;TXE designate the ratios J=p′ at critical state for triaxial compression (θ ¼ À30 o ) and triaxial extension (θ ¼ 30 o ) loading conditions, respectively. These can be calculated using (Potts & Zdravković, 1999) For the calibrated parameters, the variation of the angle of shearing resistance with the value of Lode's angle is illustrated in Fig. 2(a). Note that higher values of Z would have increased substantially the simulated strength of the material for non-triaxial loading conditions, which would be unrealistic.
The characterisation of the strength of the material requires the calibration of the parameters controlling the shape of the non-linear Hvorslev surface. These were determined by plotting the measured effective stress paths, of highly overconsolidated samples sheared under undrained and performing a non-linear regression, which gave α ¼ 0Á25 and n ¼ 0Á40. The calibrated non-linear Hvorslev surface is plotted together with the normalised stress paths in Fig. 2(b), demonstrating good agreement between the measured and simulated yield loci. The results of the oedometer tests, described in the paper by Zdravković et al. (2019), indicate that the consolidation behaviour of Cowden till, particularly that of the more superficial layers which are dominant for shorter piles, can be represented with good accuracy by ν 1 ¼ 2Á20 and λ ¼ 0Á115. The swelling response, as measured in the same tests, is characterised by κ ¼ 0Á021, thus determining the tangent bulk modulus (equation (6)). Since the latter is linearly dependent on the value of mean effective stress, a similar variation with p′ was assumed for the elastic shear modulus, meaning that m G ¼ 1Á0 in equation (4).
A summary of field and laboratory measurements of the Cowden till stiffness is presented in the paper by Zdravković et al. (2019), including results from the new seismic cone penetration tests (SCPTs), bender element (BE) tests and triaxial compression (TXC) and triaxial extension (TXE) tests using high-resolution linear variable differential transducers (LVDTs), as well as from historic cross-hole (CH) and downhole (DH) geophysics. As discussed in the paper by Zdravković et al. (2019) and shown in Fig. 3(a), the dynamic measurements SCPT, BE and DH indicate considerably higher stiffness values than those obtained using local transducers (TXC, TXE), which may be attributed to insufficient resolution of the latter to capture elastic stiffness at very small strains. Given the consistency exhibited by the dynamic measurements, a value of G 0 =p′ ¼ 1100Á0 was assumed representative of the small-strain shear modulus. As discussed in the paper by Zdravković et al. (2019), a much higher stiffness profile from the historic CH data might be attributed to soil variability at the previous Cowden test area (described in the paper by Powell & Butcher (2003)), which was substantially larger than the PISA test area and located in a different part of the Cowden site (see Zdravković et al., 2019). Adopting a reference pressure, p′ ref , of 100 kPa and discarding the effect of void ratio on the value of the elastic shear modulus due to lack of clear experimental evidence (i.e. f G e ð Þ ¼ 1Á0), implies that parameter G Ã 0 in equation (4)    The measured variations of the normalised secant shear modulus, G sec =p′, with strain level, E d , for individual triaxial tests are shown in Fig. 3(b), together with the simulated behaviour obtained using equation (4), with parameters a ¼ 9Á78 Â 10 À5 , b ¼ 0Á987 and R G;min ¼ 0Á05. As seen in Fig. 3(a), there are a number of triaxial tests for which the maximum stiffness is below that given by the assumed relationship G 0 ¼ 1100p′. Therefore, the values chosen for a, b and R G;min in most cases imply higher-than-measured stiffness for relatively small strains. To mitigate the impact of this overestimation, the non-linear elastic part of the constitutive model was calibrated to predict a considerably sharper reduction in stiffness with increasing deviatoric strain than observed in most tests, such that at moderate strains the simulated behaviour is closer to the average stiffness measured in all the triaxial tests. This is shown more clearly in the normalised G sec =G 0 degradation with deviatoric strain in Fig. 3(c), further highlighting the need to consider globally the stiffness during calibration, rather than focusing on producing, in isolation, the maximum stiffness ( Fig. 3(a)) and its normalised reduction curve (Fig. 3(c)).
The final parameters to be determined are those controlling the shape of the plastic potential function associated with the non-linear Hvorslev surface, β and m. These were evaluated using a trial-and-error approach, focusing on the rate at which overconsolidated samples approach the critical state, as discussed in Tsiampousi et al. (2013). In the present case, a simple linear form for the plastic potential function was assumed (i.e. m ¼ 1Á00), with a value of β ¼ 0Á20 appearing to reproduce the results of triaxial tests with greater accuracy, as shown in Appendix 2. The resulting set of calibrated parameters is listed in Table 1, while the performance of the model is briefly illustrated by comparisons with experimental data in Appendix 2.

NUMERICAL ANALYSIS Geometry and boundary conditions
All analyses reported herein were carried out with the Imperial College Finite Element Program (ICFEP) (Potts & Zdravković, 1999). Given the existence of a plane of symmetry in the problem of a laterally loaded pile, only half of the domain was discretised, as shown in an example FE mesh in Fig. 4. In total, 14 test piles were installed at Cowden , with four distinct geometries, as summarised in Table 2, representing medium-(D ¼ 0Á762 m) and large (D ¼ 2Á0 m)-diameter piles that were the focus of numerical predictions in the current paper. Multiple tests were conducted on the same pile geometries to confirm the repeatability of the field measurements.
The far vertical boundary was positioned at a radial distance of 30 m for analyses CM2, CM3 and CM9 (D ¼ 0Á762 m), while for analysis CL2, where a pile of larger diameter (2·0 m) was simulated, this was increased to 80 m, mitigating any possible boundary effects. In all analyses the depth to the bottom boundary was 40 m, at the interface of the Cowden till deposit with chalk. The tubular open-ended pile is discretised with 280 eight-noded shell elements (Schroeder et al., 2007), of which 18 rows of elements in the Z-direction discretise the embedded part of the pile and ten rows the stickup. The steel was assumed to be linear elastic, with a Young's modulus of E = 200 GPa and a Poisson ratio, ν = 0·30. Furthermore, the pile was assumed to be wished in place, hence disregarding any installation effects. This was considered acceptable as a laterally loaded pile derives its capacity from the mobilisation of a significant volume of soil around the pile, well beyond the pile-soil interface zone that Non-linear Hvorslev surfaceshape (Tsiampousi et al., 2013), equation (1) α ¼ 0Á25; n ¼ 0Á40 Non-linear Hvorslev surfaceplastic potential (Tsiampousi et al., 2013), equation (2) β ¼ 0Á20; m ¼ 1Á00 Virgin consolidation line ν 1 ¼ 2Á20; λ ¼ 0Á115 Non-linear elasticityswelling behaviour κ ¼ 0Á021 Non-linear elasticitysmall-strain shear modulus (Taborda et al., 2016), equation (4) Non-linear elasticityshear stiffness degradation (Taborda et al., 2016), equation (4) a  may be disturbed by installation. Consequently, the effect of the disturbed zone is thought to be less significant compared to an axially loaded pile, which relies on interface conditions for its capacity. As shown later, for lateral loading it is much more important that the interface can simulate opening of a gap around the pile, as this would affect its capacity. Consequently, the outer pile-soil interface is discretised using 180 16-noded zero-thickness interface elements (Day & Potts, 1994), the behaviour of which is represented by an elastoplastic Tresca model with elastic normal (K N ) and shear stiffness (K S ) values of 1Á0 Â 10 5 kN=m 3 and a limited tensile capacity (in the present case, assumed to be zero). The interface shear strength in compression, denoted as S u;int , is defined using an innovative approach where the interface element adopts the initial undrained shear strength of the adjacent soil element, thus seamlessly reflecting the contributions of the spatial variations in stress conditions, overconsolidation ratio (OCR) and loading direction (θ). For the modified Cam Clay model, an analytical solution for the undrained shear strength is utilised (Potts & Zdravković, 1999) where OCR is determined using σ′ vm =σ′ v0 , with σ′ vm being the maximum vertical effective stress that the soil has been subjected to and σ′ v0 being the initial vertical effective stress, K OC 0 is the current value of the coefficient of earth pressure at rest, while K NC 0 is its value associated to normally consolidated conditions and calculated as ð1 À sin ϕ′ cs Þ. Lastly, B is given by To avoid rigid body movements of the mesh, the displacements along all directions (X , Y and Z) are set to zero over the bottom boundary of the mesh (i.e. at Z ¼ À40Á0 m). Moreover, the displacements normal to the vertical cylindrical boundary are prescribed as zero. To ensure that the X-Z plane at Y ¼ 0 is a plane of symmetry, the displacements in the Y-direction over this plane are set to zero. Similarly, at the nodes of shell elements, defining the edges of the pile contained in this plane, the rotational degrees of freedom about the X-and Z-axes are also set to zero. The load applied at the top of the pile (i.e. Z ¼ h ¼ þ10 m) is simulated by imposing an incremental uniform horizontal displacement in the X-direction to all the nodes defining the perimeter of the section. The calculated reactions at each node are then summed up to give the total load H=2.

Ground conditions
A detailed description of the ground conditions is provided in the paper by Zdravković et al. (2019), while herein only information relevant to the numerical analyses is presented. In terms of initial conditions, the total vertical stress is obtained by assuming a saturated bulk unit weight, γ sat , of 21Á19 kN=m 3 , which is consistent with the values of density reported by Powell & Butcher (2003) and Ushev et al. (2015). The adopted variation with depth of the at-rest coefficient of earth pressure, K 0 , is shown in Fig. 5(a), with a profile that follows the Powell & Butcher (2003) data and imposes the shallow cut-off of 1·5, as discussed in the paper by Zdravković et al. (2019). The pore water pressures measured in situ revealed an under-drained distribution in the top 12 m, Fig. 5(b), as a result of the soil permeability, k, reducing considerably with depth and the presence of a lower aquifer with a phreatic surface below the groundwater table . The profile adopted in the analyses (solid line) is derived at the initiation of the ground conditions, using an exponential law for the permeability variation with mean effective stress (see Vaughan, 1994;Potts & Zdravković, 1999), fitted to reproduce the reported reduction in permeability while maintaining pore water pressures in equilibrium with the water table at 1 m depth and the values estimated in the sand layer. The resulting non-linear underdrained pore water pressure profile exhibits good agreement with field measurements. Below the sand layer the available data on pore pressures are limited and a hydrostatic variation is adopted. Clearly, as the value of permeability simultaneously depends on the mean effective stress and controls the pore water pressure profile, an iterative approach is required to establish the correct initial pore water pressure profile to be used in the analysis. However, it should be noted that the procedure followed was needed solely to ensure that the under-drained pore water pressure profile would represent the steady-state conditions and that the correct initial effective stress state was reproduced, as it controls not only the simulated undrained shear strength profile (equation (11)), but also the soil stiffness (equations (4) and (5)).
Since all the piles were wholly embedded in Cowden till (maximum length was 10·5 m, see Table 2), the soil profile adopted in all analyses is simplified to a single layer of this material, modelled in its entirety with the parameters listed in Table 1. The existence of the two deeper sand layers is ignored. Moreover, given the low permeability of Cowden till, the numerical analyses were performed as undrained, a condition achieved by specifying that the bulk modulus of the pore fluid was three orders of magnitude larger than that of the soil (see Potts & Zdravković (1999) for further details of this procedure).
As indicated by equation (11), the undrained shear strength simulated by the adopted constitutive model is a function of the stress state, the model parameters and the OCR. Therefore, based on the measured variation with depth of the undrained shear strength in triaxial compression, S u;TXC (Fig. 5(d); see Zdravković et al. (2019) for additional details), the OCR profile is derived, as shown in Fig. 5(c). The fact that the adopted OCR profile agrees very well with measured data is an indication of a high level of consistency between the calibrated model parameters and adopted ground conditions. This is further demonstrated in Fig. 5(e) by the very good agreement between the measured and predicted profiles of undrained strength in triaxial extension, S u;TXE , in the top 8 m of the deposit. As K 0 . 1 in this part of the ground, the initial soil state is in triaxial extension (TXE) and the predicted profile (grey line) is the result of the adopted deviatoric shape of the yield surface (Van Eekelen, 1980). The experimental data show a slightly larger scatter compared to the TXC tests in Fig. 5(d), but this is attributable to the usual issue of 'necking' failure, typical for triaxial testing in extension. The undrained strength profile at the start of 3D FE analysis is therefore a combination of the S u;TXE profile in the top 8 m of the ground and S u;TXC below 8 m.

General considerations
The influence of the geometric characteristics of the monopile on its response is illustrated in Fig. 6, where the ZDRAVKOVIĆ , TABORDA, POTTS ET AL.
normalised deformed shapes of piles CM2 and CM3, characterised by the two extreme values of L=D, are shown for a ground-level displacement, v G , in the direction of loading (X ) corresponding to 10% of their diameter (0Á1D ¼ 76Á2 mm). Note that v G corresponds to the displacement at the leading (front) edge of the pile. Throughout the PISA project, v G ¼ 0Á1D was assumed to correspond to the ultimate limit state to allow direct comparisons between the different piles. Although the two piles have the same diameter, they have different lengths, with values of L=D of 3·0 and 10·0, respectively. Pile CM2 exhibits a rigid-body rotation about a point located at a depth corresponding approximately to 70% of the pile length, while the response of pile CM3 is significantly more flexible. Also shown in Fig. 6 are the maximum depths of the gaps formed around the two piles as a consequence of the applied loading. The depth of the gap varies considerably around the pile circumference, with its maximum occurring in the plane of symmetry, on the active side of the pile. The rigid-body rotation of pile CM2 results in the formation of a deeper gap in normalised terms, L G =L ¼ 0Á67, than is the case for pile CM3, for which L G =L ¼ 0Á53, L G being depth of the gap. The significant extent of the gaps determined from the numerical model closely replicates the observed field behaviour Historic data (Robson, 1988) Piezometers (PISA) Simulated pore water pressure profile  , with measured L G =L of 0·72 and 0·48 for piles CM2 and CM3, respectively. This highlights the need to employ non-linear elasto-plastic interface elements, incorporating the possibility of gapping, in the simulation of laterally loaded piles, particularly for relatively low values of L=D. In the analyses discussed here, it is assumed that the gap remained dry after formation and therefore the tension cut-off was set to zero. However, this cut-off value can be set as depth-dependent, corresponding to a possible hydrostatic pore water pressure distribution, in order to simulate water-filled gaps, such as the ones expected in piles installed offshore.
The importance of adopting the Van Eekelen (1980) deviatoric shape for the yield surface, which is capable of reproducing the effect of loading conditions on the strength of the material, is demonstrated in Fig. 7(a). This figure shows the spatial distribution of the value of the Lode's angle, θ, obtained for pile CM2 at a ground-level displacement of v G ¼ 0Á1D. According to the adopted definition for θ (see Appendix 1), θ ¼ À30°corresponds to triaxial compression loading, whereas θ ¼ 30°characterises that in triaxial extension. As expected, the value of K 0 . 1 adopted in the top 8 m (Fig. 5(a)) implies that the material surrounding the pile is initially characterised by θ ¼ 30°. As the pile is pushed forward in the X -direction, the loading conditions applied to the soil evolve rapidly and substantially. In effect, at the chosen stage of v G ¼ 0Á1D, the superficial zone in front of the pile, where relatively large plastic strains occur, assumes values close to θ ¼ À30°in the immediate vicinity of the pile and therefore an operational strength of ϕ′ ¼ 27°. The magnitudes of θ increase with increasing radial distance in front of the pile, engaging again the initial θ ¼ 30°in the distant parts of the superficial soil domain that are unaffected by the pile loading. Furthermore, a considerable variation in the circumferential direction can be seen in Fig. 7(a), with the region immediately behind the pile being characterised by high values of θ, for which the angle of shearing resistance is 32°, caused by the lateral unloading of the material which takes place until the gap is formed.
Contours of mobilised deviatoric strains, E d , are illustrated for the v G ¼ 0Á1D stage in Fig. 7(b). Note that, in this case, a logarithmic scale is adopted for the contour lines. Clearly, in the soil surrounding the pile, the strain level, and hence the shear stiffness of the material (equation (4)), changes rapidly, further reinforcing the need for a robust modelling approach which is capable of simulating soil behaviour over the entire range of deformation levels that occur in the problem. Figure 8 compares the simulated embedded response of piles CM2 (L=D = 3) and CM3 (L=D = 10) with that interpreted from the field tests Byrne et al., 2019). Note that, for a meaningful comparison, the plots in Fig. 8 are obtained by selecting stages of the analyses where either the ground-level displacement (Figs 8(a) and 8(c)) or the ground-level moment (Figs 8(b) and 8(d)) match those of the field test. Therefore, it is important to highlight that the deformed shapes do not necessarily correspond to the same loading level, while the bending moment distributions are not calculated for an identical value of groundlevel displacement. Also shown in Fig. 8 are data, labelled 'optimised structural model', for the embedded pile response. This model was determined by obtaining a best fit between a structural model of the embedded pile and the measured strain gauge and inclinometer data. The process of determining this model is described in the paper by Burd et al. (2019).

Embedded pile response
The results in Fig. 8 demonstrate the excellent performance of the FE model, which is capable of reproducing with equal accuracy the measured deformations (Figs 8(a) and 8(c)) and bending moments (Figs 8(b) and 8(d)) for the two selected pile geometries, despite their significantly different values of slenderness ratio. More importantly, both the measurements and the model demonstrate a significant moment reaction at the toe of the rigid CM2 pile (Fig. 8(b)) and practically a zero value of this component at the toe of the flexible CM3 pile (Fig. 8(d)). The agreement between the numerical and field data confirms that a successful reproduction of a complex boundary value problem of this sort requires the accurate simulation of soil behaviour under the full range of strain levels, since ground movements appear to be concentrated in the top half of the pile, with the lower half imposing significantly smaller deformations on the soil.

Load-displacement response
The analyses results for the four chosen pile geometries in terms of load-displacement curves are compared in Fig. 9 with the observed field behaviour . As the analyses were not designed to simulate strain-rate effects, the predicted pile capacities should be compared with the end points of experimental holding stages as these correspond to negligible strain rates. Clearly, up to ground-level displacements of 0Á1D, there is a very good agreement between the FE simulations and the field results, with the loading capacity at v G ¼ 0Á1D being predicted with considerable accuracy for all pile geometries, thus suggesting that the predictive ability of the numerical model is practically unaffected by changes in pile slenderness ratio or diameter. The importance of adopting a small-strain stiffness model within the elastic region of the modified Cam Clay model is further reinforced in Fig. 10, where the initial stages of the loading tests and associated predictions are shown in greater detail, by focusing on ground-level displacements of up to v G ¼ 0Á01D. In addition to reproducing the behaviour at large displacements where the stiffness at moderate to large strains is dominant (Fig. 9), the FE model successfully captures the initial stiffness of the pile-soil system for all geometries. This excellent agreement between the FE model and the measured response demonstrates the suitability of the chosen modelling approach for the shear modulus at very small to moderate strains, which in turn provides a systematic connection with the results of the field and laboratory testing (Fig. 3).
To assess the quality with which the observed loaddisplacement response is reproduced by the FE model, the following accuracy metric is adopted where A ref is the area below the reference load-displacement curve (in the present case, that obtained in the field tests) and A diff is the area of the region delimited by the reference and simulated curves. The procedure followed for determining the reference curves is explained in the paper by Byrne et al.
. The geometric representations of A ref and A diff are shown in Fig. 11 for two distinct cases: in the first ( Fig. 11(a)), which is designated η 0Á1D , the overall performance of the developed numerical model is assessed by considering ground-level displacements up to 10% of the pile diameter, while in the second ( Fig. 11(b)), the range of displacements is limited to 0·1% of the pile diameter, to evaluate the ability of the numerical model to reproduce the observed initial stiffness of the pile. The latter metric is labelled η sd , as it focuses on the small-displacement range and its introduction is justified by the need to assess the adequacy of the numerical model to predict the response of the soil-pile system under operational conditions, which tend to be characterised by a large number of small loading cycles. A value of η ¼ 1 corresponds to a perfect reproduction of the reference curve, with any deviations from that behaviour leading to a reduction of its value (A diff is always positive and therefore η assumes values always equal to or below η ¼ 1. The interpreted values of η 0Á1D and η sd are illustrated in Figs 12(a) and 12(b), respectively, for the four considered pile geometries. The results of the FE analyses are seen to replicate with reasonably high accuracy both the overall observed pile response (average η 0Á1D . 88%) and the behaviour at small displacements (average η sd . 85%). It is also noted that the values of η 0Á1D are consistently high for all considered geometries, which include two distinct values of diameter and three values of the pile slenderness ratio, L=D, ranging from 3 to 10.

CONCLUSIONS
This paper presents a systematic and comprehensive procedure for developing a numerical model for the simulation of laterally loaded piles in overconsolidated glacial clays, which is applied in the study of laterally loaded piles installed at the PISA test site at Cowden. The 3D FE analyses were conducted before the pile tests were performed and the following were the key steps of the numerical procedure.
(a) The choice of the soil constitutive model and its calibration for Cowden till were based on the project objectives (monotonic loading, no cyclic or strain-rate considerations) and on available experimental evidence for the foundation soil. Particular emphasis was placed on the combined interpretation of field and laboratory testing results, as well as on the significant engineering judgement needed to ensure the consistency of calibration and realistic soil response under a wide range of strains. Based on experimental data, the adopted modified Cam Clay framework (Roscoe & Burland, 1968) was extended crucially with a non-linear Hvorslev surface (Tsiampousi et al., 2013), a generalised shape of the yield and plastic potential surfaces in the deviatoric plane (Van Eekelen, 1980), and a non-linear variation of the elastic shear modulus, dependent on both stress and strain levels (Taborda et al., 2016). precise prediction of the observed field behaviour in terms of load-displacement curves, pile deflections and bending moments. (c) The FE model developed in ICFEP (Potts & Zdravković, 1999) was validated by predicting accurately the response to lateral loading of four distinct PISA test piles installed at Cowden. Their geometries include two different diameters (0·76 m and 2·0 m) and three values of slenderness ratio, L=D (3·0, 5·25 and 10·0).
The close agreement of the predicted and measured responses of the four piles proved the premise of the PISA project that the new Winkler-type design method for laterally loaded piles could be derived directly from advanced, site-specific, FE modelling. As a result, the application of this FE model

Reference
Other v 0·001D Fig. 11. Graphical definition of the accuracy metric for assessing prediction methods for the (a) full load-displacement curve and (b) smalldisplacement range was key in the development of the new PISA monopile design methodology, which is outlined in the paper by Byrne et al. (2017).

ACKNOWLEDGEMENTS
The PISA 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, SSE, Vattenfall and Van Oord. The authors acknowledge gratefully the work of Socotec UK Ltd (formerly ESG) as the main contractor for the design and execution of the field testing programme.

APPENDIX 2. PLASTIC POTENTIAL ASSOCIATED WITH THE HVORSLEV SURFACE
The calibration process of the adopted constitutive model was concluded by assessing parameter β, which controls the shape of the plastic potential on the dry side. This was carried out using a trial-and-error procedure, where the simulated response for various values of β was compared to that measured experimentally. Although this was carried out for multiple tests, in Fig. 13 only the results obtained for test CR100IUC8.2an undrained triaxial compression test on a 100 mm sample collected at 8·2 m depth (see Table 3 and Zdravković et al. (2019) for test details)are shown in terms of the axial strain-deviatoric stress (ε a -q) and axial strainexcess pore pressure (ε a -Δu) diagrams. Based on the observed behaviour, a value of β = 0·20 was selected.
To demonstrate the model's performance using the final set of model parameters, listed in Table 1, two other triaxial tests carried out on 100 mm samples, collected at two depths -2·50 m and 5·00 mwere simulated. The initial conditions of the tests are listed in Table 3 (also in the paper by Zdravković et al. (2019)) and the results of the simulations, in terms of ε a -q and ε a -Δu plots, are compared with the test data in Fig. 14. Clearly, there is a good agreement between observed and simulated responses, confirming the suitability of the chosen modelling approach and respective calibration. NOTATION A diff area difference between reference and simulated curve in equation (13