Numerical modelling of the response of an unsaturated silty soil under wetting and gravitational loading processes

This paper presents the results of a numerical study aiming at simulating the response of an unsaturated fine-grained soil under wetting and gravitational loading processes. This study is based on the results of some centrifuge tests carried out to assess the influence of partial saturation on the laterally loaded pile response. The hydro-mechanical behaviour of the silty soil is described using a constitutive model adapted to unsaturated conditions. The model predictions are compared with the measurements provided by LVDTs and laser transducers in the first phases of the experimental study. Besides validating the model, the numerical study aimed at investigating the influence of the after-compaction conditions on both the displacement field and the evolution of the more significant state variables during imbibition and gravitational loading processes. Finally, an additional analysis is conducted to determine the effects of the pile installation on the soil response.


Introduction
In geotechnical engineering practice, pile foundations are designed assuming fully saturated or dry soil conditions. In many situations, however, natural soils are found in a state of partial saturation as the groundwater table is at a great depth. The experimental study conducted by Lalicata et al. [1] aimed at assessing the influence of partial saturation on the response of laterally loaded piles installed in fine-grained soils. The centrifuge tests were performed on loose and dense soils under fully and partially saturated conditions. Each test consists of three or four phases: imbibition at 1xg, inflight phase, horizontal loading step, and, eventually, rising of the water table.
The small-scale tests were numerically simulated to deeper understand the contribution of suction to the pile response under working loading [2]. The numerical results presented in this paper refer to the looser unsaturated soil layer and aim at investigating the soil response during the imbibition and in-flight phases.
The hydro-mechanical behaviour of the silty soil was described using the Modified Cam Clay Model (MCCM) extended to unsaturation conditions [3] with a modified version of the Gardner water retention curve [4]. A series of numerical analyses were performed to: • validate the model in a real case study; • determine the impact of the after-compaction conditions on both the displacement field and the evolution of the more significant hydro-mechanical variables; • assess the influence of the procedure employed for pile installation on the soil response.
The results of this study can be very helpful in all cases, where small uncertainties affect the evaluation of the after-compaction conditions. The knowledge of the effects induced by any change in the initial state variables is crucial to identify the set of initial values providing the best fitting of the measurements.
The paper is organized as follows. In the first part, the main characteristics of the centrifuge tests are briefly reminded. In the second part, the numerical model is described emphasizing the constitutive laws adopted for describing the hydro-mechanical behaviour of the silty soil. Then, the numerical results are compared with the measurements to verify the predictive capabilities of the model. Finally, the influence of the after-compaction conditions and the installation procedure on the soil response is investigated.

Experimental study
The experimental study was conducted by Lalicata et al. [1] using the geotechnical centrifuge of IFSTTAR in Nantes. The material used is the commercial B-Grade kaolin consisting of 90% silt and 10% clay. A detailed description of the centrifuge tests is reported in [5], herein only the main features are briefly reminded.
The B-Grade kaolin was statically compacted in six layers under a constant displacement rate, reaching the final height of 180 mm. Below it, a sand layer 10 mm thick surrounded by a geotextile was deposited.
After soil compaction, an instrumented model pile of diameter 12 mm and embedment length 150 mm was installed into the silty layer. Before installing the pile, a pre-hole of diameter around 13 mm was drilled using a manual screwing system. The pre-hole avoids any damages to the extensometers housed on the external lateral surface of the pile.
The experimental study consists of 4 tests conducted on loose and dense soils under fully and partially saturated conditions. The soil models were prepared at void ratios e 0 of 0.93 and 0.75, and saturation degree S r0 of 0.42 and 0.5, for loose and dense states, respectively. These values are assumed as reference despite small variations in the after-compaction conditions. Indeed, the adopted compaction procedure does not permit to know the real spatial distribution of void ratio, suction and saturation degree within the soil layer. Each test is made up of four phases. After soil compaction, the strong-box is hydraulically connected to a water tank and null pore water pressure is prescribed at the base of the soil model. The imbibition process at 1xg reduces the aftercompaction suction leading the negative pore pressures into typical ranges for practical applications. Then, the gravitational acceleration is increased 100 times while the pore water pressure at the base is raised to 120 or 190 kPa, corresponding to a water table located at 70 mm from the ground surface or coincident with it. Once attained stationary conditions, the pile is laterally loaded and, finally, the soil layer is fully saturated. The model was extensively instrumented to follow the soil-pile response during equalisation and pile loading phases. In this study, particular attention is focused on the measurements given by the LVDTs and laser transducer monitoring the soil vertical displacements and the water level in the reservoir, respectively.

