A Numerical Model to Study the Response of Piles under Lateral Loading in Unsaturated Soils

: The interaction between a laterally loaded pile and the surrounding soil is typically limited to the shallower soil layer. Often, this zone is above the water table and therefore the interaction takes place under unsaturated conditions. The available evidence is scarce but suggests that unsaturated conditions play a major role on the pile’s response. The actual mechanisms governing the soil–pile interaction under unsaturated soil conditions are not understood entirely, and this paper provides a useful insight on this topic. The analysis is carried out with a fully coupled three-dimensional numerical model, the soil behaviour is simulated with a Modiﬁed Cam Clay Model extended to unsaturated conditions. The model accounts for the increase in stiffness and strength of unsaturated soils as well as the volumetric collapse upon wetting. The constitutive model is calibrated on the laboratory data and validated against centrifuge data with satisfying agreement. The results highlight the substantial differences in the soil reaction against the pile depending on different water saturation proﬁles. The study also shows that the inﬂuence of unsaturated conditions on the pile response increases as the pile’s ﬂexibility increases. Comparing the ﬁndings with currently available design methods such as the p-y curves, it is found that these do not adequately describe the unsaturated soil reaction against the pile, which opens the door for new research in the ﬁeld. The proposed numerical model is a promising tool to further investigate the mechanisms underlying the soil–pile interaction under unsaturated soils.


Introduction
Unsaturated soil mechanics is receiving rising interest from the academic community and the geotechnical construction industry. Some of the appeal derives from the potential use of unsaturated soil mechanics as a tool to reduce the carbon footprint of various civil engineering applications [1][2][3][4][5]. For instance, Speranza et al. [6] suggested that the increase in strength resulting from the suction s may be conveniently adopted to reduce the overdesign of the retaining structures in temporary works, leading to a reduction in emissions and materials used. Dynamic centrifuge tests on a shallow foundation performed by Borghei et al. [7] indicated that groundwater fluctuations significantly affect the seismic response of geotechnical systems and that neglecting them does not necessarily lead to safer or more economical design solutions. Several research studies investigated the effects of soil suction on the pile and piles group capacity under axial forces mainly by means of small-scale tests in cohesionless soils or numerical analyses [2,8,9] and, so far, the suction effect on the bearing capacity of a single pile seems quite well established.
The influence of the unsaturated soil conditions on the response of laterally loaded piles has received less attention; despite this, it is well known that the pile behaviour is largely affected by the mechanical properties of the shallower soil layer [7], which is often above the water table. To the best of the authors' knowledge, only few studies have so far addressed the relevance of the soil suction on the response of laterally loaded piles in unsaturated soils [10][11][12][13]. Stacul et al. [10] implemented a BEM model accounting for soil suction and highlighted appreciable differences on the pile's response even in the case of a shallow water table (1 m depth).
Lalicata et al. [11] carried out two centrifuge tests on a relatively short and rigid pile under lateral loading in an unsaturated silty soil. The soil models were compacted at two different void ratios: one loose (e 0 = 0.93) and one dense (e 0 = 0.75). The water content was nominally the same, w = 15%. Two additional tests were conducted under saturated soil conditions to provide a baseline response for comparison purposes. The study quantified the increase in the lateral stiffness and ultimate capacity due to the suction as well as the reduction in maximum bending moment at small load level (see Table 1). The effects of soil partial saturation on the pile response were found to depend on the soil state (loose or dense) and on the shape of the soil water retention curve, SWRC. Table 1. Main outcomes of the centrifuge tests by Lalicata et al. [11].

Soil State
Stiffness Increase * Lateral Capacity Increase *
Centrifuge modelling was revealed to be a powerful tool to investigate soil structure interaction problems even under unsaturated soil conditions. However, due to the experimental complexities, performing a complete parametric study investigating different water table elevations and/or different types of piles was not possible. Moreover, some data such as the lateral capacity in the unsaturated soil tests, the bending moments in one test, and the soil reactions were missing or uncertain.
In order to overcome the limitations of the experimentation, a numerical model has been developed. When validated against experimental data, numerical analyses represent a useful tool to better understand the geotechnical problems and to explore different scenarios from those experimentally investigated. This paper presents the main results of the numerical study of laterally loaded piles in unsaturated soils. The key objectives of the study are:

