Load-Settlement Analysis of Axially Loaded Piles in Unsaturated Soils

: 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.


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 noncohesive 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.
Water 2024, 16 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.

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]: where τ xy is shear stress on x-plane in y direction, σ y is the total stress in y direction, u w is pressure in the water phase, τ zy is shear stress on z-plane in y direction, γ is the total unit weight of saturated soil, n w is porosity with respect to the water phase, n w γ 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 γ 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.
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.

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]: where  is shear stress on x-plane in y direction,  is the total stress in y direction,  is pressure in the water phase,  is shear stress on z-plane in y direction,  is the total unit weight of saturated soil,  is porosity with respect to the water phase,   is gravity body force for the water phase,  is porosity concerning the air phase (i.e., percentage of the surface of the element going through air), and   is gravity body force for the air phase,  is the percentage of element surface that is soil particles,  is the percentage of element surface that is contractile skin,  is the interaction force between the water and the contractile skin,  is the interaction force between the water and the soil particles,  is the interaction force between the air and the contractile skin,  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 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.

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.
nal elements of the stress state matrix [22].Although they influence value, suction invariants are not included in the linear shear stress-str unsaturated soils, and the relationship is the same as for saturated soil modulus and Poisson's ratio should be obtained for unsaturated soil co data or empirical expressions; this will be described subsequently.B mentioned assumptions and process, pile-soil system stiffness is eva form for the water table at an assumed level along the pile shaft.This w the next section.

Development of the Proposed Model
In practice, a part of the pile length may be located below the wate 2 presents a schematic view of a pile embedded in an unsaturated so saturated, deep soil deposit.The upper layer is subjected to the capilla the soil profile can be assumed to be two layers in terms of suction p there exist saturated and unsaturated zones.In this case, the variation level and its effect on the pile behavior should be investigated to avoid conservative design.In this section, a load-settlement equation for the two soil layers is developed using the above-mentioned approach.Considering the soil-pile system shown in Figure 2, for a pile with table level at nL (0 < n < 1) from the ground surface, the governing d for the pile settlement in each soil layer can be written as [12]: where w is pile settlement, z represents depth, ro is pile radius,  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]: where w is pile settlement, z represents depth, r 0 is pile radius, µ = , ν is soil Poisson's ratio, G is soil shear modulus, L is pile length, E p is pile material elasticity modulus, 4  η(1−υ) 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).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: 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: The term P t refers to the pile top load.For a compressible pile, the elasticity theory gives: 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: 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: By integrating Equation ( 9), the applied load at the pile head (z = 0) is given by: Assuming a constant soil Poisson's ratio along depth and G θ G s = α, the total pile-soil system stiffness (P t /w t ) is given by: Substituting n = 0 and α = 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 α = 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%.
Water 2024, 16, x FOR PEER REVIEW 7 o As shown in Figure 3, the pile settlement changes from 45% to 238% under cons axial load during the drying and wetting cycles.Such changes fall into a certain range a the first cycle, and they become about 33% larger than the saturated condition in the wetting cycle.This should be considered while studying the response of axially loa piles embedded in unsaturated soil, for which the effect of drying and wetting cycle the soil-pile interaction becomes of concern.When drying and wetting cycles are exp enced, the overall pile-soil stiffness range is altered by this interaction.To establish representative shear modulus and Poisson's ratio for the unsatura soil condition, the authors reviewed some correlations of shear modulus and Poiss ratio variation with the suction stress available in the literature, as shown in Table 1 seen in some studies, unsaturated shear modulus was correlated with suction by con 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.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.

Nomenclature
G unsat -unsaturated small strain shear modulus µ unsat -unsaturated Poisson's ratio G dry -small strain shear modulus at dry condition µ sat -unsaturated Poisson's ratio G sat -saturated small strain shear modulus µ d -unsaturated Poisson's ratio G 0 = the shear modulus at 1 kPa confining pressure S r -the soil saturation degree E unsat -unsaturated Young's modulus V s -shear wave velocity in the soil E sat -saturated Young's modulus G s -specific gravity OCR-over consolidation ratio K 0 -coefficient of lateral earth pressure at rest e-void ratio ρ-soil density (u a − u w )-soil suction P atm -atmospheric pressure (σ − u a )-principal net stress p r -the reference pressure (σ ′ − u a )-mean effective stress f (e) = 1 (0.3+0.7e 2 ) -the void ratio function proposed by Hardin and Black [28].(σ ′ − u a ) c -mean apparent reconsolidation stress

