Load-Settlement Analysis of Axially Loaded Piles in Unsaturated Soils

Abstract

Unsaturated soil covers a significant part of the world, and studying the behavior of deep foundations in this medium is an important step in increasing accuracy and economic efficiency in geotechnical studies. This paper presents an analytical solution to investigate the load-carrying characteristics of single piles embedded in unsaturated soils, accounting for the effect of groundwater level on the pile’s response. For this purpose, relationships for shear modulus and Poisson’s ratio for unsaturated soils were collected from the literature to consider their effects as key parameters on pile performance. A parametric study was conducted to evaluate the effect of soil moisture content on the behavior of the pile-soil system for different soil types, and the effect of pile slenderness on its load-settlement behavior was studied for varying soil moisture contents. The results indicate that the pile stiffness increases as the soil suction increases while below a critical slenderness value, hence increasing the pile load capacity. However, this improvement occurs within a limited range of soil suction that is narrower for coarse-grained soils. The pile settlement corresponding to soil failure was also evaluated by modifying the existing solutions for unsaturated soils. The developed solutions were verified against the predictions of published solutions as well as the results of finite element analysis and pile load tests. It was found that the system stiffness decreases by 50% when the water table rises from the pile toe level to the ground surface in the studied soil.

1. Introduction

In conventional design methods, the soil is either saturated or dry. However, as arid regions cover one-third of the world’s land area [1], many sites are covered with unsaturated soil and rarely—or never—become saturated along the entire depth. This makes it uneconomic to design costly structures such as piles, which are built in saturated soil for the whole life cycle, where this never happens. Besides, higher volumes of construction materials and operations used for this inefficient design would lead to more energy consumption and carbon emissions. In addition, in some cases, variations in the water table levels in the vicinity of foundations may happen for several reasons, like irregular utilization of groundwater, excavation, etc., and should be evaluated using unsaturated soil mechanics. On the other hand, pile foundations are widely used to transfer loads to stiff strata and impose considerable expenses on projects. The overall circumstances make it important to study the effect of unsaturation on pile behavior.

Experimental investigations on the pile shaft load-transfer behavior and ultimate capacity of piles under static axial load in unsaturated soil have been conducted by several studies, considering the suction effect on the shear strength of both cohesive and non-cohesive soils [2,3,4,5,6,7,8,9,10,11]. The results of the experiments agreed on pile capacity growth with the rise of suction and the different scopes of this phenomenon in coarse- and fine-grained soils.

Studies also attempted to verify the test results with analytical and numerical methods. As an analytical approach to this subject, Vanapalli and Taylan [2] modified the conventional α, β, and λ methods for evaluating pile shaft capacity to account for the effect of soil suction. However, their study did not consider pile load-settlement behavior through analytical methods. Several valuable studies evaluated pile settlement employing the classical soil mechanics approach to soil moisture [12,13,14,15]. However, there is no analytical solution to characterize the load-settlement behavior of axially loaded piles embedded in unsaturated soil subjected to the design load, i.e., to evaluate the pile serviceability limit state. Hence, in this study, two analytical solutions will be presented to account for pile load-settlement behavior under service loads and calculate the corresponding settlement to the ultimate shaft-bearing capacity. Moreover, parametric studies have been conducted to recognize the coupled effect of pile slenderness and soil suction on pile-soil system stiffness in six soil types. Further, based on the observed behavior of piles in different unsaturated soils, a critical slenderness ratio was introduced to be considered for design purposes.

Numerical simulation of piles in unsaturated soil has been mostly carried out using the finite element method by implementing Mohr–Coulomb failure criteria for unsaturated soil [3,4,16]. Abed and Vermeer [17] stepped further by considering the failure criteria according to the Barcelona Basic Model (BBM), using model parameters from literature and expert judgment. Wu and Vanapalli [18], however, used the state surface model for the simulation of piles in unsaturated soil to avoid the challenges associated with evaluating the parameter values of BBM. Previous studies have made similar statements about BBM shortcomings [19,20], which will be discussed further. Meanwhile, to verify the analytical model proposed in this study, a numerical simulation was performed using Mohr-Coulomb failure criteria regarding the loading scope of the proposed model.

Evaluation of suction effects on determinative soil parameters on the load-settlement behavior of piles in unsaturated soil has been carried out in this study as well, to have a holistic insight into the evaluation of pile behavior in unsaturated soil. The pile’s load-settlement behavior under service loads is affected by the soil’s small strain shear modulus and its Poisson’s ratio, in addition to the pile’s geometry and Young’s modulus [12]. Acknowledging the recognizable difference between soil and pile Young modulus, Schanz et al. [21] highlighted the effect of soil’s elastic shear modulus as the dominant factor in pile behavior under service loads. This research studied the effect of variations of G in wetting and drying cycles on pile settlement. Moreover, to compensate for the difficulties associated with laboratory testing of G in unsaturated soil, proposed equations for the estimation of this parameter have been collected from the literature and discussed. Then, the most applicable equation in the literature has been opted for and further studied for evaluation of its fitting parameter in different soils.

Overall, the primary purpose of this study is to evaluate the load-settlement behavior of single piles embedded in unsaturated soil under service load and ultimate bearing capacity and assess the effect of soil types and pile slenderness ratio in this regard. Subsequently, soil parameters with a major effect on pile behavior in unsaturated soil and proposed models to evaluate them are studied to help with parameter estimation for design purposes.

Therefore, in Section 2, a model will be introduced to analyze pile load settlement under service loads in contact with the water table and unsaturated soil, starting with a description of the implemented methodology in Section 2.1, presenting the model and the study on the involved parameters in Section 2.2, and the numerical verification of the model in Section 2.3. Further, Section 3 will include parametric studies. In Section 3.1, the effect of soil moisture and pile slenderness on pile-soil stiffness using the proposed model is assessed for six different soil types. Based on that, a critical slenderness ratio is introduced in Section 3.2. Section 4 focuses on the estimation of pile head settlement corresponding to the ultimate shaft bearing capacity in unsaturated soil and the effect of pile slenderness on that. To wrap up, Section 5 will explain the practical usage of the proposed models in practice. In the final section, conclusions are drawn, and further studies are suggested.

2. Load-Settlement Relationship for Pile Embedded in Two-Layered Soil

2.1. Analysis Method

To derive the pile-soil system stiffness in unsaturated soil considering an analytical solution, the equilibrium equation for forces acting on an unsaturated soil element, shown in Figure 1, is considered along the vertical direction. The corresponding equation along the y direction is given by Clifton et al. [22]:

$$\frac{\partial {\tau }_{xy}}{\partial x}+\frac{\partial ({\sigma }_y-u_w)}{\partial y}+\frac{\partial {\tau }_{zy}}{\partial z}+(\gamma -n_w{\gamma }_w-n_a{\gamma }_a)+(n_p+n_c)+\frac{\partial u_w}{\partial y}-F_{cwy}-F_{pwy}-F_{cay}-n_a\frac{\partial (u_a-u_w)}{\partial y}=0$$

(1)

where ${\tau }_{xy}$ is shear stress on x-plane in y direction, ${\sigma }_y$ is the total stress in y direction, $u_w$ is pressure in the water phase, ${\tau }_{zy}$ is shear stress on z-plane in y direction, $\gamma$ is the total unit weight of saturated soil, $n_w$ is porosity with respect to the water phase, $n_w{\gamma }_w$ is gravity body force for the water phase, $n_a$ is porosity concerning the air phase (i.e., percentage of the surface of the element going through air), and $n_a{\gamma }_a$ is gravity body force for the air phase, $n_p$ is the percentage of element surface that is soil particles, $n_c$ is the percentage of element surface that is contractile skin, $F_{cwy}$ is the interaction force between the water and the contractile skin, $F_{pwy}$ is the interaction force between the water and the soil particles, $F_{cay}$ is the interaction force between the air and the contractile skin, $u_a$ is pressure in the air phase, and dx, dy and dz are unit dimensions of the element.

In the derivation of stress matrices from Equation (1), Cliffton et al. [22] stated that “In any case, the inclusion of a porosity term (i.e., a soil property) in the description of the state of stress is not in keeping with continuum mechanics” and ignored the aforementioned terms. Thus, for establishing the load-settlement relationship for piles embedded in soil with the linear elastic response, the equilibrium of the soil element will be reduced to the conventional form.

