A numerical technique for modeling the behavior of single piles in unsaturated soils

Pile foundations are widely used in both saturated and unsaturated soils. In certain scenarios, these foundations are subjected to combined vertical and lateral loads. Conventionally, saturated soil mechanics principles are routinely used for the design of pile foundations in unsaturated soils. Such approaches contribute to unreliable estimates of the behavior of piles due to ignoring the influence of matric suction. In this paper, a comprehensive numerical technique is proposed for simulating the behavior of single piles subjected to combined vertical and lateral loads in unsaturated soils by taking account of the nonlinear behavior of shear strength and the elastic modulus of unsaturated soils. This is achieved through a subroutine that was developed for use in the ABAQUS software. The proposed numerical method provided reliable prediction of the vertical load-displacement behavior of a published model pile tested in saturated and unsaturated sands. In addition, 3D finite element analysis was extended to simulate the influence of variations in ground water table (GWT) on the vertical bearing capacity and the influence of vertical loads on lateral response of piles. The proposed numerical technique is a promising tool for implementing the state-of-the-art understanding of the mechanics of unsaturated soils into conventional engineering practice.


Introduction
Pile foundations are typically long slender axial columns embedded into the soil to carry large loads safely from the superstructure to the soil below and significantly limit the settlements. These foundations are widely used for supporting structures such as the high-rise buildings, bridges, earth-retaining structures; in certain scenarios, they are subjected to simultaneous vertical loads and lateral loads. In many situations, the entire pile or a portion of the pile is placed in unsaturated soils. However, pile foundation designs are conventionally based on principles of saturated soil mechanics, ignoring the influence of matric suction in unsaturated soils. Due to this reason, both the pile load carrying capacity and their settlement behavior may not be reliably estimated.
The interaction effect of combined vertical and lateral loads on the pile response is also of significant interest in geotechnical engineering practice applications. Several studies are reported in the literature towards understanding the combined influence of vertical and lateral loads on the pile behavior. For example, Jain et al. (1987) [5], Anagnostopoulos & Georgiadis (1993) [6] and Lu & Zhang (2018) [7] reported from laboratory tests that the vertical load contributes to a decrease in the horizontal displacement. The laterally loaded pile behavior is a nonlinear problem and cannot be reliably simulated using two-dimensional (2D) plane-strain analysis (Hazzar et al. 2019 [8]). For this reason, Karthigeyan et al. (2006Karthigeyan et al. ( , 2007 [9] [10], Zheng & Wang (2008) [11], Hussien et al. (2014) [12] and Chatterjee & Choudhury (2016) [13] performed 3D Finite Element Analysis (FEA) and found that the vertical load has significant influence on lateral load capacity of piles embedded in sandy soils. The numerical studies carried out by Hazzar et al. (2017) [14], Zormpa and Comodromos (2018) [15] showed that the influence of vertical load on the lateral pile response is rather limited. The above details suggest that the mechanism of the interaction effect of combined vertical and lateral loads is complex. Most studies in the published literature are related to combined loads action on piles only for saturated soils.
There is limited research that has been devoted towards studying the response of laterally loaded piles in unsaturated soils from experimental and numerical studies (Mokwa et al. 2000 [16]; Stacul et al. 2017 [17]; Lalicata et al. 2019 [18]). To the best of the knowledge of the authors, three-dimensional (3D) numerical modeling studies are not undertaken to predict the loaddisplacement behavior of piles in unsaturated soils taking account of the combined influence of vertical and lateral loads.
Three key objectives are focused in the present study based on numerical modeling studies to model the pile-soil interaction behavior. The first objective was directed towards predicting the behavior of single piles subjected to vertical loading in sand under saturated and unsaturated conditions taking account of the influence of matric suction. The second objective focused to investigate the influence of varying GWTs on the vertical bearing capacity of single piles in unsaturated soils. The third objective was to understand the influence of vertical loads on the lateral response of piles in unsaturated soils. The numerical studies were undertaken extending FEA with the aid of ABAQUS software via user-defined subroutine, USDFLD (User Subroutine to redefine field variables at a material point). The USDFLD facilitates to model the nonlinear behavior of shear strength and the elastic modulus of unsaturated soils and use this information in the prediction of the variation of pile load carrying capacity behavior. The information of the bearing capacity of the single pile also was derived from the load versus displacement curve. The proposed numerical technique in this study is a promising tool for the practicing engineers for implementing the state-of-the-art mechanics of unsaturated soils into practice.

Background
The bearing capacity and settlement behavior of pile foundations are significantly influenced by the shear strength and stiffness of surrounding unsaturated soils. Therefore, the shear strength and modulus of elasticity of unsaturated soils are two key properties used in the development of the finite element model in the present study. Fredlund et al. 1978 [19] provided a theoretical framework for interpreting the shear strength behavior of unsaturated soils in terms of two independent stress state variables; namely, net normal stress (n -ua) and matric suction, (ua -uw), where n is total stress, ua is air pressure, uw is pore water pressure. The variation of shear strength taking into account of matric suction is nonlinear for various soils (Escario et al. 1989 [20], Vanapalli [24]). However, determination of the shear strength of unsaturated soils from experimental studies is time consuming, complex and expensive. For this reason, several researchers proposed procedures for predicting the shear strength of unsaturated soils using the information of effective shear strength parameters and the Soil-Water Characteristic Curve (SWCC). The SWCC represents the relationship between the mass (and/or volume) of water in the soil and the soil suction. Fig. 1 shows the relationship between the shear strength of unsaturated soil and the SWCC at different stages of desaturation. The SWCC can be divided into three main zones as boundary effect zone, transition zone and residual zone. In the boundary effect zone, there is a linear increase in shear strength up to the air-entry value. The contribution of internal friction angle due to matric suction, ϕ b , is equal to internal friction angle ϕ. In the transition zone and residual zone, the contribution of matric suction starts to decrease as the wetted contact of area decreases. The shear strength increases nonlinearly as matric suction increases. ϕ b tends to be lower than the ϕ. Beyond the residual zone, the shear strength of unsaturated soils may increase, decrease or remain constant depending on the type of soil. In the present study, the variation of shear strength of unsaturated soils taking into account of matric suction was predicted using the semi-empirical model proposed by Vanapalli et al. (1996) [21] which requires the information of SWCC and the effective shear strength parameters of saturated soil Eq. 1 and Eq. 2 for predicting shear strength of unsaturated soils were introduced into ABAQUS via user subroutine USDFLD to simulate the non-linear behavior of the unsaturated soil. The ABAQUS software has features to solve complex non-linear problems and has computational capacities to deliver high performance solutions extending FEA. The subroutine USDFLD allows the user to define field variables at every integration point of an element. The material properties can be defined as functions of field variables in a material point. The USDFLD provides access to the solutions from previous increment and allows the user to update field variables for the current increment. The proposed numerical technique accounts for the contribution of matric suction towards shear strength bases on Eq. 1.
Vanapalli & Mohamed (2007) [28] conducted a series of model footing tests on unsaturated soils and concluded that the modulus of elasticity of unsaturated soils is significantly influenced by matric suction. Oh et al. (2009) [29] proposed Eq. 3 to predict the variation of modulus of elasticity of unsaturated soils (i.e. Eunsat) using SWCC and modulus of elasticity of soils under saturated conditions (i.e. Esat) where Pa = atmosphere pressure (101.3 kPa). The fitting parameter is related to the footing size. The fitting parameteris equal to 1 for coarse-grained soils and 2 for fine-grained soils.
The variation of modulus with respect to matric suction is similar to trends shown in Fig. 1 for the shear strength of unsaturated soils.

Numerical Model Simulation
Three key objectives were proposed to be addressed in this study. As mentioned earlier, the first objective was directed towards predicting the vertical load versus displacement behavior of single piles taking account of the influence of matric suction in sand. A model pile test subjected to only vertical loading in Unimin 7030 sand under saturated and unsaturated conditions by Vanapalli et al. (2018) [4] was used to validate the numerical model. Previously, Al-Khazaali et al. 2016 [30] performed FEA using PLAXIS 2D software to predict these experimental results. The average suction values of 2 and 4 kPa of unsaturated conditions were calculated at the centroid of the stress bulb zone below the pile base, which is at a depth of 1.5D. However, in the present study, the variation of shear strength of unsaturated soils corresponding to varying matric suction was considered in the model, which is a rigorous approach. The second objective was directed towards investigating the influence of varying GWTs on the vertical bearing capacity of single piles in unsaturated soils. The pile and soil dimensions were larger than those in the first objective. The ultimate bearing capacity was obtained from load versus displacement curve and was used for the third objective. The focus of third objective was to understand the influence of vertical loads on the lateral response of piles in unsaturated soils. The numerical simulation results using the proposed 3D finite element (FE) model performed on the same Unimin 7030 sand are summarized for addressing the second and third objectives. As per original plans, experimental investigations were also proposed to be undertaken; however, they could not be conducted due to COVID-19 pandemic restrictions. For this reason, only details of FEA studies are presented in this study, for addressing second and third objectives.

Objective I
Experimental studies by Vanapalli et al. (2018) [4] performed in Unimin 7030 Sand were used for validating objective 1. The soil properties of Unimin 7030 Sand are summarized in Table 1. The drying path of SWCC of Unimin 7030 Sand measured using Tempe Cell apparatus (see Fig. 2) was fitted using Fredlund & Xing (1994) [31] equation [4]   where C(ψ) is the correction function, ψ is matric suction, e is Euler's number, and a, n, m are fitting parameters with values of a, n = 4.341, m = 177.27.   [29] suggested that the fitting parameter  varies from 0.5 to 2.5 depending on the foundation size. =1 was chosen in this model. Since sand is non-plastic soil, the fitting parameter  = 1 was used in Eq. 3. An average value of Eunsat along with the pile length was used as uniform elastic modulus of soil in the model. Therefore, Eunsat used for GWT = 450 mm and 650 mm is 6000 kPa and 9100 kPa, respectively. The dilation angle of Unimin Sand was estimated to be 4.5°, which is slightly higher than 10% of DS 415, 1984 [34] The coefficient of earth pressure K0 equal 0.42 was used according to empirical relationship (i.e. K0 = 1-sin from Jaky's 1944 [35]). A surface to surface contact algorithm was used to simulate small sliding and gapping at the pile-soil interface. The pile-soil interface friction angle (´) was assumed to be ´= 2/3ϕ= 24.2°. A displacement of 15 mm was applied on the top of pile to simulate the loading procedure during model pile tests.
Comparisons between the measured data (from Vanapalli et al. 2018 [4]), FEA results using PLAXIS 2D (from Al-Khazaali et al. 2016 [30]) and FEA results using ABAQUS are presented in Fig. 3. It can be seen that numerical results from Al-Khazaali et al. 2016 [30] based on average matric suction values overestimate the load versus displacement behavior of single piles in unsaturated sands. However, reliable comparisons between the numerical simulations and experimental results were obtained by taking account of the influence of shear strength and modulus of elasticity with varying matric suction in ABAQUS software with USDFLD. capacity reduces as the soil desaturates and approaches the residual zone. It can be seen that the calculated bearing capacity using analytical methods underestimates the vertical bearing capacity of single piles. This may be attributed to the theoretical equations that are based only on shear strength parameters do not take into account of the settlement behavior which is strongly related to the modulus of elasticity of unsaturated soils. Fig. 5 shows a schematic diagram of single pile subjected to vertical and lateral loads in sand. Five different GWTs varying from the bottom of pile to the ground level at 0 mm, 150 mm, 300 mm, 450 mm and 600 mm were considered in the FE model. Matric suction profile above the GWT was assumed to be linearly distributed assuming hydrostatic conditions (i.e. ua -uw = -w·hunsat, where hunsat is the height above the GWT). The geometry and mesh of the FE model are shown in Fig. 6. Considering the symmetry of pile-soil system and load direction, half of the 3D FE model was established to simulate the load versus displacement behavior of single piles under vertical and lateral loads. The embedded length (L) and diameter (D) of the pile were 600 mm and 60 mm, respectively. The soil volume with radius of 500 mm and height of 1200 mm was simulated to eliminate the boundary effects. The soil cylinder was discretized with 9648 eight-node stress-pore-pressure coupled brick elements (C3D8P), while the pile was discretized with 528 eight-node brick elements (C3D8). Local mesh refinement was applied to take into account of the high stress variation near the pile shaft and tip. The pore water pressure of all the nodes at GWT surface was set to zero. The Unimin 7030 Sand and pile were modeled following the same procedures discussed in the first objective. The unsaturated modulus of elasticity Eunsat of the sand was determined using Eq. 3. The fitting parameter  and  are both equal to 1. The values of Eunsat for GWT at 150, 300, 450 and 600 mm were 2370, 3760, 5100 and 6720 kPa, respectively. All degrees of freedom at the model base were constrained. The displacements at the lateral boundaries in the radial directions were fixed.

Objective II & III
For addressing the second objective, only vertical load was applied to the top of pile. The ultimate vertical load capacity (Pu) of a single pile can be evaluated from load versus displacement curve using double tangent method as summarized in Prakash & Sharma 1990 [37]. The value of Pu is defined as the vertical load at intersection point of the initial tangents to the load versus settlement curve and the tangent to or the extension to the final portion of the curve.
For addressing the third objective, the loads were applied in two stages. In the first stage, various vertical loads V = 0 Pu, 0.25 Pu, 0.5 Pu and 0.75 Pu were applied on the top of pile. In the second stage, a horizontal displacement of 6 mm (10%D) was applied on the top of pile, while the vertical load was kept constant. Reaction forces of the nodes in the horizontal direction were used to calculate lateral loads.   Fig. 8 shows the relationship between the variation of ultimate vertical bearing capacity pult and matric suction. The matric suction at the centroid of soil layer above GWT was marked as representative matric suction values. Higher matric suction contributes to an increase in soil stiffness and shear strength of unsaturated soils, thereby leading to higher vertical bearing capacity than that of saturated soils with lower stiffness and shear strength.    9.3%, respectively. From these results, it can be concluded that the vertical load has only limited influence on the lateral response of piles. It can also be observed that an increase in the lateral resistance of the pile in unsaturated soils is lower in comparison to saturated soils. The positive effect of vertical load in the lateral response of piles in unsaturated soils with higher strength and stiffness is comparatively less than in saturated soils with lower strength and stiffness.

Summary
In this paper, a comprehensive numerical technique is proposed to simulate the load-displacement behavior of piles in saturated and unsaturated soils under both vertical and lateral loads. This is achieved by introducing a user subroutine USDFLD into ABAQUS to describe the nonlinear behavior of shear strength and modulus of elasticity of unsaturated soils. The required information for modeling includes only the saturated shear strength parameters, saturated modulus of elasticity and the SWCC.
FEA was performed using the proposed numerical technique to predict the load-displacement behavior of single piles that were only subjected to vertical loading. Comparisons were made between the FEA results and a model pile load test from the literature. There is a good agreement between the measured data and the predicted FEA results. In addition, 3D numerical modeling was conducted to investigate the influence of variation of GWTs on the vertical bearing capacity and the influence of vertical loads on the lateral response of piles in unsaturated soils. Experimental studies were planned to provide comparisons with the numerical results; however, they could not be conducted due to recent pandemic restrictions. The numerical results show that the vertical load significantly increases as the GWT depth decreases from the natural ground surface level. It was also found that the vertical load has limited influence on the lateral response of piles in unsaturated sands.
The interaction effect of combined vertical load and lateral load is complex due to the influence of many factors, such as the soil type and pile slenderness ratio. Therefore, further studies experimental studies and numerical investigations are required for better understanding the behavior of piles under combined loads in various types of unsaturated soils.