Validated for SM-ML soil
Cho and Santamarina [31] [ is the equivalent effective stress due to capillarity at a given S r

Proposed equation Description
Manusco et al. [32] (Gunsat 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.
β 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 G sat 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 A, n and κ are fitting parameters Fitting parameters were evaluated for 9 cohesive and cohesionless soils.

Biglari et al. [34]
-A 0 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 n 0 is a stiffness coefficient that accounts for the effect of P ′′ on stiffness in small strain range -OCR p ′′ 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(S r ) is a non-dimensional function of the degree of saturation defined above, where a ′ and b ′ are fitting parameters

Khosravi and McCartney [35]
λ p is the slope of the water retention curve -A g , n g , m g , k g and γ are fitting parameters Fitting parameters evaluated for completely decomposed tuff

Proposed equation Description
Zhou et al. [41] Tested on silty fine sand 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: where subscripts d and w stand for dry and wet conditions, respectively, θ 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 θ 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 2 /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 2 /gr and 2.04 m 2 /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. 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.
Poorly graded sand SP 6.8 0.602 17.2 0.584 34.5 0.565 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 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.

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 23 .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].
keep it consistent with the Mohr-coulomb model, the soil above groundwater was div into different layers, and the corresponding unsaturated soil elasticity modulus and sion to each suction level was used for them.A similar approach to modeling unsatu 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 c fied as low plasticity clay (CL) according to the Unified Soil Classification System (U The soil shear modulus changes with suction, as presented in Figure 5, as reporte Sawangsuriya et al. [25].The corresponding elastic modulus for each suction level was obtained from Figure 5, using Equation ( 16) and ν = 0.3.
For this soil, the saturated cohesion and the friction angle were 10 kN/m 2 and ϕ = 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: 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.: where κ = −0.0016I 2 p + 0.0975I p + 1 and δ ′ is the interface friction angle, which is assumed to be 0.8ϕ 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.
step, the pile material is substituted in the hole, and the model is analyzed under g loading.Finally, the prescribed displacement is applied to the pile head.7 shows vertical displacement vectors for the FE analysis of a pile emb in the ground with the water table level at −3 m from the ground surface.Figures 8 present vertical displacement sections for the soil around the pile shaft and below th base at the groundwater table levels of −3 and −10 m.As seen, when the groundwate is at −10 m, the soil around the pile shaft is fully mobilized, and the displacement proximately uniform along the pile shaft.This is due to the large shear stress in thi dition compared with the state when the water table level is −3 m.Moreover, for th below the pile base, it takes a longer distance for the vertical displacement to reac because of the large load concentration.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.Figures 8 and 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.
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 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 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 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.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.Variation of pile-soil system stiffness with nL using Equation ( 14) for cases in Figure 11.

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 highplasticity 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  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.Figure 12.Variation of pile-soil system stiffness with nL using Equation ( 14) for cases in Figure 11.

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 highplasticity 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 Variation of pile-soil system stiffness with nL using Equation ( 14) for cases in Figure 11.

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.
Water 2024, 16, 337 20 of 29 pores resists more due to high surface tension and strong molecular As dr proceeds, the suction stress increment decreases for coarse-grained soils because of t larger pore sizes and sparser pore water compared with fine-grained ones.On the o hand, the suction stress increment is more uniform in fine-grained soils, which can to a wider range of shear modulus improvements.The data for soils were obtained from tests conducted employing an apparatus could apply two stress state variables, the net confining pressure and matric suction, ing 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 ( ), the slenderness ratio (L/ro), and the soil moisture content.The elastic modulus of the conc pile is 3 × 10 7 kPa, and its radius is 0.5 m.This figure shows that  generally decre as the saturation degree increases.This means that a higher soil saturation degree l to an increased settlement at the pile head under the same applied load.However, the of change in  with suction diminishes at a certain limiting value of low saturation centage for each soil type.The desaturation has a limited range, affecting the mechan properties of sandy soils compared with those of fine-grained soils.For the SP soil con ered in the present study, the pile-soil stiffness is not improved by further desatura below  = 91%.Figure 14c shows that a decrease in saturation degree from 100% to for CL soil leads to an improvement of pile-soil stiffness of about 44%, and the chang saturation from 100% to 77% results in an increase of pile-soil stiffness of 230%.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 0 ), and the soil moisture content.The elastic modulus of the concrete pile is 3 × 10 7 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.

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.
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.

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, λ =  / , 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  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  against L/ro, decreases after it reaches a peak value.The critical slenderness ratio (L/ro)cr for which the pile stiffness reaches its maximum may be estimated from: It is observed from Figure 15 that the slope of the curve k t against L/r 0 , decreases after it reaches a peak value.The critical slenderness ratio (L/r 0 ) cr for which the pile stiffness reaches its maximum may be estimated from: L r 0 cr = −3e −5 λ 2 + 0.124λ + 14.4 for λ < 1000 (19) Figure 15 shows that for λ > 1000, the (L/r 0 ) 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 0 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.

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 0 = 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.

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/ro = 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 Sr.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/(1 − ), 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/(1 − ) 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 Pt/Gwtro calculated from Equations ( 14) and (20) for ML soil using the data given in Figure 14b.To compensate for neglecting the last term in the denominator of Equation ( 14), the first term in the numerator, 4/η(1 − υ), 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/η(1 − υ) 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 t r 0 calculated from Equations ( 14) and (20) for ML soil using the data given in Figure 14b.Equation ( 21) is modified such that the shear stress is incorporated as: where  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.Equation ( 21) is modified such that the shear stress is incorporated as: where τ ′ 0 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, 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: 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: where c ′ a is the interface adhesion, β is the Burland-Bjerrum coefficient equal to k 0 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 2 p + 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 0 , is obtained from Equation ( 21) for various slenderness ratios, and G/τ 0 is plotted in Figure 18 for µL = 0.8.It is noted from Figure 18 that the w t /r 0 value decreases with increasing G/τ 0 This decrease has the same gradient for various slenderness ratios.This means that the variation of w t /r 0 with G/τ 0 is independent of the pile slenderness ratio.Moreover, for large G/τ 0 values, there is no notable difference between w t /r 0 of piles with L/r 0 greater than 40. Figure 18 shows that for a pile with a constant radius and decreasing G/τ0 with depth, the value of wt increases, resulting in a larger pile settlement.It is also observed that, for a constant value of G/τ0 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.

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 Figure 18 shows that for a pile with a constant radius and decreasing G/τ 0 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/τ 0 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.

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.

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 0 value; • The variation of w t /r 0 with G/τ 0 for various L/r 0 values shows that in unsaturated soils for constant values of G/τ 0 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.

Figure 2 .
Figure 2. Pile embedded in the two-layered soil profiles.

Figure 2 .
Figure 2. Pile embedded in the two-layered soil profiles.

Figure 3 .
Figure3.Elastic pile settlement under dying and wetting path calculated using Equation (14) corresponding saturation degrees.

Figure 3 .
Figure3.Elastic pile settlement under dying and wetting path calculated using Equation (14) and corresponding saturation degrees.
reflecting the inherent soil properties in the ij plane, with dimensions m/s -(σ i − u a ) and σ 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] (G o, us ) hh = A f (e)(σ 0 − u a ) n + Cθ κ (u a − u w ) -A, C and κ are fitting parameters -