To obtain a relationship between the soil shear stress acting on the pile shaft and the corresponding deformations, the equilibrium of the soil element is integrated, and the linear shear stress-strain relationship is considered, as stated by Randolph and Wroth [12]. It should be noted that for unsaturated soil, the suction invariants exist only in the diagonal elements of the stress state matrix [22]. Although they influence the shear modulus value, suction invariants are not included in the linear shear stress-strain relationship of unsaturated soils, and the relationship is the same as for saturated soils. Hence, the shear modulus and Poisson’s ratio should be obtained for unsaturated soil conditions using test data or empirical expressions; this will be described subsequently. Based on the abovementioned assumptions and process, pile-soil system stiffness is evaluated in a general form for the water table at an assumed level along the pile shaft. This will be explained in the next section.

2.2. Development of the Proposed Model

In practice, a part of the pile length may be located below the water table level. Figure 2 presents a schematic view of a pile embedded in an unsaturated soil layer overlying a saturated, deep soil deposit. The upper layer is subjected to the capillary effect, and thus the soil profile can be assumed to be two layers in terms of suction pressure, for which there exist saturated and unsaturated zones. In this case, the variation of the water table level and its effect on the pile behavior should be investigated to avoid failure or an overly conservative design. In this section, a load-settlement equation for the pile embedded in two soil layers is developed using the above-mentioned approach.

Considering the soil-pile system shown in Figure 2, for a pile with length L and water table level at nL (0 < n < 1) from the ground surface, the governing differential equation for the pile settlement in each soil layer can be written as [12]:

$$\frac{{\partial }^{2}w}{\partial z^2}=\frac{2}{r_0^2\zeta {\lambda }_1}w(z)$$

(2)

$$w(z)=A{e}^{{\mu }_{1}z}+B{e}^{-{\mu }_{1}z}$$

(3)

where w is pile settlement, z represents depth, r₀ is pile radius, $\mu =\sqrt{\frac{2}{\xi \lambda r_0^2}}$, λ = E_p/G, $\xi =\mathit{ln}\frac{2.5l(1-\upsilon )}{r}$, ν is soil Poisson’s ratio, G is soil shear modulus, L is pile length, E_p is pile material elasticity modulus, $\frac{4}{\eta \left(1-\upsilon \right)}$ denotes the soil reaction at the pile base with η typically taken as 1, and A and B are constants that are determined using appropriate boundary conditions. Subscripts 1 and 2 denote soil layers 1 and 2, respectively (Figure 2).

Equation (3) is first solved for the upper part of the pile length embedded in the unsaturated layer. Applying the first boundary condition for the upper part (nL), and following a procedure described by Randolph and Wroth [12], Equation (3) is solved, and the pile settlement is obtained from:

$$w_2(t)=w_1(z=nL)$$

(4)

$$w_2(t)=\frac{P_b}{r_0{G}_{sat}}\left\{\frac{\eta (1-{\nu }_2)}{4}\cosh ({\mu }_{2}l-{\mu }_{2}nl)+\frac{1}{\pi r_0{\lambda }_2{\mu }_2}\sinh ({\mu }_{2}l-{\mu }_{2}nl)\right\}$$

(5)

where P_b is the soil reaction force transmitted by the pile toe, G_sat is the saturated shear modulus, and w_t is the pile head settlement.

Applying the second boundary condition for axial forces substituting the pile top load gives:

$$P_2(t)=P_1(z=nL)$$

(6)

$$P_2(t)=\frac{2\pi P_b}{{\zeta }_2{r}_0}\left\{\frac{\eta (1-{\nu }_2)}{4{\mu }_2}\sinh ({\mu }_{2}l-{\mu }_{2}nl)+\frac{{\zeta }_2{r}_0}{2\pi }\cosh ({\mu }_{2}l-{\mu }_{2}nl)\right\}$$

(7)

The term P_t refers to the pile top load.

For a compressible pile, the elasticity theory gives:

$$\frac{\partial w}{\partial z}=\frac{-P(z)}{\pi r_0^2{\lambda }_1{G}_{\theta }}$$

(8)

$$\frac{\partial P(z)}{\partial z}=\frac{-2\pi G_{\theta }}{{\zeta }_1}w(z)$$

(9)

where G_θ is the shear modulus corresponding to the volumetric water content θ in soil layer 1.

Differentiating Equation (3) and substituting into Equation (8), and substituting Equation (5) into Equation (3), give:

$$A{e}^{{\mu }_{1}nl}+B{e}^{-{\mu }_{1}nl}=\frac{P_b}{r_0{G}_{sat}}\left\{\frac{\eta (1-{\nu }_2)}{4}\cosh ({\mu }_{2}l-{\mu }_{2}nl)+\frac{1}{\pi r_0{\lambda }_2{\mu }_2}\sinh ({\mu }_{2}l-{\mu }_{2}nl)\right\}$$

(10)

$$A{e}^{{\mu }_{1}nl}-B{e}^{-{\mu }_{1}nl}=\frac{-2P_b}{{\mu }_1{\zeta }_2{r}_0^3{\lambda }_1{G}_{\theta }}\left\{\frac{\eta (1-{\nu }_2)}{4{\mu }_2}\sinh ({\mu }_{2}l-{\mu }_{2}nl)+\frac{{\zeta }_2{r}_0}{2\pi }\cosh ({\mu }_{2}l-{\mu }_{2}nl)\right\}$$

(11)

where A and B in Equation (3) are constants that are obtained by solving Equations (10) and (11) simultaneously.

The pile settlement at depth z is then computed from:

$$\begin{array}{c}w(z)=\frac{P_b}{r_0{G}_{sat}}\left\{\frac{\eta (1-{\nu }_2)}{4}\cosh ({\mu }_{2}l-{\mu }_{2}nl)+\frac{1}{\pi r_0{\lambda }_2{\mu }_2}\sinh ({\mu }_{2}l-{\mu }_{2}nl)\right\}\times \cosh ({\mu }_{1}nl-{\mu }_{1}z)+\\ \frac{2P_b}{{\mu }_1{\zeta }_2{r}_0^3{\lambda }_1{G}_{\theta }}\left\{\frac{\eta (1-{\nu }_2)}{4{\mu }_2}\sinh ({\mu }_{2}l-{\mu }_{2}nl)+\frac{{\zeta }_2{r}_0}{2\pi }\cosh ({\mu }_{2}l-{\mu }_{2}nl)\right\}\times \sinh ({\mu }_{1}nl-{\mu }_{1}z)\end{array}$$

(12)

By integrating Equation (9), the applied load at the pile head (z = 0) is given by:

$$P_t=\frac{-2\pi G_{\theta }}{{\zeta }_1}\times \left[\begin{array}{l}\frac{P_b}{r_0{G}_{sat}}\left\{\frac{\eta (1-{\nu }_2)}{4}\cosh ({\mu }_{2}l-{\mu }_{2}nl)+\frac{1}{\pi r_0{\lambda }_2{\mu }_2}\sinh ({\mu }_{2}l-{\mu }_{2}nl)\right\}\times \frac{\sinh ({\mu }_{1}nl)}{{\mu }_1}+\\ \frac{2P_b}{{\mu }_1{\zeta }_2{r}_0^3{\lambda }_1{G}_{\theta }}\left\{\frac{\eta (1-{\nu }_2)}{4{\mu }_2}\sinh ({\mu }_{2}l-{\mu }_{2}nl)+\frac{{\zeta }_2{r}_0}{2\pi }\cosh ({\mu }_{2}l-{\mu }_{2}nl)\right\}\times \frac{\cosh ({\mu }_{1}nl)}{{\mu }_1}\end{array}\right]$$

(13)

Assuming a constant soil Poisson’s ratio along depth and $\sqrt{\frac{G_{\theta }}{G_s}}=\alpha$, the total pile-soil system stiffness (P_t/w_t) is given by:

$$k_t=G_{sat}{r}_0\frac{\frac{2\pi }{\zeta }{\alpha }^2\frac{L}{r_0}\frac{\tanh ({\mu }_{1}nL)}{{\mu }_{1}L}+\alpha \frac{4}{\eta (1-\nu )}\tanh ({\mu }_{1}nL)\tanh ({\mu }_{2}L-{\mu }_{2}nL)+\frac{2\pi }{\zeta }\frac{L}{r_0}\frac{\tanh ({\mu }_{2}L-{\mu }_{2}nL)}{{\mu }_{2}L}+\frac{4}{\eta (1-\nu )}}{1+\frac{4}{\eta (1-\nu )}\frac{1}{\pi {\lambda }_2}\frac{L}{r_0}\frac{\tanh ({\mu }_{2}L-{\mu }_{2}nL)}{{\mu }_{2}L}+\frac{1}{\alpha }\tanh ({\mu }_{1}nL)\tanh ({\mu }_{2}L-{\mu }_{2}nL)+\frac{1}{\pi {\lambda }_1}\frac{1}{{\alpha }^2}\frac{L}{r_0}\frac{\tanh ({\mu }_{1}nL)}{{\mu }_{1}L}}$$