Numerical modelling
The numerical study presented in this paper aims at investigating the soil response during the imbibition and in-flight phases. All small-scale tests were numerically simulated, achieving satisfactory outcomes. Herein, only the results related to the centrifuge test performed on the loose unsaturated soil layer are analysed.
The test was numerically simulated using the commercial finite element code Abaqus/Standard at the model scale. Since both the water flow and soil strains develop along the vertical direction, an axial-symmetric model was adopted. The soil column is 190 mm high and consists of 78 8-node axial-symmetric elements CAX8P.
The mechanical behaviour of the B-Grade kaolin was numerically simulated using the Modified Cam Clay model extended to unsaturated conditions [3]. This model is formulated in terms of Bishop's effective stresses and well-reproduces the response of fine grained soils, as demonstrated by several authors, among others: [6], [7], [8], [9]. The model is fully formulated using few constitutive parameters that can be easily determined with laboratory tests. They are: the slope in the semi-log compressibility plane of the unloading-reloading line κ; Poisson ratio ν, the slope in the stress invariants plane of the critical state line M , the specific volume at p'=1 kPa on the normal consolidation line in the semi-log compressibility plane N and, finally, its slope λ. The model is calibrated against oedometric tests carried out on compacted soil samples keeping constant water content and suction ( [2]).
The hydraulic behaviour of the B-Grade kaolin is described adopting a modified version of the Gardner model ( [4]) as Soil Water Retention Curve (SWRC): where s is the suction, n controls the slope of the curve, and a is related to the air-entry value. The SWRC is calibrated against the results of wetting-drying cycles performed under oedometric conditions (see Fig. 1). The permeability k of the B-Grade kaolin varies with S r according to the power law ( [10]): where k sat is the saturated permeability and β controls the shape of the curve. The values for all constitutive parameters are reported in Table 1.
The initial conditions are considered uniformly distributed over all the soil column accounting for the small dimensions of the model. The initial effective stresses are assumed isotropic and equal to suction scaled by saturation degree. This assumption is rather reasonable since the suction is significantly greater than the gravity-induced stresses. In the MCCM, the preconsolidation pressure p c ' is obtained once known the initial conditions in terms of effective stresses, void ratio, and saturation degree. In the compressibility plane, p c ' corresponds to the abscissa of the intersection point between the unloading reloading line passing through the initial stress state and the normal consolidation line related to S r0 [11]. As expected, the soil results overconsolidated.
The imbibition phase is modelled by imposing zero pore water pressure at the base of the soil domain. After 860 min the gravity loading is built-up 100 times, and the pore pressure at the base is increased up to 120 kPa.
Accounting for the uncertainties in the aftercompaction conditions, five simulations were performed starting from the initial conditions summarized in Table  2. In T6, a series of gap elements was introduced into the model to simulate the pre-hole. This analysis was performed to verify to what extent the procedure employed for pile installation affects the soil response.