1.
Validate a numerical model against experimental data to investigate the behaviour of piles under lateral loading in unsaturated soils.

2.
Provide an insight on the experimental results by Lalicata et al. [11] with a focus on the flexural behaviour of the pile, on the soil reactions and on the strength mobilisation at higher loads.

3.
Present the capability of the model by extending the experimental findings to piles with different bending stiffness, which is known to be one of the key parameters of the pile response [15][16][17].
In the first section, the main characteristics of the centrifuge tests are briefly summarised. Next, the constitutive laws adopted for describing the hydro-mechanical soil behaviour, as well as their calibration against laboratory data, are presented. In the third section, the numerical model is presented and validated against the experimental data obtained in the centrifuge. The influence of soil partial saturation on the response of the pile is investigated, with an emphasis on the soil-pile reaction profiles. In conclusion, the last section focuses on the influence of soil suction on different types of piles.

Centrifuge Tests
Four centrifuge tests on a single free head pile were carried out at the IFSTTAR research centre of Nantes. The models were accelerated at 100× g. The soil used was a B-Grade kaolin (90% fine silt and 10% clay). Soil models were statically compacted on the dry side of the Proctor curve in a tub 300 mm in diameter. The final height was 18 cm. At the base of the container, a 10 mm thick sand layer, sandwiched between geotextile, acted as drainage boundary. The top of the model was sealed with a plastic film during the tests. The initial conditions of the centrifuge tests are summarised in Table 2. The model pile was a closed-end aluminium tube (diameter 12 mm, thickness 1 mm). The embedment length was 150 mm. The load was applied 35 mm above the ground surface. With the exception of T05, in the other tests the pile was instrumented with ten levels of strain gages. At the prototype scale (N = 100 g) the model reproduced a 1.2 m diameter 15 m long pile with a flexural rigidity of 3.9 × 10 6 kN·m 2 .
After the sample preparation, the key phases of the centrifuge tests were: 1.
Phase 0: Pile installation at 1× g. The bore was slightly larger (13 mm) than the diameter of the pile to avoid any damage to the strain gages due to the high firmness of the compacted soil samples. 2.
Phase 1: 1× g wetting. A zero-pore pressure was applied at the model bottom to reduce the post-compaction suction (Figure 1a). 3.
Phase 2: 100× g flight and equalisation. The model was accelerated at N·g. The water table level was controlled by a water reservoir supplied continuously from outside the centrifuge and connected to the base of the model (Figure 1b).

4.
Phase 3: Pile loading. The pile was pushed laterally at a slow and constant displacement rate of 0.003 mm/min to avoid the development of excess pore pressure ( Figure 1c). The main features of the tests are summarised in Table 2. Further details can be found in Lalicata et al. [11,18,19].