(14)

Substituting n = 0 and $\alpha =1$ in Equation (14), represents the water table level being at the ground surface when the whole soil is fully saturated. Considering n = 1 and $\alpha =1$ displaying a fully unsaturated situation where the soil below the pile base is also unsaturated and the water table is well below the pile toe, to ignore its effect on the soil-pile response. It is worth mentioning that this model is better suited to be used for bored piles or with due consideration of the installation effect on soil elastic behavior than other methods.

As mentioned earlier, unsaturated soil parameters in Equation (14) should be obtained for the proper suction value using laboratory tests or empirical equations. A contributing attempt in this research study is the collection and evaluation of almost all empirical expressions giving the shear moduli and Poisson’s ratios for unsaturated soils. This attempt and the required considerations for evaluating the soil shear modulus and Poisson’s ratio in unsaturated soil will be described subsequently.

Earlier studies demonstrated that the soil stiffness increases as its suction increases, and therefore, Poisson’s ratio decreases [23,24]. Thota and Vahedifard [24] studied the variation of Poisson’s ratio with saturation degree for 23 soils in the literature and their own tests. The variation of Poisson’s ratio with suction (from dry to wet conditions) was found to be from 4% to 1390% for different soil types, with the least variation belonging to sandy soils and the highest variation for lean clay. The overall effect of this variation on pile load settlement depends on simultaneous shear modulus changes with suction.

The soil shear modulus can be correlated to the suction stress using the soil water characteristic curve (SWCC) [25,26] or other involved parameters. It should be considered that the variation of shear modulus with suction and the SWCC follow different paths in drying and wetting, and the number of cycles affects their variation pattern. The higher the cycle number, the more the SWCC converges to a limited criterion compared to its distance from the first cycle [27]. This behavior differs for different soils; nevertheless, it should be considered while studying pile-soil behavior in unsaturated soils.

To further analyze this issue, the pile settlement variation under drying and wetting cycles can be calculated using Equation (14) for a single-layered soil, as illustrated in Figure 3. A 10 m pile with a 0.5 m radius and a 3000 kN axial load is considered in a mixture of 75% fine sand and 25% kaolin clay classified as SC soil, whose variation of shear modulus under drying and wetting cycles is studied by Ngoc et al. [27]. The Poisson’s ratio is assumed to be constant with suction changes and equal to 0.3 due to the lack of measured or reported data. The sensitivity of the settlement to 66% of the change in Poisson’s ratio was about 3%.

As shown in Figure 3, the pile settlement changes from 45% to 238% under constant axial load during the drying and wetting cycles. Such changes fall into a certain range after the first cycle, and they become about 33% larger than the saturated condition in the first wetting cycle. This should be considered while studying the response of axially loaded piles embedded in unsaturated soil, for which the effect of drying and wetting cycles on the soil-pile interaction becomes of concern. When drying and wetting cycles are experienced, the overall pile-soil stiffness range is altered by this interaction.

To establish representative shear modulus and Poisson’s ratio for the unsaturated soil condition, the authors reviewed some correlations of shear modulus and Poisson’s ratio variation with the suction stress available in the literature, as shown in

Table 1. Shear modulus and Poison’s ratio relations with matric suction.