Reference analysis
In this section, the results of T1 are analysed and compared to the experimental data; this test is the one providing the best fitting prediction of the soil response. Figure 2a shows the temporal evolution of both the vertical displacements at the soil surface and the absorbed water volume during the imbibition phase. The measurements provided by LVDTs lie within the region coloured by pale grey. The wetting process induces the development of vertical displacements directed upwards attaining final values in the range (2.8, 3.4) mm. The water infiltration into the soil layer entails an increase in water volume by around 2 l. Globally, the model wellreproduces the soil response even if the initial displacement rate is slightly overestimated whereas the final one is underestimated. Moreover, the model predicts absorbed water volumes greater than the measured ones with maximum differences of around 25%. Figure 2b shows the temporal evolution of the vertical displacements at the upper boundary during the in-flight phase. The increase in the gravity loading leads to immediate vertical displacements of about 7 mm increasing in time up to 10 mm. Also in this case, the model predictions match well the measurements given by LVDTs. Figure 3 shows the evolution over time of the profiles with depth of the pore water pressures u, vertical σ v ', radial σ r ', circumferential σ θ ' effective stresses, and specific volume v in the most significant moments of the test. Obviously, under oedometric conditions, radial and circumferential effective stresses coincide.  The imbibition process entails a significant reduction in the after-compaction suction that reaches final values ranging from zero (at the base) to 80 kPa (at the top). It is important to note that this phase ends before attaining stationary conditions, corresponding in this case to pore water pressures close to zero in the entire soil layer (Fig.  3a, line 2). No significant change in total vertical stresses is observed in the imbibition phase. Therefore, σ v ' evolves coherently to the variations in suction scaled by saturation degree (Fig. 3b, line 2) and horizontal stresses σ r ' and σ θ ' vary according to the imposed oedometric conditions (Fig. 3c, line 2). The wetting-induced stress unloading induces positive volume changes, especially in the lower part of the silty layer where the most significant variations in u occur (Fig. 3d, line 2).
In the in-flight phase, the pore water pressure at the base is fixed to 120 kPa. The g-level increase provokes the sudden attainment of full saturation in the lower part of the soil column, where positive pore water pressure excesses ∆u develop (Fig 3a, line 3). In the upper part, the material remains in unsaturated regime and the u profile modifies increasing its slope. Upon ending the consolidation-equalization process, the u profile is hydrostatic with a water table located 7 cm below the soil surface (Fig 3a, line 4). The increase in the gravity loading induces positive variations in all effective stress components and, consequently, reductions in volume of the solid skeleton (Figs 3b, c, d, line 4). In the upper part, the soil compression occurs instantaneously with the gravity loading since soil conditions are substantially drained. This explains the occurrence of the immediate settlements shown in Fig. 2b. In the lower part, the volume changes develop progressively over time because they are associated to the dissipation of ∆u (Fig.   2b). Finally, it is worth to notice that in the upper part the soil response is controlled by elastic stiffness, because the small increase in effective stresses is not sufficient to lead the stress state on the current yield surface. Contrary, the volume changes in the lower part of the soil domain (z>3 cm) are mostly irreversible.

Influence of after-compaction conditions
As mentioned above, the numerical parametric study consists of five analyses conducted under oedometric conditions by modifying the initial values for void ratio and saturation degree (see Table 2). For each variable, the selected range includes the initial reference values assumed in the experimental work. Figure 4a shows the variations in time of both the vertical displacements at the soil surface and water content during the imbibition phase in T1, T2, and T3. These tests differ only for the initial void ratios ranging from 0.91 to 0.95. Keeping unchanged both the effective stresses and saturation degree, the increase in e 0 causes a reduction in p c '. In the compressibility plane, the current state moves upwards approaching the normal consolidation line related to S r0 . Consequently, the ratio of p c ' to mean effective stress p' decreases from 1.7 (T1) to 1.0 (T3).
In all tests, the soil swells upon wetting, but the magnitude of the vertical displacements decreases for increasing initial void ratios. This result is primarily due to the differences in the position of the stress state with respect to the yield surface and, secondarily, to the variations in elastic stiffness modulus. In wetting processes, the stress unloading is always accompanied by the shrinkage of the elastic domain. When the yield surface captures the current stress state, positive plastic volumetric strains start to develop. These deformations partially balance the elastic swelling induced by the stress unloading. Thus, it follows that the vertical displacements are to be lower in magnitude in T3, where the soil is already in the elastic-plastic regime.
The evolution of the absorbed water volume is not affected by the initial void ratio (Fig. 4a). This result is not surprising since both the SWRC and permeability curve are independent of the current void ratio, and the initial conditions in terms of S r0 and s 0 are the same in all three tests. Figure 4b shows the evolution over time of the vertical displacements at the soil surface in T1, T2, and T3 during the in-flight phase. In all cases, the soil settles as soon as the g-level is increased. Since the S r0 profiles coincide in the all three tests, the small differences in the short-term settlements are to be attributed to the different soil compressibility. The latter affects also the final stationary values that increase passing from T1 to T3 according to the increase in void ratio.
As revealed by Table 2, T1, T4, and T5 differ for the initial saturation degree and thus the suction, whereas the void ratio remains unchanged. The reduction in S r0 entails an increase in both effective stresses and preconsolidation pressures, which vary keeping quite constant the ratio p c '/p'.   Figures 5a and b show the numerical results of the three tests during the imbibition and in-flight phases, respectively. When zero pore water pressure is imposed at the base, the gradient of the hydraulic head increases passing from T1 to T5. Consequently, the absorbed water volume grows up more quickly and attains higher final values in T5 than in the other tests (Fig. 5a). Because the u w distributions is mostly the same at the end of the imbibition phase, the higher values of the vertical displacements at the top are related to the initial higher value of suction. Finally, the soil response during the in-flight phase is poorly affected by the initial conditions in terms of S r0 (Fig. 5b).