Constitutive Modelling
The hydraulic behaviour of B-grade kaolin is described using a modified version of the Gardner model [20]: where n controls the slope of the curve and −a/n is related to the air-entry value. The soil water retention curve, SWRC in (Equation (1)), is calibrated against the results of wetting-drying cycles performed using a suction-controlled oedometer cell [21] reported in [11]. Here, only the wetting branch is reported, as it is representative of the hydraulic paths imposed in the centrifuge tests. The experimental data, in terms of saturation degree versus suction, show a marked dependency on initial void ratio (Figure 2a), which is taken into account in the SWRC model (Equation (1)  The permeability of B-grade kaolin evolves with the saturation degree according to the following power law [22]: here, the saturated permeability k sat (depending on voids ratio) was determined by oedometer and falling-head permeability tests. It assumes the values of 4.5 × 10 −9 m/s for e 0 = 0.93 and 2.0 × 10 −9 m/s for e 0 = 0.75 [23]. The parameter α, set equal to 2.0, was calibrated against the experimental data of the infiltration phase, by means of preliminary sensitive analyses. The mechanical behaviour of the solid skeleton is described with a Modified Cam Clay Model (MCCM) extended to unsaturated conditions [24]. The MCCM well-reproduces the response of fine-grained soils under a broad range of initial conditions during wetting and loading processes, as demonstrated for single elements by Casini [25] and Casini et al. [26] as well as by other researchers [27][28][29] for boundary value problems. The model is characterised by a relatively simple mathematical formulation described in detail in Rotisciani et al. [30]. Here, only the main features are briefly recalled emphasising the main differences from the original model by Roscoe and Burland [31].
The model is formulated within the framework of critical state soil mechanics in terms of Bishop's effective stresses [32], defined as: σ ij = σ ij + χ s δ ij , with χ = S r [33]. The yield surface has an elliptical shape, symmetric with respect to the isotropic axis, whose evolution is controlled by a double-hardening mechanism: where . ε p v and . S r are respectively the increments of the plastic volumetric strain and saturation degree. In Equation (3), the first term coincides with the hardening law of the MCCM typically employed for saturated soils and the latter represents the hydro-mechanical coupling term introduced to capture some peculiar aspects of the hydro-mechanical behaviour of partially saturated soils (i.e., the changes in shear strength with suction and the volumetric collapse upon wetting). More specifically: (i) under drying paths, the elastic domain expands, leading to an increase in shear strength for heavy over-consolidated samples; (ii) under wetting processes, the yield surface shrinks keeping the current stress state in the elasto-plastic regime, and the solid skeleton collapses reducing its volume.
In the MCCM, the material response in the elastic domain is hypoelastic and the flow rule is associated. The values for constitutive parameters and their meaning are listed in Table 3. The calibration of the model was performed back-analysing a set of oedometer and triaxial tests carried out on saturated and unsaturated soil samples prepared with the same procedure adopted in the centrifuge tests and reported in [18]. Figure 2b shows the comparison between model predictions and measurements in oedometer tests performed on a fully saturated soil sample and under constant water content. The model captures all the main aspects of the soil response experimentally observed and provides predictions in agreement with experimental data in all the tests performed.
In the MCCM, the elastic shear stiffness, Equation (4), and the constant volume shear strength, Equation (5), are stress dependent: where p is the mean effective stress. As the latter is expressed in terms of Bishop's effective stress, the model accounts for the increase in stiffness and strength under unsaturated soil conditions.

Numerical Model
The centrifuge model is related to the full-scale prototype following the scaling laws [34]. A centrifuge model accelerated of N times the Earth's gravity (g), replicates a full-scale prototype N times larger. The relevant scaling laws for this study are summarised in Table 4. Unless otherwise stated, the simulations are performed at the prototype scale. The pile response is studied by means of three-dimensional coupled analyses via the finite element code Abaqus/Standard. The constitutive soil model is implemented in a user-defined subroutine [35,36]. The hydro-mechanical behaviour of the kaolin is described with the Modified Cam Clay Model extended to unsaturated conditions and the modified Gardner SWRC previously introduced.

Three-Dimensional Model
The 3D model simulates half of the cross section of the pile and the soil model taking advantage of the symmetry along the vertical plane containing the force, Figure 3. The finite element volume, 18 m thick and 30 m wide, is discretised with three-dimensional solid continuum porous elements C3D8P for both the soil and the pile. A mesh refinement is adopted near the pile. The latter (embedded length L = 15 m, diameter D = 1.2 m) is installed in the centre of the soil layer, after the initial step in which the stress state from the axisymmetric analyses is applied. Quadratic beam elements B32 with null stiffness (Young Modulus E = 1 kPa) are constrained to the central nodes of the pile, in order to efficiently obtain the bending moment profile along the pile length. The model restricts the vertical displacements at the base. At the lateral boundaries, the displacements are constrained in the normal direction. The pile's head is free to rotate as in the centrifuge tests. The rotation along the X-axis of all nodes belonging to the structural elements on the symmetry plane is restrained. The water table position is fixed by imposing the corresponding pore pressure value at the bottom of the model. The load is applied at a master node, rigidly connected to the nodes of the pile's head. A concentrated moment is applied to the same node to take into account the eccentricity of load with respect to the ground level. As in the centrifuge tests, the loading rate is slow enough to ensure a drained response of the soil.
The pile is modelled as a linear elastic material with a Young's modulus E p = 38.3 GPa and a Poisson's ratio υ = 0.2. The model pile in centrifuge tests was a hollow aluminium tube. For simplicity, the pile is modelled as a solid element. The same flexural rigidity, E p I p = 3.9 × 10 6 kN·m 2 (where I p is the moment of the inertia of the pile), between the simulation and the tests, is guaranteed thanks to the adoption of an equivalent Young' modulus E p rather than that of the aluminium.
The soil-pile interface is described with a purely frictional Mohr-Coulomb law assuming a friction angle at the interface δ = ϕ /2, where ϕ is the critical state friction angle.

Stress History of the Model
The mechanical response of fine soils depends on the history of the material. In the centrifuge tests, prior to the loading of the pile, the soil model is subjected to mechanical and hydraulic processes. These influence the final state of the soil and, consequently, the response of the pile upon loading.
Taking advantage of the symmetry of the model and of the perturbations applied, an axisymmetric analysis has been set up to reproduce the soil stress history prior to the loading of the pile. The details of this procedure are presented in Rotisciani et al. [37]. Here, only the main aspects are reported. The results are at the model dimensions.
In the centrifuge tests, the pre-hole was slightly larger (13 mm) than the pile's diameter (12 mm) to avoid damaging the strain gages. This small tolerance of 0.5 mm was taken into account in the axisymmetric model. When the clearance between the nodes of the pile and the nodes of the soil vanishes, the nodes interact with the same contact law used in the three-dimensional model.
The initial conditions follow from the measurements reported in Table 5. Owning the small dimensions of the model, all state variables are considered uniformly distributed within the soil layer. The initial effective stresses are assumed to be isotropic and equal to S r ·s; consequently, the soil results in a slight over-consolidated state [38,39]. The wetting process is simulated by imposing zero pore pressure at the base of the model (Phase 1). After a time period which varies depending on the specific test under consideration, the gravity loading exerted on the solid skeleton is increased 100 times (Phase 2). At the same time, the pore pressure at the base is increased to the desired value (190 kPa in T05, T09; 120 kPa in T06, T08) corresponding, at equilibrium, to a water table coinciding to the soil surface or located at 7 cm of depth.
The simulations are compared with the gathered data (surface movements, pore pressures, and volume water absorbed) with satisfying results. For the sake of brevity, the comparisons will not be shown here but they are thoroughly discussed elsewhere [19,37]. Figure 4 shows the profiles of the stresses and state variables during the analysis. The results refer to the test T06 (looser unsaturated test, LU). The key findings are: 1.
During the wetting event at 1× g, the soil swells due to the effective stress reduction and the increase in degree of saturation. Although the elastic swelling prevails, this is an elasto-plastic process because the contraction of the yielding surface following the S r increase is faster than the effective stress reduction.

2.
The fast g-increase, and thus the increase in the vertical total stress, entails a positive variation of the pore pressures. In the upper part of the model (z < 10.5 cm), the soil remains unsaturated, the vertical effective stress increases and the void ratio reduces.
In the lower part (z > 10.5 cm) the combination of high values of degree of saturation and the high vertical stress leads to instantaneous saturation of the soil and thus to the developments of positive excess pore pressures.

3.
In the equalisation following the g-increase, the water table gradually moves to 7 cm of depth. The dissipation of the excess pore pressures is coupled with the development of positive volumetric strains in the lower part of the model. Volumetric collapse is observed in the unsaturated zone between 3 and 10.5 cm of depth.

4.
The void between the pile and the soil represents a free boundary for the soil as long as the void exists. During the analysis, the soil surrounding the pile deforms not only vertically, but also laterally, until closing the void and restoring, at same point, the oedometer conditions. 5.
At the end of the equalisation in flight, the gap was completely closed for all the analyses. Therefore, this was not considered in the three-dimensional model.

Initial Conditions for the Three-Dimensional Model
The profiles (with depth) of the main variables at the end of the axisymmetric analyses, and thus the initial state of the three-dimensional model, are presented in Figure 5: The vertical effective stress profile varies with the elevation of the water table. The small difference in the two unsaturated cases depends on the slightly different degree of saturation with voids ratio, as in the SWRC.

2.
Although the current stress state belongs to the yield surface in all the tests (i.e., the soil is in an elasto-plastic state), the ratio between the horizontal and the vertical effective stress, i.e., the K 0 values, depends on the initial void ratios. 3.
In the loose soil cases (LS and LU) K 0 is, practically everywhere, lower than one and close to the normal consolidated value predicted by the model [40]. The soil is then close to a normal consolidated state.

4.
For the dense soil cases, K 0 reduces with depth but it is higher than one for the firsts 6 m and 10 m of depth in DU and DS, respectively. The soil is in extension conditions, and thus, during the pile's loading, the stress path initially moves inside the yielding surface. The soil response is initially elastic and similar to that of an over consolidated soil, at least in the region involved in the soil-pile interaction.

5.
The void ratio profiles give rise of the different state of the material. In LS and LU, e reduces with depth while in DU and DS is practically constant with it.

Model Validation
The comparison of the load-deflection relationships between experimental data (CE) and numerical predictions (FEM) is shown at the prototype scale in Figure 6a,b for the loose and dense soil respectively. The model predictions are satisfying for the looser state over the entire range of displacements (Figure 6a), while the model tends to overestimate the initial stiffness and underestimate the soil-yielding threshold load for the denser layer (Figure 6b). From a qualitative point of view, the model correctly reproduces the increase in both stiffness and strength of the soil-pile system due to the presence of an unsaturated zone above the water table. As experimentally observed, the model points out that the effects of partial saturation reduce as void ratio decreases (i.e., with the increasing in stiffness of the saturated soil).   The model fairly reproduces the bending moment profiles at the different loading stages under unsaturated and saturated conditions with more satisfactory results for the normal consolidated case, Figure 7a. At 0.3 MN of load, the model is more rigid of the data and the prediction of the shape of the bending moment curve is unsatisfactory, even though the value of the maximum moment is correct. The predictions improve as the external load increases. The maximum bending moment and its position well agree with the experimental data in this case.
Because of the overestimation of soil-pile initial stiffness for the dense soil (Figure 6b), the model slightly underestimates the value of the maximum bending moment and its position along the pile depth, Figure 7b,c.

Detailed Analysis of the NC Model
As may be expected, the model provides a better agreement with data in the loose soil cases, where the soil is essentially in a normal consolidated state. This is also the case where the effects of partial saturation are magnified. Consequently, a detailed analysis of the pile response is presented for these cases only (LS and LU) in the following section.

Global Pile's Response
The pile response is further analysed by means of normalised profiles of deflection y/D, bending moment M/Hē (whereē is the eccentricity of H from the ground level), and the soil reaction p/p max . p max is the maximum positive reaction calculated at y/D = 25% and identified as the failure condition. The comparison between the saturated and unsaturated conditions is presented in Figure 8 for the load of 0.4 MN. While the deflection and the bending moment are direct output of the analysis, the soil reaction is obtained through the double derivation of the interpolating function of the bending moment. The accuracy of the result strongly depends on the class of the interpolating function and the number of points available for the calculation [41]. In this study, a 5-intic spline is used to interpolate the bending moment profile over a 40-node beam. The validity of the calculation was verified by comparing the latter with the one resulting from the integration of the forces acting on the soil-pile interface for several loads, founding negligible differences (<2%) only for high loads. Moving from saturated to unsaturated conditions, the maximum bending moment reduces of approximately 32% and its position moves upward from 6 to 3.2 m from the ground level. This indicates a strong reduction in the significant soil volume involved in the soil-pile interaction. The soil reaction against the pile is very different in the two analysed scenarios. Indeed, under saturated conditions, the soil reaction assumes a null value at the ground level and it increases in depth up to the position of the maximum bending moment. Below this depth, it starts to decrease, changing sign at the rotation point at 10 m. Similar to the moment profile, the soil reaction involves the whole length of the pile, as expected for a very stiff pile in soft soils. On the contrary, in partial saturation regimes p/p max is maximum at the ground level and then plummets to zero moving downward, as a consequence of effective stress states far from being null at the ground surface. The deflected shape of the pile is affected by the position of the water table. When the shallower layer is in unsaturated conditions, the horizontal displacements at the pile tip are small and the deflection occurs in the upper part of the pile. Conversely, in fully saturated conditions, the rigid rotation (at a depth of approximately 11.5 m) must be added to the lateral deflection giving values of displacement greater than those obtained under unsaturated conditions. The partial saturation increases the yielding threshold of the material; thus, at the same external load, plastic strains and mobilised shear strength are smaller compared to those obtained under saturated conditions. A representative case is reported in Figure 9, in terms of isocontours of the mobilised strength q/q max (q is the deviatoric stress and q max is the maximum deviatoric stress attainable at the current p ) with depth z and distance x/D for a horizontal load H = 0.4 MN. The mobilised strength in the upper part of the soil layer under unsaturated conditions is negligible, whereas, under saturated conditions, a significant concentration of q/q max is depicted near the ground surface that propagates towards the rotation point.
The pile response at large displacements (y/D~0.25) is analysed in Figure 10. In both saturated and unsaturated conditions, the pile undergoes almost a rigid body motion. Indeed, the change in the strength distribution into the soil due to partial saturation is not able to modify the failure mechanisms at elevated horizontal displacements, as demonstrated by the contours of the mobilised strength in Figure 11. It is worth noting that the external force is very different in the two cases; indeed, under unsaturated soil conditions, this is approximately three times higher than that obtained in the saturated case. It is imperative to properly take into account the significant differences in the soil reaction. Under unsaturated conditions, p/p max is almost constant and equal to the maximum value up to z = 6 m and then reduces with depth, identifying a null value at a depth of~11 m.
For z w = 0 the reaction is similar to that observed in Figure 8 with H = 0.4 MN, which was very close to failure, as shown in Figure 6.

p-y Curves>
The numerical p-y curves at different depths (z = 2, 4, and 6 m) are presented in Figure 12 for the saturated and unsaturated soil conditions. For any given displacement, the soil reaction is higher under unsaturated soil conditions than that obtained for z w = 0. When z w = 7 m, the curves yield for higher displacement and loading values compared to the saturated case, at any depth z. In the displacement's range investigated, the soil reaction keeps increasing in the unsaturated soil while it is close to the ultimate values in the other case.

Comparison with Literature Solutions
Historically, p-y curves represent the most used method to compute the pile response under lateral loading. Assuming a Winkler formulation for the soil, several authors provided the definition of the subgrade modulus K, the ultimate soil resistance p u , or the shape of the p-y relationships to be used depending on the characteristics of the soil surrounding the pile [42][43][44][45][46][47]. In this study, a simple hyperbola was used to fit the p-y curves as suggested by several researchers [48,49]. Employing this procedure, K and p u are automatically provided for every elevation of the pile.
Results are presented in Figures 13a,b and 14a,b as function of the water table position, together with the solutions usually adopted in the literature. As the loading was a drained event, the comparison is made with the solutions formulated in terms of effective stresses parameters.  As may be expected, for the saturated soil the distributions of K and p u linearly increase with depth starting from null values at the ground surface. The unsaturated soil exhibits significant stiffness and strength values even at z = 0, which slightly increase with depth. Differences between the two trends tend to reduce with depth as suction decreases. In terms of stiffness, the results for z w = 0 lie within the range suggested by Terzaghi [42], Figure 13a. In the unsaturated case, z w = 7 m in Figure 13b, the experimental results fall in the range identified by ibid for all the depth presented although in the shallower layer, z ≤ 2 m, the fitting trend is different. However, it is worth noticing that instead the stiffness distribution is very different from the one provided by ibid. Terzaghi [42] adopted two linear trends increasing with the total stress applied. The result obtained within this work shows a small increase in stiffness with depth because the total stresses increment is balanced by the progressive increase in the degree of saturation with z, approaching the water table.
For z w = 7 m (Figure 14b), adopting the classical assumption of dry soil above the water table (unit weight γ = 18 kN/m 3 ), the literature solutions fail to predict the limiting pressure both in qualitative and quantitative terms. In particular, they underestimate the soil resistance in the shallower depth and overestimate it approaching the position of the water table. Adopting the Bishop's effective stress, the capability of the literature solutions improves even in unsaturated conditions; although, the gradient is still too high.

Influence of Partial Saturation on Different Type of Piles
The numerical model is able to capture the main effects of the soil partial saturation on the pile's response experimentally observed. The comparison is particularly satisfying for normal consolidated soils. The centrifuge tests were, however, limited to one type of pile and one position of the water table. With the confidence gained in the numerical model, and without the ambition to present a complete parametric study which is beyond the scope of this work, an additional set of analyses is carried out to extend the discussion on the influence of partial saturation on the pile's response.
The geometry, the soil properties and state (loose soil: e 0 = 0.93) as well as the hydraulic conditions are the same of the centrifuge tests described above. The flexural rigidity of the pile, E p I p , is varied from 10 4 to 10 7 kN·m 2 . This very wide interval covers the most common types of piles used in practice from the tubular steel piles to the concrete solid piles [52,53], and includes also the pile tested in the centrifuge by Lalicata et al. [12] that had a flexural rigidity of 3.9 × 10 6 kN·m 2 .
The response of the pile is summarised using the following parameters: 1.
H/y, the stiffness of the load-deflection curve; 2. M max , the maximum bending moment M max ; 3.
z max , the position of M max ; 4.
L c , the critical length.
All these parameters are evaluated at a small displacement level, such as y/D = 1%, which may represent the working load conditions, Price and Wardle [54].
The secant stiffness calculated for the saturated and unsaturated soil conditions is plotted in Figure 15a against the flexural rigidity of the pile. In both cases, the stiffness increases with E p I p and it is always higher under unsaturated soil conditions than under saturated conditions. For z w = 0 m, the soil-pile stiffness is progressively less affected by the E p I p variations when this is higher than 1 × 10 6 kN·m 2 , denoting the tendency to a rigid behaviour of the pile; this threshold may move forward to E p I p > 6 × 10 6 kN·m 2 for z w = 7 m. The increase in secant stiffness (ISS) due to partial saturation may conveniently be expressed by: where the subscripts "unsat" and "sat" are for unsaturated and saturated soil conditions, respectively. ISS values, presented in Figure 15b, clearly highlights how the effects of partial saturation increase for decreasing values of E p I p , i.e., long piles. When E p I p > 6 × 10 6 kN·m 2 and the pile behaviour switches from long to short, the increase in the secant stiffness is minimum, although still remarkable (+200%), and it is independent from E p I p . It is worth noting that the pile tested in the centrifuge by Lalicata et al. (2019) represents a lower bound of the beneficial effects of the partial saturation on the pile response. The normalised maximum bending moment is reported in Figure 16 in function of E p I p . The trend is somehow similar to that of the secant stiffness in Figure 15: the maximum bending moment increases with E p I p and the gradient reduces when the pile behaves as a short pile. For z w = 7 m the maximum bending moment is always smaller than that measured for z w = 0. The reduction in the maximum moment due to partial saturation, defined likewise as the increase in the secant stiffness, is practically constant and equal to 20% in the explored E p I p range. The critical length L c is the portion of the pile that reacts to the perturbation applied at the pile head, thus identifying the significant volume of soil involved by the pile deflection. Randolph (1981) defines L c as the depth where the bending moment is zero. Consequently, a flexible (or infinitively long) pile is a pile that has a total length, or more correctly a slenderness ratio L/D, higher than the critical slenderness of the pile L c /D. Using L c following Randolph's definition, Figure 17 shows that under unsaturated soil conditions the critical slenderness ratio is always lower than that observed for z w = 0 and that the pile behaves as a flexible pile for a wider range of E p I p . The normalised position of the maximum bending moment, z max /D, varies accordingly with M max and, under unsaturated soil conditions it moves toward the ground surface of 1-2 diameters.

Concluding Remarks
A three-dimensional numerical model is developed to study the behaviour of piles under lateral loading in unsaturated soil conditions. The soil stress-strain response is described by the modified Cam Clay model extended to unsaturated conditions and the water retention properties are represented with the modified Gardner's model. Both models are calibrated against laboratory data in saturated and unsaturated conditions at different water contents and void ratio.
The model is validated against the centrifuge tests results of Lalicata et al. [12]. The comparison between model predictions and experimental data is satisfactory for normal consolidated soils under saturated and unsaturated conditions. With the confidence gained by the model validation, the following findings may be highlighted with reference to normal consolidated states: 1.
The model correctly captures the increase in stiffness of the pile response due to unsaturated conditions. 2.
The value of the maximum bending moment is in agreement with the data. At small load levels, the computed bending moment profile denotes a slightly stiffer response than that experimentally measured. The prediction improves as the external load increases.

3.
The magnitude and the shape of the soil reaction profiles against the pile are substantially influenced by the unsaturated zone. 4.
The parameters of the p-y curves, namely, the initial stiffness and the ultimate capacity, depend on the suction level of the soil. Under unsaturated soil conditions, those parameters start from high values at the ground level and then slightly increase with depth.
Besides the increased stiffness and strength of the unsaturated soil, in this case, the improved performance of the pile also depends on the plastic strains and mobilised strength levels. For any external load applied, these factors are always smaller under unsaturated conditions.
The impact of partial saturation on the pile response increases for flexible long piles. When the pile is short and almost rigid, the increment in the secant stiffness of the load deflection curves is minimum but still notable (+200%).

2.
The reduction in the maximum bending moment (20%), due to unsaturated conditions, is less affected by the flexural rigidity of the piles. 3.
The unsaturated shallower layer reduces the critical length and the depth of the maximum bending moment.
This study demonstrated the validity of the numerical model in capturing and detailing the response of laterally loaded piles in unsaturated soils. This model can be then adopted to carried out a more complete parametric study varying, for instance, the water table depth, the water retention properties of the soil, the pile head' fixity and its stiffness.
The study also demonstrated that the methods currently used to design the piles (i.e., p-y curves) turn out to be inadequate when trying to describe the unsaturated soil reaction against the pile. Consequently, additional experimental and numerical research is needed to fill this gap in the actual knowledge. Funding: This research received no external funding.

Data Availability Statement:
The data that support the findings of this study are available from the corresponding author upon reasonable request.
Acknowledgments: This research was performed in the framework of the GEOTRANSALP-PILE-UNSAT agreement in which the authors are involved. The institutions joining the agreement are greatly acknowledged.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations b
hydro-mechanical coupling parameter of the MCCM extended to unsaturated soil conditions D outer diameter of the pile E p Young' modulus of the pile E p I p flexural rigidity of the pilē e eccentricity of the applied load on the pile e, e 0 void ratio, initial void ratio g gravity acceleration H lateral load I p moment of the inertia of the pile K subgrade modulus K 0 horizontal to vertical stress ratio k, k sat soil permeability, saturated permeability L embedded pile length M, M max bending moment, maximum bending moment M CSL Slope of the critical state line in (q, p ) plane N scaling factor in the geotechnical centrifuge N 0 e at p = 1 kPa on the saturated NCL n, a Gardener's parameters p, p max , p u soil reaction, maximum soil reaction, soil resistance p mean effective stress p c hardening parameter . p c incremental variation of the hardening parameter q, q max deviatoric stress, maximum deviatoric stress attainable at the current p S r , S r0 degree of saturation, initial degree of saturation . S r incremental variation of the degree of saturation s suction t time u pore pressure w, w 0 gravimetric water content, initial water content x horizontal distance from the centre line of the pile y lateral deflection of the pile z, z w depth, depth of the water table