-
A and n are experimental parameters p ′ c 0= initial mean apparent pre-consolidation stress -∆e p = plastic change in void ratio λ and κ the slopes of the loading and unloading branches of the isotropic compression curve -K and K ′ = 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 soilOh and Vanapalli[36]

θ w −θ d m -
w and d indices are related to the wet and dry states of soil θ = volumetric water content m is an experimental fitting parameterFitting parameter evaluated for 16 cohesive and cohesionless soils Georgetti et al.[37] an experimental fitting parameter having the same units as the small-strain shear modulus -S e = effective degree of saturation γ 0 and β are experimental fitting parameters Fitting parameters were evaluated for 29 cohesive and cohesionless soils.Cao et al. [44] V s = (a−e) 2 1+e (σ c −u a )−(u w −u a )( uw −ua uw −u a0 ) −0.55 1kPa n G uns = ρV s a and n are fitting parameters Fitting parameters evaluated for copper tailing Liu et al. [45] G uns = AF(e)F(w c ) σ ′ P a n F(w c ) = (100w c ) −α w c is water content in decimal -A, α and n are fitting parameters Fitting parameters evaluated for Yan'an loess Khosravi et al. [46]

-
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 ≥ 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

Figure 5 .
Figure 5. Variation of shear modulus with degree of saturation for CL soil[25].