Nomenclature
Gunsat—unsaturated small strain shear modulus		μunsat—unsaturated Poisson’s ratio
Gdry—small strain shear modulus at dry condition		μsat—unsaturated Poisson’s ratio
Gsat —saturated small strain shear modulus		μd—unsaturated Poisson’s ratio
G0 = the shear modulus at 1 kPa confining pressure		Sr—the soil saturation degree
$E_{unsat}$—unsaturated Young’s modulus		Vs—shear wave velocity in the soil
$E_{sat}$—saturated Young’s modulus		Gs—specific gravity
OCR —over consolidation ratio		K0—coefficient of lateral earth pressure at rest
e —void ratio		ρ—soil density
$(u_a-u_w)$—soil suction		Patm—atmospheric pressure
$\left(\sigma -u_a\right)$—principal net stress		$p_r$—the reference pressure
$\left({\sigma }^{\prime }-u_a\right)$—mean effective stress		$f\left(e\right)=\frac{1}{(0.3+0.7e^2)}$—the void ratio function proposed by Hardin and Black [28].
${\left({\sigma }^{\prime }-u_a\right)}_c$—mean apparent reconsolidation stress
Reference	Proposed equation		Description
Wu et al. [29]	$\begin{array}{l}{G}_{unsat}=\left[1+H(S_r)\right]G_{dry}\\ H(S_r)=\left\{\begin{array}{cc}(a-1)\sin (\frac{\pi S_r}{2b})& S_r\le b\\ (a-1)F(S_r)& S_r>b\end{array}\right.\\ F(S_r)=\frac{1}{2}{(\frac{S_r-b}{100-b})}^2·\sin \left\{\frac{\pi }{100-b}\times \left[S_r+50-\frac{3}{2}b\right]\right\}\end{array}$ - a is the maximum value of $G_{unsat}/G_{dry}$ - b is the optimum saturation		Tested on fine-grained non-cohesive soil
Pereira and Fredlund [30]	${\mathsf{\mu }}_{unsat}=0.3+\frac{{\mu }_{sat}-0.3}{\left[1+{\left(\frac{(u_a-u_w)}{c}\right)}^b\right]}$ - $c=c_{1}x{\left(\sigma -u_a\right)}^2+c_{2}x\left(\sigma -u_a\right)+c_3$ - $c_1=9.39\times {10}^{-4} c_2=7.46\times {10}^{-2} c_3=4.07$ - $b=b_{1}x{(\sigma -u_a)}^{b2 } b_1=49.01 b_2=-6.10\times {10}^{-1}$		Validated for SM-ML soil
Cho and Santamarina [31]	$\begin{array}{l}{[V_s]}_s={[V_s]}_{s=1}{\left[1+\frac{2{\sigma }_{eq}^{\prime }}{(1+K_0){\sigma }_v^{\prime }}\right]}^{\beta }\sqrt{\frac{e+G_s}{eS+G_s}}\\ G_s=\rho {[V_s]}_s\end{array}$ - σ′eq is the equivalent effective stress due to capillarity at a given Sr		Validated for 4 cohesionless soil mixtures
Manusco et al. [32]	$\frac{{\left(G_{unsat}\right)}_{s\le S_{ev}}}{P_a}=A{\left[\frac{{(p-u_a)}_C+s}{P_a}\right]}^n{OCR}^m (u_a-u_w)<air entry value$ - ${(p-u_a)}_C$ is the mean net stress of the drying path under consideration - A is the stiffness index - n and m are the stiffness coefficients representing the stiffness of the material under the reference pressure (at an over consolidation ratio (OCR) of 1) and the sensitivity of the stiffness to stress state and history, respectively. $G_{unsat}={\left(G_{unsat}\right)}_{{(u_a-u_w)}^{*}}\left\{\left[1-r\right]e^{-\beta (s-s^{*})}+r\right\} \left(\right(u_a-u_w)-{(u_a-u_w)}^{*})\ge 0$ - β is the parameter that controls the rate of increase of soil stiffness with increasing suction and is associated with the soil sensitivity to suction changes - r is the ratio between shear stiffness at $(u_a-u_w)$* and the threshold value of Gsat for increasing suction and is related to the stiffness increase due to menisci water. - Molding water content seems to cause some significant differences in soil response. This variable affects the ratio between the saturated value and the unsaturated threshold values of shear stiffness, and the suction values, called ${(u_a-u_w)}^{*}$.		Tested on silty sand
Ng et al. [33]	$G_{unsat\left(ij\right)}=C_{ij}^{2}F\left(e\right){\left[\frac{({\sigma }_i-u_a)}{p_r}\frac{({\sigma }_j-u_a)}{p_r}\right]}^n\times {\left[1+\frac{(u_a-u_w)}{p_r}\right]}^{2b}$ - $C_{ij}$ is a constant reflecting the inherent soil properties in the ij plane, with dimensions m/s - $({\sigma }_i-u_a)$ and $({\sigma }_j-u_a)$ are principal net stresses in the ij plane - n is empirically derived constant - b = empirical exponent reflecting the influence of matric suction on the shear modulus.		Tested on ML soil
Sawangsuriya et al. [25]	${\left(G_{o, us}\right)}_{hh}=Af\left(e\right){\left({\sigma }_0-u_a\right)}^n+C{\theta }^{\kappa }(u_a-u_w)$ - A, C and κ are fitting parameters - n = 0.5		Fitting parameters were evaluated for 9 cohesive and cohesionless soils.
Sawangsuriya et al. [25]	${\left(G_{o, us}\right)}_{hh}=Af\left(e\right){{[\left({\sigma }_0-u_a\right)}^n+{\theta }^{\kappa }(u_a-u_w)]}^n$ - A, n and κ are fitting parameters		Fitting parameters were evaluated for 9 cohesive and cohesionless soils.
Biglari et al. [34]	$\begin{array}{l}{p}^{″}=(p-u_a)+S_r(u_a-u_w)\\ OC{R}_{p^{″}}=\frac{p_{\max }^{″}(\zeta )}{p^{″}(\zeta )}\\ \zeta =f(s)·(1-S_r)\\ h(S_r)=1-a^{\prime }[1-\exp (b^{\prime }\frac{\zeta }{f(s)})]\\ G_0=A_0{p}_{atm}{(\frac{p^{″}}{p_{atm}})}^{n_0}{(OC{R}_{p^{″}})}^{m}f(e)h(S_r)\end{array}$ - A0 is a dimensionless stiffness index in a small strain range - $P^{″}$ is the average skeleton stress - $P_{max}^{″}$ average skeleton stress at the yielding point - n0 is a stiffness coefficient that accounts for the effect of $P^{″}$ on stiffness in small strain range - OCRp″ is the over consolidation ratio in p″:ξ plane - ξ is the bonding variable - f(s) is a function that depends on the size of the particles and the value of the water surface tension - m accounts for the effect of OCR on the stiffness - h(Sr) is a non-dimensional function of the degree of saturation defined above, where a′ and b′ are fitting parameters		Fitting parameters evaluated for Zenoz kaolin
Khosravi and McCartney [35]	$G_{unsat}=A{P}_a{\left[\frac{p_c^{\prime }0}{\left(\sigma -u_a\right)}\exp \left(\frac{∆e^p}{\lambda -\kappa }\right)\right]}^K{\left[\frac{\left(\sigma -u_a\right)}{\left({\sigma }^{\prime }-u_a\right)}\exp \left(b\left|S_{e0}-S_e\right|\right)\right]}^{K^{\prime }}{\left(\frac{\left({\sigma }^{\prime }-u_a\right)}{P_{atm}}\right)}^n$ - A and n are experimental parameters - $p_c^{\prime }0$= initial mean apparent pre-consolidation stress - $∆e^p$= plastic change in void ratio - $\lambda and \kappa$ the slopes of the loading and unloading branches of the isotropic compression curve - K and $K^{\prime }$= hardening parameters - b = double hardening parameter referring to the effects of changes in soil saturation - $S_e$ residual degree of saturation - $S_{e0}$ initial residual degree of saturation		Tested on low plasticity soil
Oh and Vanapalli [36]	$G_{\max (uns)}=G_{\max (sat)}\left[1+\zeta \left(\frac{(u_a-u_w)}{P_{atm}/101.3kPa}\right)(S^{\xi })\right]$ - ξ and ζ are fitting parameters		Evaluated for sandy soil behavior and fitting parameters evaluated for 5 sandy soils
Lu and Kaya [26]	$G_{unsat}=G_d+(G_w-G_d){\left(\frac{\theta -{\theta }_d}{{\theta }_w-{\theta }_d}\right)}^m$ - w and d indices are related to the wet and dry states of soil - θ = volumetric water content m is an experimental fitting parameter		Fitting parameter evaluated for 16 cohesive and cohesionless soils
Georgetti et al. [37]	$G_{uns}=G_{sat}+\frac{u_a-u_w}{m+n·(u_a-u_w)}$ - m and n are fitting parameters		Fitting parameters evaluated for clayey medium to fine sand of colluvial origin
Adem and Vanapalli [38]	$\frac{E_{unsat}-E_{sat}}{P_{atm}}=10.74{\left(X\right)}^{2.62}$ - X is an experimental constant		Evaluated for inflatable soil
Wong and Mašín [39]	$\begin{array}{l}{(\sigma -u_a)}^{\prime }=(p-u_a)+\chi (u_a-u_w)\\ {\lambda }_p=\frac{\gamma }{\ln {\chi }_0}\ln \left[({\chi }_0{}^{\frac{{\lambda }_{p}0}{\gamma }}-{\chi }_0){(\frac{e}{e_0})}^{(\gamma -1)}+{\chi }_0\right]\\ {\chi }_0={(s_{e0}/s)}^{\gamma }\\ G_0=p_r{A}_g{e}^{-mg}{\left(\frac{{(\sigma -u_a)}^{\prime }}{p_r}\right)}^{n_g}{S}_r{}^{\frac{-k_g}{{\lambda }_p}}\end{array}$ - λp is the slope of the water retention curve - Ag, ng, mg, kg and γ are fitting parameters		Fitting parameters evaluated for completely decomposed tuff
Hoyos et al. [40]	$\log \left(G_{unsat}\right)=\log \left(G_0\right)+\log {(p-u_a)}^B+\log {(u_a-u_w)}^H$ - H and B are experimental parameters - G0 = the shear modulus at 1 kPa confining pressure - $(p-u_a)$ is net confining pressure - $(u_a-u_w)$ is matric suction		Obtained for SM soil
Zhou et al. [41]	$\begin{array}{l}\zeta =f(u_a-u_w)(1-S_r)\\ f(u_a-u_w)=\frac{3T_s}{Rs}\frac{(\sqrt{9+8Rs/T_s}-3)(\sqrt{9+8Rs/T_s}+1)}{16}\\ K=\frac{\upsilon p^{*}}{\kappa }\\ G_{unsat}=C_0(1+e^{-3})\left({(\frac{{(\sigma -u_a)}^{\prime }}{p_{atm}})}^{0.5}+n_s{\zeta }^{0.5}\right)\end{array}$ - Ts is the surface tension coefficient of water that is equal to 72.8 mN/m at 20 °C - R is the radius of spherical particles - p* is the soil effective stress - κ is the slope of the isotropic unloading-reloading line (URL) in the $\nu -\mathrm{l}\mathrm{n}\left({\sigma }^{\prime }-u_a\right)$ plane C0 and ns are experimental parameters		Fitting parameters evaluated for completely decomposed Tuff
XU and Zhou [42]	$G_{unsat}=C_0{(1+e)}^{-3}[{\left(\frac{\left({\sigma }^{\prime }-u_a\right)}{p_r}\right)}^{0.5}+C_s{\xi }^{0.5}]$ - C0 is a constant reflecting the effect of the inherent soil structure on Gunsat. - Cs describes the effect of additional normal forces between soil particles provided by water menisci (ξ) on Gunsat - ξ = f(s)(1 − Sr), f(s) is presented against suction for different r (radius of the spherical particles) values in a graph		Validated for compacted clayey silt and sand
Dong et al. [43]	$G_{unsat}=A_0{\left(\frac{1}{S_e}\right)}^{\beta }·{\left(\frac{\left({\sigma }^{\prime }-u_a\right)}{P_{atm}}+1\right)}^{{\gamma }_0}$ - $A_0$= an experimental fitting parameter having the same units as the small-strain shear modulus - $S_e$= effective degree of saturation - ${\gamma }_0$ and $\beta$ are experimental fitting parameters		Fitting parameters were evaluated for 29 cohesive and cohesionless soils.
Cao et al. [44]	$\begin{array}{l}{V}_s=\frac{{(a-e)}^2}{1+e}{\left(\frac{({\sigma }_c-u_a)-(u_w-u_a){(\frac{u_w-u_a}{u_w-u_{a0}})}^{-0.55}}{1kPa}\right)}^n\\ G_{uns}=\rho V_s\end{array}$ - a and n are fitting parameters		Fitting parameters evaluated for copper tailing
Liu et al. [45]	$\begin{array}{l}{G}_{uns}=AF(e)F(w_c){\left(\frac{{\sigma }^{\prime }}{P_a}\right)}^n\\ F(w_c)={(100w_c)}^{-\alpha }\end{array}$ - wc is water content in decimal - A, α and n are fitting parameters		Fitting parameters evaluated for Yan’an loess
Khosravi et al. [46]	$G_{unsat}=A{P}_a{\left[\frac{p_c^{\prime }0}{p_n}\exp \left(\frac{∆e^p}{\lambda -\kappa }\right)\right]}^K{\left[\frac{p_n}{p^{\prime }}\exp \left(b\left|S_{e0}-S_e\right|\right)\right]}^{K^{\prime }}{\left(\frac{p^{\prime }}{P_a}\right)}^n$ - All the parameters are the same as in Khosravi and McCartney [35], except $S_e$ which is obtained by $S_e=S_{app}-S_a$ - $S_{app}$ is apparent saturation calculated by a formula presented in the reference paper and $S_a$ is effective saturation of trapped air		Validated for 4 sets of silty and SP-SM soil
Pham [47]	$G_{unsat}=p_{r}f\left(e\right)\left[b_1{\left(\frac{(u_a-u_w)S_r}{p_r}\right)}^{m_1}+c{\left(1-S_r\right)}^k\right] for S_r\ge S_{rs}$$G_{unsat}=p_{r}f\left(e\right)\left[b_2{\left(\frac{(u_a-u_w)S_r}{p_r}\right)}^{m_2}+c{\left(1-S_r\right)}^k\right] for S_r<S_{rs}$ - $S_{rs}$ = degree of saturation corresponding to the shrinkage limit and can be determined from the results of the shrinkage test - $b_1$ and $m_1$ are model parameters for evaluating the contribution of matric suction to $G_{unsat}$ when $S_r\ge S_{rs}$, ${,b}_2$ and $m_2$ are model parameters for evaluating the contribution of matric suction to $G_{unsat}$ when $S_r<S_{rs}$ - c and k are model parameters for evaluating the combined contribution of van der Waals attractions and electric double-layer repulsions		Fitting parameters evaluated for 4 cohesive soils
Mahmoodabadi and Bryson [48]	$G_{unsat}=G_{sat}{[1+{\left(\frac{e_i^{{\psi }_i^{\prime }}(u_a-u_w)}{{\omega }_i}\right)}^{n_i}]}^{m_i}$ - ${\psi }_i^{\prime },n_i,{\omega }_i\mathrm{a}\mathrm{n}\mathrm{d}{m}_i$ are modified SWCC parameters developed in the reference paper. - $e_i$ is the void ratio of the soil corresponding to the original net normal stress		Validated for 13 cohesive and cohesionless soils
Thota et al. [24]	${\mu }_{unsat}={\mu }_{dry}+({\mu }_{sat}-{\mu }_{dry})\left[1-{\left\{1+{\left(\frac{(u_a-u_w)}{1.04S_r}\right)}^{0.19n+2.71}\right\}}^{\frac{1}{0.19n+2.71}-1}\right]$$n$ is the pore size distribution parameter		Validated for 5 cohesive and cohesionless soils
Yan et al. [49]	$\begin{array}{l}{G}_{\max }=\left[1+H(S_r)\right]G_{dry}\\ H(S_r)=\left\{\begin{array}{cc}(a-1)\sin (\frac{\pi S_r}{2b})& S_r\le b\\ (a-1)F(S_r)& S_r>b\end{array}\right.\\ F(S_r)={(1-\frac{S_r-b}{100})}^2·\sin \left\{\frac{\pi }{100-b}\times \left[S_r+50-\frac{3}{2}b\right]\right\}\end{array}$ - a is the maximum value of $G_{unsat}/G_{dry}$ - b is the optimum saturation		Tested on silty fine sand

. As seen in some studies, unsaturated shear modulus was correlated with suction by considering the simultaneous effects of confining pressure. In some other studies, other factors, including loading plane, inherent soil structure, overconsolidation ratio, and hardening parameters, were considered.

According to Table 1, most equations contain parameters that are not usually measured in the laboratory. For example, Manusco et al. [32] correlated unsaturated shear modulus using the Barcelona Basic Model’s parameter β and stiffness index and coefficients from the literature. Khosravi and McCartney [35] presented a model using unsaturated soil’s hardening parameters. Adem and Vanapalli [38] presented a model dedicated to inflatable soil and relationships developed by [40,41,42,43,47]. Khosravi and McCartney [35] included two or more fitting parameters that limit their application scope.

In general, most of the equations gathered for estimating unsaturated shear modulus require complicated tests and calculations that involve several experimental fitting parameters. In addition, some equations and their curve-fitting parameters have been evaluated based on limited soil data. Hence, they may not apply to practical foundation design. Meanwhile, the correlation proposed by Lu and Kaya [26] for evaluating the shear modulus of unsaturated soil is:

$$G=G_d+(G_w-G_d){\left(\frac{\theta -{\theta }_d}{{\theta }_w-{\theta }_d}\right)}^m$$

(15)

where subscripts d and w stand for dry and wet conditions, respectively, $\theta$ is volumetric water content, obtained by multiplying the degree of saturation (S_r) by porosity (n), and m is a fitting parameter evaluated for the abovementioned 16 soils. In Equation (15), characters G and $\theta$ with no subscript are unsaturated values around the average of the dry and wet shear moduli.

Equation (15) has only one curve fitting parameter, which has a clear conformity with the soil’s specific surface area. This correlation was verified for 16 different soils with reasonable results. These soils included low-plasticity clay (CL), high-plasticity clay (CH), poorly graded sand (SP), clayey sand (SC), low-plasticity silt (ML), and organic silt (OL) according to the Unified Soil Classification System (USCS), with a range of shear modulus from 6.15 to 269.6 MPa, a plasticity index from zero to 73, and a specific surface area of 2.04 to 122.16 m²/gr [26]. It is noted that Ghazavi et al. [50] used Equation (15) for the determination of the reaction modulus coefficient and soil-pile stiffness system for laterally loaded piles in unsaturated soils successfully and obtained reasonable results.

In Equation (15), the fitting parameter “m” shows how sensitive the soil’s mechanical properties are to changes in its moisture content. Greater “m” means more sensitive soils, with low sensitivity for cohesive soils and high sensitivity for cohesionless soils. For example, m = 0.01 and 0.2 for the Ottawa and Esperance sands, respectively. Both soils are classified as poorly graded sand (SP) and have almost the same characteristics except for their specific surface areas of 11.25 m²/gr and 2.04 m²/gr for Esperance sand and Ottawa sand, respectively. Meanwhile, silty soils have greater m values, e.g., m = 0.84 for BALT silt and 0.8 for Bonny silt. Three tested clayey soils showed highly non-linear variation behavior for Young’s modulus with moisture content variation and larger “m” values compared to other soil types, i.e., 0.96 for Denver bentonite, 1.37, and 1.63 for low plasticity clays [26].

To further assess the “m” value, this study gathered data from the literature on the relationship between shear modulus and suction for 14 distinct soils. The “m” value was then computed for each soil, as shown in

Table 2. “m” value for different soil types.

Reference	Soil Type	USCS Classification	Ip	Net Confining Stress (kPa)	m
Hoyos et al. [40]	Silty sand	SM	5	400	0.790
				200	1.270
				100	1.507
				50	1.529
Espitia et al. [51]	Lower Hauterivian Barremian claystone	claystone	10	0	0.763
Zhan [52]	Compacted expansive soils from Zao-Yang, China	CH	31	50	0.925
				200	0.992
Miao et al. [53]	Compacted expansive soils from Nanyang, China	CH	31.8	25	1.070
				62.5	1.094
				112.5	1.838
Miao et al. [54]	Compacted expansive soils from Guangxi, China	CH	31.1	50	0.725
				200	0.725
Khosravi and McCartney [35]	Bonny Silt	ML	4	100	0.624
				150	0.352
				200	0.24
Miller and Muraleetharan [55]	Minco Silt	CL	8	39	2.439
Ng et al. [33]	Completely decomposed tuff from Hong Kong	ML	14	110	0.893
				300	0.61
Zhou et al. [41]	Completely decomposed tuff from Hong Kong	ML	14	100	1.566
Khosravi et al. [46]	Bonny Silt	ML	4	40	1.298
				200	1.055
				600	0.806
	Ardabil silt	ML	8	40	0.437
Ngoc et al. [27]	A mixture of 75% fine sand and 25% kaolin clay	SC	9.8	20	1.238
Khosravi et al. [56]	F-75 blast furnace sand	SW		3.50	0.553
				12	0.423
				20	0.308
Khosravi et al. [57]	Bonny Silt	ML	4	125	1.350
				175	0.795
				225	0.782
Takkabutr [58]	High plasticity clay from southeast Arlington	CH	37	0	1.755
				6.8	1.488
				17.2	1.350
				34.5	1.342
	Poorly graded sand	SP		6.8	0.602
				17.2	0.584
				34.5	0.565

. Given that Equation (15) does not account for the influence of the confining pressure on the unsaturated shear modulus, the findings clearly indicate that the “m” value relies on the net confining pressure and diminishes with increasing confining pressure. Greater confining pressures result in a reduction in the soil’s sensitivity to variations in suction. Therefore, there is no singular “m” value that represents the relationship between shear modulus and suction for each soil type, and the change pattern differs depending on the applied confining pressures.

Figure 4 shows the variation of “m” with confining pressure for each soil. Figure 4a,b,d,e show the decreasing trend of “m” with the increment of net confining pressure for silty sand, high-plasticity clay from southeast Arlington, completely decomposed tuff from Hong Kong, F-75 blast furnace sand, and poorly graded sand. As an exception, the tests performed by Zhan [52], Miao et al. [53], and Miao et al. [54] on expansive soils show that the “m” value does not decrease with increasing confining pressure, and for the Nanyang soil with the highest free swelling of 74%, the “m” value has a sharp increase with an increment in confining pressure.

Figure 4f illustrates how the calculated or measured SWCC can affect the calculated “m” value for a single soil. As seen, the variation of “m” for bonny silt from three studies results in different curves and slopes. This is due to the differences between the presented SWCCs in these studies and, consequently, the different saturation degrees and volumetric water contents used in Equation (15). This means that the effect of drying and wetting cycles should also be considered while choosing the volumetric water content for using Equation (15). As the SWCC diagram has different paths for drying and wetting cycles, the proper cycle should be chosen regarding the soil and pile conditions.

It should be noted that in Equation (15), a single “m” value is used for predicting the Young and shear moduli [26]. This means that the variation pattern of both parameters with suction is assumed to be the same, and the effect of simultaneous variation of Poisson’s ratio with suction is neglected. In addition, this relationship should be used in the monotonic range of variation of shear modulus with matric suction. This limits its application in sandy soils and thus requires further investigation.

2.3. Verification of Proposed Model for Two-Layered Soil with Numerical Data

To validate the developed approach and the predictions of Equation (14), finite element analyses were performed employing a computer program coded into Plaxis 2D v8.5 [59]. For this purpose, an axisymmetric finite element model was established for a solid concrete pile with 10 m length and 1 m diameter. Three water table levels were considered in the analysis: −3, −6, and −10 m from the ground surface. The variation of suction above the water table is assumed to be linear with zero air pressure.

The decision on choosing the material model and the extent of simulating unsaturated soil behavior was made according to the scope of the required outputs in this study and the findings of previous studies. The numerical simulation in the current study is used to verify the applicability of Equation (14)—which analyzes pile behavior under service loads in an elastic region—for a pile in contact with a water table.

There are studies which have used unsaturated soil material models for the simulation of pile behavior. Abed and Vermeer [17] used the Barcelona Basic Model (BBM), derived the model parameters from the literature, and filled the gaps with expert judgment. Regarding the difficulties associated with parameter evaluation for BBM, the discourse on this subject has had remarkable results. D’onza et al. [19] performed laboratory tests in seven distinct research teams to evaluate BBM parameters for low-plasticity clay. They reported significant variation between parameter values derived by each team regarding the interdependency of parameters and their coupled effect on unsaturated soil behavior. For instance, the variation of the values attained for the reference pressure ($p^c$) was 2 × 10²³. This extremely wide difference is rooted in the dependency of this parameter on other estimated parameters and the method for doing so. Another unsaturated soil model used for the simulation of pile behavior is the state surface model, conducted by Wu and Vanapalli [18]. They simulated the hydro-thermo-mechanical behavior of piles, which falls outside the scope of this research, focusing on the estimation of pile behavior under service loads.

Therefore, the implementation of an elastoplastic material model while considering the suction effect on the mechanical properties of soil was considered to fulfill the expectations. This assumption was aligned with previous research outputs. The Mohr-Coulomb criterion has been used by several studies simulating piles in unsaturated soil [3,4,16]. Al-khazaali and Vanapalli [4] highlighted that this model may not be suitable for simulating soil hardening behavior. However, the present research focuses on service loads with small settlements for piles. Hence, the Mohr–Coulomb criterion was utilized to simulate the soil material in the finite element model, and the pile was considered a linear elastic material.

To consider the contribution of suction variation with depth on soil properties and to keep it consistent with the Mohr-coulomb model, the soil above groundwater was divided into different layers, and the corresponding unsaturated soil elasticity modulus and cohesion to each suction level was used for them. A similar approach to modeling unsaturated soil with the Mohr–Coulomb model was implemented in previous simulations [4].

The soil considered in this analysis is known as Minnesota lean clay, which is classified as low plasticity clay (CL) according to the Unified Soil Classification System (USCS). The soil shear modulus changes with suction, as presented in Figure 5, as reported by Sawangsuriya et al. [25].

The corresponding elastic modulus for each suction level was obtained from Figure 5, using Equation (16) and ν = 0.3.

$$G=\frac{E}{2\left(1+\nu \right)}$$

(16)

For this soil, the saturated cohesion and the friction angle were 10 kN/m² and $\varphi =28°$, respectively. The cohesion of unsaturated soil (c_(us)) is calculated considering the cohesion of saturated soil (c_(sat)) using the soil water characteristic curve (SWCC) given in [2] as:

$$c_{(us)}=c_{(sat)}\left(1+\frac{(u_a-u_w)}{(P_a/101.3)}(S_r{}^{a1})/a2\right)$$

(17)

where (u_a − u_w) is matric suction, S_r is saturation degree, P_a is the atmospheric pressure, a1 is a fitting parameter equal to 2 for fine-grained soils and soils with a plasticity index PI > 15 (PI = 24% for the soil considered in this example).

The contribution of suction to the soil strength is added to the unsaturated cohesion value as an apparent cohesion [60], i.e.:

$${\tau }_{us}=(u_a-u_w)(S_r^{\kappa })\tan {\delta }^{\prime }$$

(18)

where $\kappa =-0.0016I_p^2+0.0975I_p+1$ and ${\delta }^{\prime }$ is the interface friction angle, which is assumed to be $0.8\varphi$ in this study.

Figure 6 shows the generated finite element mesh with 15-noded triangle elements for the model with a water table at −6 m from the ground surface. The model shown in Figure 6 has 399 elements. The mesh has been divided into 3 major clusters with distances equal to 3 and 5 times the pile diameter, as suggested by Schmüdderich et al. [61]. The mesh setting was the same for the two other models, with water tables at −3 and −10 from the ground surface. Based on a sensitivity study, this mesh was found to be fine enough to reflect the model’s state of stress due to the negligible difference with a finer level of mesh coarseness. The vertical boundaries were prevented from horizontal displacement, and the bottom boundary was fixed in all directions (i.e., no displacements). The boundaries were extended to 20 m in the horizontal and 10 m in the vertical directions to eliminate the effect of the boundaries on the calculated response.

In this study, the pile-soil interface shear resistance was studied for similar soils in the literature. Hamid and Miller [62] studied the shear strength of the rough and smooth interface in an unsaturated lean clay with mechanical characteristics similar to the soil modeled in this study. The maximum peak-to-valley height (R_max) they considered was 0.38 mm for the rough surface and 0.0025 mm for the smooth counterface. As the roughness of the surface of the concrete bored pile in this study is closer to the tested rough surface by Hamid and Miller [62], the interface reduction factor in this study was considered the same as for the rough surface in the test, equal to unity. Wehnert and Vermeer [63] also adopted the same approach to choosing the interface reduction factor.

While installing a pile, extra residual stresses may arise from soil displacement, specifically deformations occurring in the installation of displacement piles. The phenomenon of soil setup, leading to enhanced capacity, can also be noted due to the dissipation of pore pressure during soil reconsolidation [64]. Therefore, as the model is used to verify Equation (14), the pile is assumed to be a bored pile in this model. The simulation stages were chosen according to the study by Schmüdderich et al. [61]. Starting with applying the gravity forces to the model. Next, the borehole is simulated by deactivating soil material in the hole and providing support to the hole surface. The same approach to providing support for voids in numerical simulation is adopted by Liu et al. [65] as well. In the next step, the pile material is substituted in the hole, and the model is analyzed under gravity loading. Finally, the prescribed displacement is applied to the pile head.

Figure 7 shows vertical displacement vectors for the FE analysis of a pile embedded in the ground with the water table level at −3 m from the ground surface. Figure 8 and Figure 9 present vertical displacement sections for the soil around the pile shaft and below the pile base at the groundwater table levels of −3 and −10 m. As seen, when the groundwater level is at −10 m, the soil around the pile shaft is fully mobilized, and the displacement is approximately uniform along the pile shaft. This is due to the large shear stress in this condition compared with the state when the water table level is −3 m. Moreover, for the soil below the pile base, it takes a longer distance for the vertical displacement to reach zero because of the large load concentration.

Figure 10 shows the numerical predictions for the total stress for the water table levels at −3, −6, and −10 m from the ground surface. It is observed that the total stress at the pile toe has nearly twice the value for the water table level at −3 m compared to the same value for the water table level at −10 m from the ground surface. This agrees with the results for pile capacity in these groundwater conditions.

Figure 11 shows the comparison of the load-settlement curves for the example pile calculated from Equation (14) and from the numerical simulation for the three water table levels considered in the analysis. As seen in Figure 11, there is an acceptable agreement between the numerical and analytical results calculated using Equation (14) within the elastic range. In addition, Figure 11 demonstrates that there are significant improvements in the pile stiffness and its load-carrying capacity as the water table elevation drops.

Figure 12 presents the variation of the pile-soil system stiffness with water table elevation, represented by nL. It is observed in Figure 12 that the stiffness of the pile-soil system more than doubles as the water table is at the pile tip level compared with when it is at the ground surface. The good agreement with the finite element results verified the present solution capability (Equation (14)) and its ability to simulate the pile performance for both saturated and unsaturated soil conditions.

3. Parametric Studies

3.1. Effect of Soil Moisture and Pile Slenderness on Pile-Soil Stiffness

In this section, Equation (14) from the developed solution is employed to investigate the variation of pile-soil system stiffness with the soil moisture and pile slenderness ratio for different soil types, considering a fully unsaturated situation (n = 1). To perform these analyses, it is necessary to consider the shear modulus and Poisson’s ratio changes with the soil moisture contents and consequently with the suction stress. To find out the variation of shear modulus with the soil moisture content for poorly graded sand (SP), silty sand (SM), clayey sand (SC), low-plasticity silt (ML), low-plasticity clay (CL), and high-plasticity clay (CH) soils. The authors collected data reported in previous studies [25,40,66] as presented in Figure 13. This figure represents the variation of shear modulus with suction in the drying path. At the beginning of the drying path, the existing water in pores resists more due to high surface tension and strong molecular bonds. As drying proceeds, the suction stress increment decreases for coarse-grained soils because of their larger pore sizes and sparser pore water compared with fine-grained ones. On the other hand, the suction stress increment is more uniform in fine-grained soils, which can lead to a wider range of shear modulus improvements.

The data for soils were obtained from tests conducted employing an apparatus that could apply two stress state variables, the net confining pressure and matric suction, during the measurements of the soil shear modulus with bender elements [25,40,66].

Figure 14 illustrates the variations in the stiffness of the pile-soil system ($k_t$), the pile slenderness ratio (L/r₀), and the soil moisture content. The elastic modulus of the concrete pile is 3 × 10⁷ kPa, and its radius is 0.5 m. This figure shows that $k_t$ generally decreases as the saturation degree increases. This means that a higher soil saturation degree leads to an increased settlement at the pile head under the same applied load. However, the rate of change in $k_t$ with suction diminishes at a certain limiting value of low saturation percentage for each soil type. The desaturation has a limited range, affecting the mechanical properties of sandy soils compared with those of fine-grained soils. For the SP soil considered in the present study, the pile-soil stiffness is not improved by further desaturation below $S_r=91\%$. Figure 14c shows that a decrease in saturation degree from 100% to 96% for CL soil leads to an improvement of pile-soil stiffness of about 44%, and the change in saturation from 100% to 77% results in an increase of pile-soil stiffness of 230%.

Figure 14e shows that for SP soil, the effect of a decrease in S_r from 100% to 91% on the pile-soil stiffness is two times the effect of a decrease in S_r from 91% to 25%. For CL soil, the change of S_r from 100% to 97% on the pile-soil stiffness is the same as the change from 97% to 87%. These effects of soil unsaturation on pile-soil system stiffness and, consequently, pile capacity may be different for fine-grained and coarse-grained soils. Thota and Vahedifard [67] observed similar effects on the shaft capacity of piles installed in unsaturated silt and clay. They intended to investigate the effect of temperature variation and, consequently, suction change on pile capacity. For this purpose, they entered the temperature-dependent matric suction into the relationship for pile capacity determination in unsaturated soil proposed by Vanapalli and Taylan [2]. According to the findings of Thota and Vahedifard [67], for a hypothetical 10 m pile, the pile capacity increased with soil drought in a limited range of matric suction in coarse-grained soils, but this capacity growth had a more even trend and larger range in fine-grained soils. Hence, it should be noted that in practice, under service loading, pile-soil system stiffness should be evaluated carefully considering different potential values of moisture contents.

3.2. Critical Slenderness Ratio for Pile Design

Figure 14 shows that the pile-soil system stiffness increases as the slenderness ratio increases but at a slower rate for larger slenderness ratios. The sensitivity analysis demonstrates that this trend is controlled by the pile-soil stiffness ratio, λ = $E_p/G_s$, in Equation (14). Since the shear modulus variation contributes to both shaft and pile base responses, Figure 14 shows that for greater soil shear modulus values, longer piles (with a larger slenderness ratio) generally experience greater settlement due to their elastic shortening. However, this observation was not consistent for all λ values in this dataset. In weaker soils in the present dataset, like SP and SM soils (Figure 14e,f), the $k_t$ value increases with increasing the pile slenderness ratio. This is because of the relatively negligible pile compressibility and weak base resistance. To explore this effect further, the pile-soil system stiffness was studied using Equation (14) for a fully unsaturated situation (n = 1) for an assumed pile with a 0.5 m radius, an elastic modulus of 200 GPa, and varying slenderness ratios. The soil shear modulus was chosen in a manner to change the λ value within the desired range. Figure 15 shows the behavior of the assumed pile in soil with a Poisson’s ratio of 0.3.

It is observed from Figure 15 that the slope of the curve $k_t$ against L/r₀, decreases after it reaches a peak value. The critical slenderness ratio (L/r₀)_cr for which the pile stiffness reaches its maximum may be estimated from:

$$({\frac{L}{r_0})}_{cr}=-3e^{-5}{\lambda }^2+0.124\lambda +14.4 \mathrm{for}\lambda <1000$$

(19)

Figure 15 shows that for λ > 1000, the (L/r₀)_cr values fall within an impractical range, and $k_t$ are monotonically ascending. Accordingly, choosing longer piles with a constant radius does not help to improve the pile-soil system stiffness and settlement control after a critical slenderness ratio value. For example, for a pile with Young’s modulus of 30 GPa installed in soil with a shear modulus of 60 MPa, the λ value is 500. Therefore, according to Figure 15, lengthening the pile after a L/r₀ of around 70 will no longer improve the pile-soil system stiffness. In this example, for a pile having a diameter of 0.5 m, considering a pile longer than 35 m for a smaller settlement at service loads will not lead to the desired result.

4. Estimation of Pile Head Settlement Corresponding to Ultimate Load in Unsaturated Soil

To calculate the pile settlement corresponding to the axial pile capacity in unsaturated soil, the pile load-settlement relationship should be adjusted by incorporating the average ultimate shear stress acting on the pile. For this purpose, the last term in the denominator of Equation (14) is evaluated for a pile with L/r₀ = 40, considering the pile is embedded in four soil types with various saturation degrees. The results are shown in Figure 16. As observed, the value of this term increases with decreasing S_r. However, this variation shows a small value ranging from 0.05 to 0.2 and thus may be ignored. Considering the fully unsaturated condition as proposed, the denominator in Equation (14) will be equal to unity.

To compensate for neglecting the last term in the denominator of Equation (14), the first term in the numerator, $4/\eta (1-\upsilon )$, representing the pile base contribution is ignored and thus Equation (20) for a fully unsaturated condition is obtained. As the soil becomes stiffer due to desaturation, Poisson’s ratio decreases and the term $4/\eta (1-\upsilon )$ is lower compared with the saturated condition. Overall, this simplification leads to overestimating the pile settlement in saturated soil by about 10%, and this value was found to be around 15% for unsaturated soils in the present study. Figure 17 shows the difference between P_t/Gw_tr₀ calculated from Equations (14) and (20) for ML soil using the data given in Figure 14b.

$$\frac{P_t}{G{r}_0{w}_t}=\frac{2\pi }{\xi }\rho \frac{L}{r_0}\frac{\tan \mathrm{h}\left(\mu L\right)}{\mu L}$$

(20)

Equation (21) is modified such that the shear stress is incorporated as:

$$\frac{w_t}{r_0}={\left(\frac{{\tau }_0^{\prime }}{G}\right)}_{ave}\frac{\xi \mu L}{\tan \mathrm{h}\left(\mu L\right)}$$

(21)

where ${\tau }_0^{\prime }$ denotes the average shear stress along the pile shaft and can be substituted with the soil shear strength for the appropriate pore water pressure dissipation condition. It should be noted that, as this equation uses average shear stress and does not consider pile driving forces, this equation is better used for bored piles or by having accurate characteristics of soil around the pile.

The validity of the above approach for unsaturated soil is evaluated in this section. As mentioned above, Equation (21) does not consider the soil resistance at the pile base. In tests listed in

Table 3. Soil and pile characteristics in pile load tests.

No.	Suction (kPa)	ω (%)	E (kPa)	ν	C (kPa)	Φ	Estimated Settlement (mm) from Equation (21)	Measured Settlement (mm) [2]	Ultimate Elastic Load (kN)
1	0	31	2500 (U)	0.49	11.5 (U)	0	0.31	0.40	0.10
2	0	31	7000 (D)	0.37	9 (D)	23	0.28	0.27	0.16
3	55	18	7500 (U)	0.49	58 (U)	0	0.50	0.40	0.48
4	55	18	17,500 (D)	0.37	21 (D)	23	0.22	0.30	0.62
5	110	16	8000 (U)	0.49	84 (U)	0	0.53	0.24	0.57
6	110	16	18,500 (D)	0.37	21 (D)	23	0.20	0.32	0.80
7	205	13	15,000 (U)	0.49	46 (U)	0	0.34	0.30	0.63
8	205	13	33,000 (D)	0.37	78 (D)	23	0.23	0.25	0.83

Note: (D) and (U) imply drained and undrained conditions, respectively.

, Vanapalli and Taylan [2] considered a 20-mm gap under the pile tip to eliminate the base resistance contribution. Hence, these tests can be used confidently to examine Equation (21) for unsaturated soil. Vanapalli and Taylan [2] used piles with a diameter of 2 cm and a length of 20 cm. Moreover, they installed the piles in unsaturated lateritic soil and glacial till, classified as low-plasticity clay (CL) soil by the Unified Soil Classification System (USCS), in a borehole. So, it is considered a bored pile, and the tests can be verified with the analytical solution. They obtained Young’s modulus from drained tests performed using ring shear tests and conducted unconfined compression tests to obtain Young’s modulus for undrained conditions. Young’s modulus for test No. 9 was estimated from the ratio of $E_{drained}/E_{undrained}$ of other tests [68]. The unsaturated shear modulus was calculated using Equation (16).

For undrained tests 1, 3, 5, and 7 listed in Table 3, the average shear stress in Equation (21) is substituted using the modified α method presented by Vanapalli and Taylan [2]. For unsaturated soil, this gives:

$$Q_{f(us)}=\alpha c_u{}_{(sat)}\left(1+\frac{(u_a-u_w)}{(P_a/101.3)}(S^{\nu })/\mu \right)\pi dL$$

(22)

For tests 2, 4, 6, and 8 in Table 3, a loading rate of 0.012 mm/min was applied to simulate the drained condition [2]. The modified β method based on Equation (23) presented by Vanapalli and Taylan [2] was used to examine Equation (21) for the drained condition. The modified β method consists of the conventional shear strength term and a term related to the shear strength due to the suction effect, which is given by:

$$Q_{f(us)}=\left[c_a^{\prime }+\beta ({\sigma }_z^{\prime })+(u_a-u_w)(S^k)(\tan {\delta }^{\prime })\right]\pi dL$$

(23)

where $c_a^{\prime }$ is the interface adhesion, β is the Burland–Bjerrum coefficient equal to k₀ tan δ, S is the soil saturation degree, δ is the pile-soil interface friction angle, $(u_a-u_w)$ is soil suction, and κ is a fitting parameter equal to $-0.0016I_p^2+0.0975I_p+1$ and I_p and represents the plasticity index.

Table 3 compares measured and computed pile head settlements with a mean relative error of 28%. The involved parameters in the α, β, or λ methods and Equation (21) and ensuring the drainage condition are important factors in obtaining proper results from the presented approach. A deviation of about 20% in suction measurement leads to a 9% error in the pile head settlement. Also, a 20% underestimation of the shear modulus causes an overestimation of 25% in the pile head settlement.

The variation of pile settlement/radius ratio, w_t/r₀, is obtained from Equation (21) for various slenderness ratios, and G/τ₀ is plotted in Figure 18 for μL = 0.8. It is noted from Figure 18 that the w_t/r₀ value decreases with increasing G/τ₀ This decrease has the same gradient for various slenderness ratios. This means that the variation of w_t/r₀ with G/τ₀ is independent of the pile slenderness ratio. Moreover, for large G/τ₀ values, there is no notable difference between w_t/r₀ of piles with $L/r_0$ greater than 40.

Figure 18 shows that for a pile with a constant radius and decreasing G/τ₀ with depth, the value of w_t increases, resulting in a larger pile settlement. It is also observed that, for a constant value of G/τ₀ and pile radius, as expected, short piles would experience less settlement compared to long ones. The shear stress in unsaturated soil depends on suction, and its variation with depth is totally dependent on environmental conditions. This should be considered carefully when increasing the pile length in the foundation design.

5. Implementation of the Proposed Method in Practice

As mentioned earlier, in many cases, the maximum saturation degree of the soil along the pile shaft is below a certain value for the entire foundation life cycle. Thus, designing foundations in saturated soil in such regions, with all the over-design and waste of resources it causes, is not in any way environmentally sustainable.

As there are no standard regulations for designing piles in unsaturated soils yet, the methods proposed in this study can be used based on engineering judgment. Engineers can study the behavior of piles with a water table along their shaft using Equation (14). However, for practical pile design, it is not safe to implement unsaturated soil mechanics when there is an active water table along the pile shaft. For arid regions with deep water tables, the maximum degree of saturation of soil in the intended life cycle can be considered. Then, the corresponding shear modulus can be obtained using Equation (15) (if the soil shear modulus for the dry and wet conditions is available) or common laboratory tests for obtaining shear modulus, and then by using Equation (15).

Having the unsaturated shear modulus and soil Poisson’s ratio, the pile head’s elastic settlement can be estimated using Equation (14) by substituting proper soil properties in both layers regarding the soil profile. Pile head settlement corresponding to the ultimate load in unsaturated soil can be calculated by Equation (21), using the unsaturated shear modulus and shear strength properties required in Equations (22) and (23). The test procedure for obtaining unsaturated shear strength properties is described in the respective references.

In many cases, like the installation of piles in unsaturated embankments, rehabilitation, and ground improvement projects, or when existing piles encounter special moisture conditions, having a detailed knowledge of pile behavior in different types of unsaturated soil is necessary. The findings described in sections “Parametric studies” and “Estimation of pile head settlement corresponding to ultimate load in unsaturated soil” provide suitable data in this regard.

6. Conclusions

An analytical solution was developed to evaluate the effect of unsaturated soil on the pile axial load-settlement behavior. The pile was assumed to be embedded in homogeneous and two-layered soil profiles, considering the soil unsaturation effect. A comprehensive contribution to the current research was to collect and analyze developed equations for unsaturated shear modulus and Poisson’s ratio in the literature to help with the evaluation of pile elastic settlement in this medium. The developed solution was then extended to approximate the ultimate soil-pile shear resistance and ultimate load carried by the pile at the desired settlement. Data from pile load tests in unsaturated soil available in the literature were used to examine the prediction accuracy of the developed model. The solution predictions were found to be in reasonable agreement with the experimental results.

Moreover, parametric studies on the effect of suction on pile-soil system stiffness were performed for six soil types, and critical slenderness was introduced, where the system behavior in unsaturated soil stops improving with the increment of pile slenderness. The following conclusions may be mentioned based on the findings of this research:

• In general, soil desaturation has a positive effect on soil-pile system stiffness, leading to improved pile performance;
• The effect of soil desaturation on pile-soil stiffness is significant over a wide range of saturation degrees for fine-grained soils, while for coarse-grained soils, the effective range is limited to saturation degrees close to 100%;
• The pile slenderness ratio has a small or no effect on the soil-pile stiffness at lower saturation degrees after a certain slenderness ratio, here called the critical slenderness ratio, which can be calculated using Equation (19). Thus, in unsaturated soils, increasing the pile length with a constant radius may not achieve better load-carrying characteristics above the critical L/r₀ value;
• The variation of w_t/r₀ with G/τ₀ for various L/r₀ values shows that in unsaturated soils for constant values of G/τ₀ and pile radius, shorter piles experience less settlement.

Overall, this paper contributes to clarifying and justifying the use of unsaturated soil mechanics in foundation engineering, especially pile foundations, and characterizes the load-settlement behavior of piles embedded in unsaturated soils. In situations where deep foundation settlement is critical, the presented analytical solution can be used for design. Furthermore, determining the load-carrying characteristics of pile groups by extending the solution is straightforward and can be addressed in further studies.