Influence of pile installation procedure
Before installing the pile, a pre-hole 150 mm long and 13 mm wide was drilled using a manual screwing system. The axial-symmetric numerical analysis T6 was performed to verify to what extent the soil-pile response was influenced by the installation procedure.
The geometry of the finite element mesh was modified by increasing its width up to 150 mm, corresponding to the strong-box radius (Fig. 6). The whole soil volume was discretized by means of CAX8P elements, whose dimensions were kept unaltered with respect to the previous simulations. The lateral boundaries prevent any horizontal displacements as expected along the symmetry axis (left boundary) and the contact surface between the soil and the strong-box (right boundary). The nodes located at the base are constrained in both directions.
GAPCYL elements were introduced to simulate the pre-hole. These elements allow the nodes to be in contact  More in details, the nodes start to interact only when the separation distance vanishes and the gap contact elements are closed. In this case, the separation distance was set equal to 0.5 mm, coinciding with the difference in the radii of the pre-hole and pile. These elements were installed up to a depth of 150 mm, corresponding to the pile length, at a distance of 6 mm from the left boundary coincident to the pile radius.
The constitutive laws employed for describing the hydro-mechanical behaviour of silty and sandy soils remain unchanged. The numerical analysis consists of two phases again. The initial conditions coincide with those assumed in T1 since it provides the best fitting prediction of the experimentally observed soil response.
No significant differences were found in the displacements induced by the wetting process and the glevel increase; they were substantially coincident to those observed in T1. Contrary, the numerical results pointed out a remarkable effect of the installation procedure on the stress distribution over the entire layer. Figure 7 shows the profiles of u, σ v ', σ r ' and σ θ ' in the most significant moments of the test along a vertical axis adjacent to the gap elements.
As the small dimension of the model, the stress variables depicted are uniformly distributed on the soil layer under geostatic conditions. The introduction of the gap elements does not affect the distribution of both pore pressures and vertical effective stresses, they evolve similarly to what observed in Fig 3a and b. While, the installation procedure strongly affects the development of radial and circumferential stresses, especially during the imbibition phase. The lack of the initial lateral confining allows the soil elements to swell along the radial direction. Therefore, for z lower than the pile length, the changes in σ r ' mainly depend on the variations in suction; for increasing depths, the soil response is under oedometric conditions (see Fig. 7c, line 2). Moreover, the presence of the gap generates an arching-effect in the surrounding soil, which limits the reduction in σ θ ' upon wetting. It, thus, follows that the circumferential stresses are greater than the radial stresses at any depths (see Fig. 7d, line 2). In the in-flight phase, the differences between the two horizontal stress components lessen and vanish below 5 mm (Figs 7c, d, lines 4). This result is not surprising because at this stage the gap is closed and the oedometric conditions are restored over the entire soil layer. From a qualitative point of view, the horizontal stress profiles related to this phase are similar to those reported in Fig. 3. From a quantitative point of view, they differ in the upper part where the changes in effective stresses due to the g-level increase are negligible.
It is worth to notice that the above-mentioned differences in the horizontal stress distributions could markedly affect the pile response. Indeed, the significant soil volume interested by pile foundation under lateral loading is typically limited at shallow depths from the ground level.

Concluding remarks
This paper presents some results of a parametric numerical study devised to reproduce the centrifuge tests carried out by Lalicata et al. [1] on a laterally loaded pile. The attention is focused only on the first two phases, during which the soil model is subjected to an imbibition process at 1xg and, then, to an increase in the gravity loading up to 100xg. The hydro-mechanical behaviour of the B-Grade kaolin is modelled using the MCCM with a modified version of the Gardner retention curve. In this paper, only the results related to the test performed on a loose unsaturated soil are presented and analysed. The main aim of the work was to verify to what extent the soil response is affected by the aftercompaction conditions and the pile installation procedure.
The results of this study point out that: • the MCCM well-reproduces the soil response in terms the evolution over time of the vertical displacements and water inflow; • the uncertainties on the after-compaction conditions appreciably affect the soil response, especially in the imbibition phase. The increase in e 0 entails the reduction in the wetting-induced vertical displacements, and the increase in the soil settlements during the in-flight phase. The changes in S r0 influence only the soil response upon wetting. The reduction in S r0 affects the growth of vertical displacements and water inflow; • the procedure employed for pile installation modifies the final distribution of the stress variables in the upper part of the soil layer. This can deeply affect the response of laterally loaded piles since it mainly depends on the mechanical properties (strength and stiffness) of the soil at shallow depths. Therefore, the results of T6 will be considered as the starting point for the next numerical simulations aimed at reproducing the pile response upon statically loading.