Figure 6 .
Figure 6.Generated fine mesh for the model with the water table at −6 m.The pile under pre fined displacement is shown on the left and the colors show different saturation levels in the

Figure
Figure7shows vertical displacement vectors for the FE analysis of a pile emb in the ground with the water table level at −3 m from the ground surface.Figures 8 present vertical displacement sections for the soil around the pile shaft and below th base at the groundwater table levels of −3 and −10 m.As seen, when the groundwate is at −10 m, the soil around the pile shaft is fully mobilized, and the displacement proximately uniform along the pile shaft.This is due to the large shear stress in thi dition compared with the state when the water table level is −3 m.Moreover, for th below the pile base, it takes a longer distance for the vertical displacement to reac because of the large load concentration.

Figure 7 .
Figure 7. Vertical displacements of soil for water table level at −3 m from the ground surface

Figure 6 .
Figure 6.Generated fine mesh for the model with the water table at −6 m.The pile under predefined displacement is shown on the left and the colors show different saturation levels in the soil.

Figure 7 .
Figure 7. Vertical displacements of soil for water table level at −3 m from the ground surface.Figure 7. Vertical displacements of soil for water table level at −3 m from the ground surface.

Figure 7 .Figure 8 .Figure 9 .
Figure 7. Vertical displacements of soil for water table level at −3 m from the ground surface.Figure 7. Vertical displacements of soil for water table level at −3 m from the ground surface.Water 2024, 16, x FOR PEER REVIEW 19 of 32

Figure 9 .
Figure 9. Vertical displacement of soil: (a) below pile base (left) and (b) soil along pile shaft (right) for water table at −10 from the ground surface.

Figure 10 .
Figure 10.Distribution of total stresses in soil corresponding to the pile head displacement of 2 cm and for groundwater table level at (a) −3 m, (b) −6 m, (c) −10 m, from the ground surface.

Figure 10 .
Figure 10.Distribution of total stresses in soil corresponding to the pile head displacement of 2 cm and for groundwater table level at (a) −3 m, (b) −6 m, (c) −10 m, from the ground surface.

Figure 12 .
Figure 12.Variation of pile-soil system stiffness with nL using Equation (14) for cases in Figure11.

Figure 11 .
Figure 11.Variation of load-pile head settlement results for numerical and analytical analyses for water table levels at −3, −6, and −10 m.

Figure 11 .
Figure 11.Variation of load-pile head settlement results for numerical and analytical analyses for water table levels at −3, −6, and −10 m.
Figure 12.Variation of pile-soil system stiffness with nL using Equation (14) for cases in Figure11.

Figure 14 .
Figure 14.Variation of pile-soil system stiffness with slenderness ratio (L/ro) and saturation degree for different soil types.

Figure 15 .
Figure 15.Variation of pile-soil system stiffness with slenderness ratio (L/ro) for different λ values.

Figure 15 .
Figure 15.Variation of pile-soil system stiffness with slenderness ratio (L/r 0 ) for different λ values.

Figure 16 .
Figure16.The value of the last term in the denominator of Equation (14) for different soil types with various saturation degrees.

Figure 16 .
Figure16.The value of the last term in the denominator of Equation (14) for different soil types with various saturation degrees.
Water 2024, 16, x FOR PEER REVIEW 26 of 32

Figure 17 .
Figure 17.Variation of load-settlement ratio versus Sr, calculated from Equations (14) and (20) for L/ro = 40 and ML soil data given in Figure14b.

Figure 17 .
Figure 17.Variation of load-settlement ratio versus S r , calculated from Equations (14) and (20) for L/r 0 = 40 and ML soil data given in Figure14b.

Table 1 .
Shear modulus and Poison's ratio relations with matric suction.
and ζ are fitting parameters 5 -T s 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 ν − ln(σ ′ − u a ) plane C 0 and n s are experimental parameters [35]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
ω i and 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 unsat = µ dry + µ sat − µ dry 1 − 1 + (ua−uw)

Table 2 .
"m" value for different soil types.

Table 3 .
Soil and pile characteristics in pile load tests.
Note: (D) and (U) imply drained and undrained conditions, respectively.