A Hyperelastic Bounding Surface Plasticity Model for Unsaturated Granular Soils

: In this paper, a state-dependent, bounding surface plasticity model that simulates the behavior of unsaturated granular soils is presented. An unsaturated, soil mechanics-compatible elastoplastic response is adopted in which no part of the response occurs in a purely elastic fashion. To create an appropriate hydro-mechanical coupling, a newer generation stress framework, consisting of the Bishop-type effective stress and a second stress variable, is used in conjunction with a soil-water characteristic curve function. Details regarding the model development, parameter estimation


Introduction
Environmental conditions such as precipitation, infiltration, and evaporation influence the location of the groundwater table and, therefore, the thickness of the unsaturated soil zone.In this zone, the soil is characterized by variations in moisture content and in the associated water pressure.Ignoring conditions in the tightly adsorbed water layer in the unsaturated zone, the pore-water pressure is negative.Depending on the type and permeability of the soil, the magnitude of the pore pressure can vary over six orders of magnitude.Reports indicate that, in the USA alone, problems associated with unsaturated soils (e.g., volume change, pressure generation, and moisture transport) annually inflict over $15 billion in damage to buildings, roads, and pipelines, e.g., [1].
Constitutive modeling of unsaturated soils has drawn attention since the 1990s.Research in this field was pioneered by the likes of Alonso, et al. [2], who extended the modified Cam-Clay model [3] to one based on the concept of two independent stress state variables.In this concept, the behavior of unsaturated soils is governed by matric suction and net stress [4,5].Based on the concept of two stress state variables and using the extended critical state concept, several constitutive models were subsequently developed and used to simulate the behavior of unsaturated soils, e.g., [6][7][8][9][10].However, Gens, et al. [11] postulated that the framework of two independent stress state variables is successful only in cases where one of these variables is kept constant.According to Nuth and Laloui [12], separating the mechanical stress from the hydraulic stress prevents the proper simulation of the hydro-mechanical coupling phenomena observed in unsaturated soils.Sheng [13] indicated that, in the two-stress state variable approach, a smooth transition between saturated and unsaturated states was missing.Using energy principles, Houlsby [14] rigorously demonstrated that two sets of stress-strain variables (e.g., the effective stress and matric suction as stress variables and the strain of the solid skeleton and degree of saturation as kinematic variables) are adequate to describe unsaturated soil behavior.
The Bishop-type [15] effective stress has been extensively used in describing unsaturated soil behavior, as evident in the works of [16][17][18][19][20][21][22][23][24][25][26][27][28][29][30].The constitutive relationships presented in the aforementioned studies fall within the framework of standard rate-independent elastoplasticity.They predict path-dependent elastic response for stress states within a yield surface, which is counter to the framework of elastoplasticity.In such formulations, the elastic response should be path-independent.Additional details regarding these observations will be provided in the Strain Decomposition section of this paper.
The effective stress framework allows the direct incorporation of the soil water characteristic curve (SWCC) and the intrinsic coupling between the mechanical and hydraulic behavior.In addition, when the state of the soil changes to full saturation, the associated effective stress reduces to Terzaghi's effective stress.This characteristic makes such a framework a viable choice for implementing a constitutive model for unsaturated soils into new or existing finite element or finite difference computer programs, as most of these are written in terms of effective stress for saturated soils.Consequently, a natural transition to unsaturated states can be realized by simply replacing Terzaghi's effective stress with Bishop's definition [12].
In elastoplastic models for unsaturated soils employing the effective stress framework, the yield surface can be made a function of the matric suction and/or the degree of saturation [19,24,26,31,32].This gives such models the ability to correctly simulate the influence of unsaturated conditions on plastic straining, including the possible occurrence of collapse compression during wetting.
The concept of a bounding surface in stress space was independently introduced by Krieg [33] and Dafalias and Popov [34] to more accurately simulate the cyclic response of metals.The basic idea underlying bounding surface models is a smooth transition from elastic to plastic states.In such models, there is a yield surface whose size and location may vary within an outer or "bounding" surface.The plastic stiffness decreases as the material state approaches the bounding surface.This beneficial characteristic has led to the development of several bounding surface models for saturated soils [35][36][37][38][39][40].
Russell and Khalili [32] first applied the bounding surface concept to simulate the rate-independent behavior of unsaturated soils.They combined the framework of critical state soil mechanics with Bishop's effective stress definition for such soils.Their model also accounted for possible particle crushing at high stresses.The capabilities of this model were assessed by simulating the axisymmetric triaxial response of both clean sand and pure clay.
By extending the bounding surface plasticity model of Bardet [36], which was originally developed for saturated sands, Morvan, et al. [41] proposed a bounding surface-based constitutive model for unsaturated soils.More recently, Patil, Hoyos, Morvan and Puppala [28] refined this model, taking into account the hysteretic behavior of a water retention curve and its dependence on the void ratio.A common trait of the aforementioned bounding surface constitutive models is that they all assume a hypoelastic response.However, as pointed out by Lashkari and Golchin [42], the lack of conservation of energy in hypoelasticity poses challenges in proving the accuracy of such models.
In the present work, a state-dependent, critical state-compatible constitutive model that accounts for the effect of matric suction on the hydro-mechanical response of unsaturated granular soils is described.The model is formulated in a rate-independent bounding surface plasticity framework and incorporates a newer generation of effective stress definition.The model represents an enhanced extension of the hyperelastic bounding surface model for fully saturated granular soils originally proposed by Lashkari,et al. [43].In the unsaturated, soil mechanics-compatible elastoplastic response associated with the present model, no part of the response occurs in a purely elastic fashion.The elastic portion of the elastoplastic response is computed from a Gibbs free energy function suitable for unsaturated granular soil, thus satisfying the law of conservation of energy.Using a single set of parameter values, this model simulates the hydro-mechanical response of unsaturated granular soils with various ranges of density, normal stress, and matric suction.The model's predictive capabilities are assessed by comparing its simulations with experimental data for clean sand, as reported by Russell and Khalili [32], and silty sand, as reported by Patil, Hoyos and Puppala [29].The model is also compared with the bounding surface-based models of Russell and Khalili [32], Morvan, Wong and Branque [41], and Patil, Hoyos, Morvan and Puppala [28].To rigorously assess the model's performance, the statistical technique of Li and Zhu [44] is employed, wherein the normalized residual sum of squares (RSS mono ) serves as the benchmark for comparison.
The authors deeply appreciate the importance of evaluating the model performance with varying stress paths, other test types, and plastic volumetric strains induced by wetting and drying cycles under isotropic loading conditions.However, due to page limit issues, the details pertaining to these issues must be deferred to future publications.

Model Development
The following sections present the development of the model.For simplicity, the model is formulated in axisymmetric triaxial space in terms of the mean normal effective stress (p') and the deviator stress (q).Although implementing the model into general research and commercial computer programs would necessitate a more general, three-invariant framework, such development is beyond the scope of the present paper.

Effective Stress Definition
Following the principles of saturated soil mechanics, the effective stress vector (σ ′ ) is written as follows: where p' and q are defined above.By adopting Bishop's relationship, p ′ is defined as follows: where p net = (p t − u a ) is the net confining stress, p t is the total confining stress, u a is the pore air pressure, s = (u a − u w ) is the matric suction, and u w is the pore water pressure.In this work, matric suction signifies the component of the total suction, excluding the osmotic fraction.In Equation (2), χ is Bishop's effective stress parameter, which varies between a value of zero (corresponding to dry conditions) and one (corresponding to water-saturated conditions).Several functional relationships have been proposed for χ, most of which relate it to the degree of saturation (S r ).Schefler [45] suggested that χ be equal to S r .The use of this definition, in conjunction with experimental data for unsaturated soils with different values of matric suction, showed that a unique critical state line (CSL) was, however, not achieved in p'-q space.This observation is evident in Figure 1a, where it can be seen that for each value of matric suction (s), a different critical state line must be defined.
To avoid the need for defining multiple CSLs in p'-q space, Manzanal, Pastor and Merodo [24] suggested that the effective degree of saturation, S re , be used in the definition of χ as follows: where S r is the current degree of saturation and S r0 is the residual state degree of saturation, which can be estimated from the SWCC.The symbols 〈 〉 represent Macaulay brackets; for some quantity x, they give ⟨x⟩ = x if x > 0 and ⟨x⟩ = 0 if x ≤ 0. As evident from Figure 1b, Equation (3) eliminates the need for defining multiple CSLs in p'-q space.

The Idealization of the Soil Water Characteristic Curve (SWCC)
The constitutive relationship between matric suction and the water content (or deg of saturation) of soils is commonly described using the SWCC.A number of functio forms have been suggested for SWCCs, (e.g., [47][48][49][50][51][52][53][54]).The SWCC of a given soil va depending on the direction and history of wetting considered.This variation in SWCC referred to as hysteresis [55]; its physical description was pioneered by Poulovassilis [ Past investigations have shown that the SWCC of natural soil differs from the SWCC the same soil recompacted to the same density at the same water content [57,58].Sim studies have reported that additional factors such as stress state, loading and wetting tory, and soil structure affect the SWCC.In addition, many studies have investigated dependence of the SWCC on a soil's initial density, (e.g., [6,32,[57][58][59][60][61][62][63]).For example, fr Figure 2, it is evident that different SWCCs are obtained for different densities, as qua fied by the void ratio (e).
Following the work of Wheeler [62] but adopting a simplified approach, Gallip Wheeler and Karstunen [63] proposed a normalization technique that yields a uni SWCC for different values of density.In this technique, the degree of saturation is co lated with not only the matric suction but also the density (i.e., specific volume).Adopt such a technique, Manzanal, Pastor and Merodo [24] defined a modified suction (s*) t

The Idealization of the Soil Water Characteristic Curve (SWCC)
The constitutive relationship between matric suction and the water content (or degree of saturation) of soils is commonly described using the SWCC.A number of functional forms have been suggested for SWCCs, (e.g., [47][48][49][50][51][52][53][54]).The SWCC of a given soil varies depending on the direction and history of wetting considered.This variation in SWCCs is referred to as hysteresis [55]; its physical description was pioneered by Poulovassilis [56].Past investigations have shown that the SWCC of natural soil differs from the SWCC of the same soil recompacted to the same density at the same water content [57,58].Similar studies have reported that additional factors such as stress state, loading and wetting history, and soil structure affect the SWCC.In addition, many studies have investigated the dependence of the SWCC on a soil's initial density, (e.g., [6,32,[57][58][59][60][61][62][63]).For example, from Figure 2, it is evident that different SWCCs are obtained for different densities, as quantified by the void ratio (e).
In the present model, the relationship of Fredlund and Xing [47] for the SWCC is adopted and combined with the aforementioned normalization technique of Gallipoli, Wheeler and Karstunen [63].The resulting normalized SWCC function thus takes the following form: In Equation (5),  ,  , and  are model parameters, and  and  are the residual values of the degree of saturation and matric suction, respectively.The value of 10 6 corresponds to the suction of 1 GPa, which is the value at which the degree of saturation,  , reaches the residual value,  , in the SWCC expression proposed by Fredlund and Xing [47].

Strain Decomposition
The strain vector () that is work-conjugate to the effective stress vector given in Equation ( 1) is as follows: where  and  are, respectively, the volumetric and distortional components of the strain vector.
Assuming infinitesimal displacements and displacement gradients, the total strain rate vector is additively decomposed into an elastic and plastic part as follows: where the superscripts e and p denote the elastic and plastic parts, respectively, of .A superposed dot indicates a material time derivative or rate, which is used here in lieu of incremental notation as a convenience since the formulation is rate-independent.Within a conventional elastoplasticity framework, the response associated with stress states within the elastic region is, by definition, independent of the stress path.However, Zhang and Lytton [64] and Sheng, et al. [65] showed that for unsaturated soils, the "elastic" volumetric change is stress-path dependent.As Figure 3 illustrates, a closed loop of net stress and suction changes within the elastic region does not necessarily lead to a Following the work of Wheeler [62] but adopting a simplified approach, Gallipoli, Wheeler and Karstunen [63] proposed a normalization technique that yields a unique SWCC for different values of density.In this technique, the degree of saturation is correlated with not only the matric suction but also the density (i.e., specific volume).Adopting such a technique, Manzanal, Pastor and Merodo [24] defined a modified suction (s*) that is a function of the void ratio according to the following: where Ω is a curve-fitting parameter.
In the present model, the relationship of Fredlund and Xing [47] for the SWCC is adopted and combined with the aforementioned normalization technique of Gallipoli, Wheeler and Karstunen [63].The resulting normalized SWCC function thus takes the following form: In Equation (5), a v , n v , and m v are model parameters, and S r0 and S 0 are the residual values of the degree of saturation and matric suction, respectively.The value of 10 6 corresponds to the suction of 1 GPa, which is the value at which the degree of saturation, S r , reaches the residual value, S r0 , in the SWCC expression proposed by Fredlund and Xing [47].

Strain Decomposition
The strain vector (ε) that is work-conjugate to the effective stress vector given in Equation ( 1) is as follows: where ε v and ε q are, respectively, the volumetric and distortional components of the strain vector.
Assuming infinitesimal displacements and displacement gradients, the total strain rate vector is additively decomposed into an elastic and plastic part as follows: where the superscripts e and p denote the elastic and plastic parts, respectively, of ε.A superposed dot indicates a material time derivative or rate, which is used here in lieu of incremental notation as a convenience since the formulation is rate-independent.Within a conventional elastoplasticity framework, the response associated with stress states within the elastic region is, by definition, independent of the stress path.However, Zhang and Lytton [64] and Sheng, et al. [65] showed that for unsaturated soils, the "elastic" volumetric change is stress-path dependent.As Figure 3 illustrates, a closed loop of net stress and suction changes within the elastic region does not necessarily lead to a closed loop of effective mean stress changes.This is because of the material state dependency of the effective stress [13].The lack of a closed-loop means that, even for seemingly elastic stress states, the response is stress-path dependent.A model that exhibits stresspath-dependent elastic behavior is inconsistent with classical elastoplasticity theory and thermodynamic considerations.
Geosciences 2024, 14, x FOR PEER REVIEW 6 of closed loop of effective mean stress changes.This is because of the material state depen ency of the effective stress [13].The lack of a closed-loop means that, even for seeming elastic stress states, the response is stress-path dependent.A model that exhibits stres path-dependent elastic behavior is inconsistent with classical elastoplasticity theory an thermodynamic considerations.According to Collins and Houlsby [66], for materials in which elastic parameters a dependent on internal parameters due to elastic-plastic coupling, the traditional elast plastic strain decomposition is not valid.Following the work of Maier and Hueckel [6 Collins and Houlsby [66] postulated that instead of using "elastic" and "plastic", the term "reversible" and "irreversible" should instead be used when describing strain rates.their definition, the elastic strain was assumed to be affected by the plastic strain; coupl elastic-plastic strain rates were thus added to the elastic strain rates.Collins and Houls [66] thus decomposed the elastic strain rate into reversible and irreversible (coupled) po tions.Table 1 summarizes their strain decomposition.In this table, the superscripts " and "i" denote the reversible and irreversible parts, respectively, of the infinitesimal stra vector.According to Collins and Houlsby [66], for materials in which elastic parameters are dependent on internal parameters due to elastic-plastic coupling, the traditional elasticplastic strain decomposition is not valid.Following the work of Maier and Hueckel [67], Collins and Houlsby [66] postulated that instead of using "elastic" and "plastic", the terms "reversible" and "irreversible" should instead be used when describing strain rates.In their definition, the elastic strain was assumed to be affected by the plastic strain; coupled elasticplastic strain rates were thus added to the elastic strain rates.Collins and Houlsby [66] thus decomposed the elastic strain rate into reversible and irreversible (coupled) portions.Table 1 summarizes their strain decomposition.In this table, the superscripts "r" and "i" denote the reversible and irreversible parts, respectively, of the infinitesimal strain vector.ε i (Adopted from Collins and Houlsby [66]).
In the present model, due to the aforementioned path dependence of the elastic strains in unsaturated soils (Figure 3), it is assumed that no purely elastic strains can develop at any time; the elastic strain is thus inherently coupled with the plastic strain.In light of the above discussion, the term "irreversible elastic strain vector" is adopted for the elastic strain that is coupled with the plastic strain.

Definition of the Elastic Response
The constitutive relation used by Golchin and Lashkari [68] and Alipour and Lashkari [69] is likewise adopted herein, i.e.: .
where K, G, and J are the hyperelastic moduli, which are non-linear functions of p'.Following the approach of Golchin and Lashkari [68], expressions for these moduli are calculated by partial differentiation of the Gibbs free energy function proposed by Einav and Puzrin [70], which guarantees the conservation of energy in any arbitrary closed loop.As derived by Golchin and Lashkari [68], this function is as follows: where K o and G o are model parameters, P re f is a reference pressure that is normally assumed to be equal to the atmosphere pressure (P atm = 101 kPa), p 0 and q 0 are the values of p' and q at zero elastic strain, respectively, and F(e) is a function of the void ratio (e) and particle shape.Lashkari, Karimi, Fakharian and Kaviani-Hamedani [43] proposed the following expression for granular soils with angular particles: The comparable expression for granular soils with well-rounded particles is as follows: In Equation ( 9), the elastic variable X, which evolves with plastic hardening, is defined as follows: where values of X max = 0.95, X min = 0.50, and B = 0.10 are adopted from Lashkari, Karimi, Fakharian and Kaviani-Hamedani [43].In Equation (12), M is the slope of the CSL in p'-q space that takes on values of M c and M e for axisymmetric triaxial compression and extension, respectively.The quantity η = q/p' is the stress ratio, n b is a model parameter (to be discussed in the next section), and ψ is the Been and Jefferies state parameter [71] that represents the state dependency of the proposed model by relating a given state of soil with its critical state according to the following: where e c is the value of the void ratio (e) at critical state for the current value of p'.
As mentioned above, the moduli K, G, and J appearing in Equation ( 8) are obtained by suitable partial differentiation of the Gibbs free energy function, giving the following For the special case of isotropic loading, η = 0.The above moduli then reduce to the following functional forms, originally proposed by Hardin and Richart [72]:

Definition of the Bounding Surface
In p'-q space, the bounding surface is represented by a straight line with the slope of M b (see Figure 4).M b is related to the slope of the CSL, M, according to the following: where n b is a model parameter by which the peak shear stress is simulated (also see Equation ( 12)).
Geosciences 2024, 14, x FOR PEER REVIEW 8 of 25 For the special case of isotropic loading,  = 0.The above moduli then reduce to the following functional forms, originally proposed by Hardin and Richart [72]:

Definition of the Bounding Surface
In p'-q space, the bounding surface is represented by a straight line with the slope of M b (see Figure 4).M b is related to the slope of the CSL, M, according to the following: where n b is a model parameter by which the peak shear stress is simulated (also see Equation ( 12)).

Definition of the Dilatancy Surface
Similar to the bounding surface, the dilatancy surface is defined as a straight line in p'-q space with a slope of M d (Figure 4).M d is also a function of the slope of the CSL according to the following: Geosciences 2024, 14, 148 9 of 25

Definition of the Dilatancy Surface
Similar to the bounding surface, the dilatancy surface is defined as a straight line in p'-q space with a slope of M d (Figure 4).M d is also a function of the slope of the CSL according to the following: where n d is a model parameter by which the volume phase transformation is simulated.

Definition of the Critical State Void Ratio
Li and Wang [73], Lashkari and Yaghtin [74], and Lashkari, Karimi, Fakharian and Kaviani-Hamedani [43] used the following relationship to define critical state void ratio for fully saturated soils: where e 0 , λ, and ζ are model parameters and p ref is as previously defined.One important attribute of Equation ( 22) is that, for fully saturated granular soils with various initial densities and stress states, it yields a unique CSL in the space of e versus For unsaturated states, critical state experimental data in this space are, however, dependent on the values of matric suction.Consequently, as shown in Figure 5a, a unique line cannot fit all of the experimental data [10].Indeed, different CSLs need to be defined in this space for different levels of matric suction.Referring to the theoretical expression derived by Fisher [76], the maximum value for () is 1.5 (for infinite suction), while its minimum value is 1.0 (for zero suction).In the present model, the following relationship, proposed by Lashkari and Kadivar [77], is used for (): Figure 6 shows the variation of () with suction s.
Recalling that Equation ( 22) applies for fully saturated conditions, in light of Equation (23), this expression is modified for unsaturated granular soils, giving the following: Gallipoli, et al. [75] proposed a normalization approach that brings critical state data with different values of matric suction into a unique line.Their approach, which is adopted in the present model, consists of the following normalization, which relates the critical state effective stress for saturated conditions to its value for unsaturated conditions according to the following: where p In Equation ( 23), the function g(ξ) is given by the following: where a and b are model parameters and ξ is a "cementation" parameter.The latter is defined by the following relationship [75]: where S r is again the degree of saturation.The term (1 − S r ) takes into account the number of water menisci per unit volume of the soil.The suction-dependent function f (s) represents the magnitude of the inter-particle contact force due to the matric suction.
Referring to the theoretical expression derived by Fisher [76], the maximum value for f (s) is 1.5 (for infinite suction), while its minimum value is 1.0 (for zero suction).In the present model, the following relationship, proposed by Lashkari and Kadivar [77], is used for f (s): Figure 6 shows the variation of f (s) with suction s.

Definition of the Flow Rule, Hardening Rule, and Dilatancy Function
The plastic shear strain rate is obtained using the non-associative flow rule of Dafalias and Manzari [39], which is written as follows: where  is a scalar loading index, η is again the stress ratio, and  is the dimensionless plastic-hardening modulus given by the following: where ℎ and  are dimensionless model parameters,  is given by Equation (20), and  is the elastic variable defined in Equation ( 12).The total irreversible strain vector is then calculated from the following equation: Recalling that Equation ( 22) applies for fully saturated conditions, in light of Equation ( 23), this expression is modified for unsaturated granular soils, giving the following: where, e 0 , λ, and ζ are again model parameters.Unlike Equation (22), Equation (27) provides a unique CSL in e − p ′ p re f (1+g(ξ)) ζ space for unsaturated soils at different initial conditions.Figure 5b illustrates the normalization effect of Equation (27).It is evident that this normalization transforms the scattered experimental data shown in Figure 5a more nearly into a unique line.

Definition of the Flow Rule, Hardening Rule, and Dilatancy Function
The plastic shear strain rate is obtained using the non-associative flow rule of Dafalias and Manzari [39], which is written as follows: .
where Λ is a scalar loading index, η is again the stress ratio, and K p is the dimensionless plastic-hardening modulus given by the following: where h 0 and c h are dimensionless model parameters, M b is given by Equation (20), and X is the elastic variable defined in Equation ( 12).The total irreversible strain vector is then calculated from the following equation: where ∂q∂X , and X ,η = ∂X ∂η , with Γ being the Gibbs free energy function given by Equation (9).
For isotropic loading (J = 0), by imposing the consistency condition and combining it with Equations ( 8) and ( 30), the scalar loading index is calculated from the following equation: which has the form of the loading index originally given by Dafalias and Manzari [39].In Equation (31), d is the following dilatancy function: In Equation (32), A d is a model parameter, and c = M e /M c , where M c and M e are again the slopes of critical state lines in axisymmetric triaxial compression and extension, respectively [39].According to Lashkari, Karimi, Fakharian and Kaviani-Hamedani [43], for practical problems, a value of c = 0.7 provides reasonable results.As such, this value is adopted in this work.Knowing that d = .
and considering Equation (30), the plastic volumetric strain rate is given by the following: .

Determination of Model Parameters
Associated with the most general form of the present model are 14 parameters (K 0 , G 0 , e 0 , λ, a, b, ζ, M c , M e , n b , n d , A d , h 0 , and c h ) whose values must be determined for a particular soil.There are also four parameters with default values of p ref = 101 kPa, X max = 0.95, X min = 0.50, and B = 0.10 that can be used for all soils.Six additional parameters (a v , Ω, n v , m v , S r0 , and s 0 ) characterize the SWCC function that was adopted in the present model (recall Equations ( 4) and ( 5)).The values of these six parameters are obtained by regression-based curve fitting of the SWCC experimental data for a particular soil.
The following sections describe how values for the aforementioned 14 model parameters are determined from experimental data.

Elastic Parameters
The parameters K 0 and G 0 enter the formulation through the Gibbs free energy function given by Equation (9).They are used in the definition of the elastic moduli K, G, and J as given by Equations ( 14)- (19).A value for the parameter K 0 can be determined from an isotropic compression test in which p' increases but q = 0, thus rendering J = 0 from Equation (16).From this test, the bulk modulus K is equal to the slope of the line tangent to the initial portion of the mean effective stress (p') versus volumetric strain curve, from which K 0 can be determined by employing Equation (17).The parameter G 0 is best determined from a test in which q is increased from zero at constant p'.Since the test starts at q = 0, η is small initially, thus rendering small values of J from Equation (16).From this test, the value of G is determined from the initial slope of the curve of deviator stress versus distortional strain.In such a test, the value of G 0 is approximately equal to 3G.Values of K 0 and G 0 can also be determined from the results of resonant column or bender element tests [42].23)- (27).Values of e 0 and λ correspond to the slope and intercept, respectively, of the aforementioned straight line.Finally, the values of M c and M e correspond to the slopes of the CSLs that best fit the ultimate states of undrained stress paths, plotted in p'-q space, for axisymmetric triaxial compression and extension tests, respectively.

Parameters Controlling State-Dependency
The parameters n b and n d control the state dependency.At peak shear stress, the state of the soil is on the bounding surface.Here the stress ratio (η) is equal to M b .Using Equation (20), n b is computed as follows: At the point of phase transformation from compressive to dilatational response, the stress point lies on the dilatancy surface, where the dilatancy function given by Equation ( 32) is zero.
Using Equation ( 21), n d is computed as follows:

Dilatancy Parameter
The dilatancy function, d, is obtained by dividing the irreversible volumetric strain increment by the irreversible distortional strain increment.Following Li and Dafalias [38], it is assumed that the irreversible strain increments are approximately equal to the total strain increments.Therefore d = .
, which can be obtained from the results of drained axisymmetric triaxial compression tests [68].With d, M d , and h known, the value of A d is then computed from Equation (32), giving the following:

Hardening Parameters
Finally, suitable values of the parameters h 0 and c h , which enter the formulation through Equation ( 29), are determined by fitting model simulations to deviator stress versus distortional strain data from axisymmetric triaxial tests.

Assessment of Predictive Capabilities
The predictive capabilities of the model are next assessed by comparing its simulations with two sets of axisymmetric triaxial test results.Included in these comparisons are simulations obtained using three other bounding surface plasticity models that were previously developed for unsaturated soils.

Simulation of the Behavior of Clean Sand
Russell and Khalili [32] performed suction-controlled axisymmetric triaxial compression tests on Kurnel sand, which is clean with no fines.For this sand, the specific gravity of the soil particles (G s ) is 2.65, the mean grain size (D 50 ) is 0.31 mm, the uniformity coefficient (C u ) is 1.83, and the minimum and maximum void ratios are 0.60 and 0.92, respectively.Experimental data for the SWCC, along with the best-fit curve given by Equation ( 5), are presented in Figure 7. CSLs in the space of q-p' and e −  2.
increment by the irreversible distortional strain increment.Following Li and Dafalias [38 it is assumed that the irreversible strain increments are approximately equal to the tota strain increments.Therefore  = ≈ , which can be obtained from the results o drained axisymmetric triaxial compression tests [68].With d, M d , and h known, the valu of A d is then computed from Equation (32), giving the following:

Hardening Parameters
Finally, suitable values of the parameters h0 and ch, which enter the formulatio through Equation ( 29), are determined by fitting model simulations to deviator stress ver sus distortional strain data from axisymmetric triaxial tests.

Assessment of Predictive Capabilities
The predictive capabilities of the model are next assessed by comparing its simula tions with two sets of axisymmetric triaxial test results.Included in these comparisons ar simulations obtained using three other bounding surface plasticity models that were pre viously developed for unsaturated soils.

Simulation of the Behavior of Clean Sand
Russell and Khalili [32] performed suction-controlled axisymmetric triaxial compres sion tests on Kurnel sand, which is clean with no fines.For this sand, the specific gravit of the soil particles (Gs) is 2.65, the mean grain size (D50) is 0.31 mm, the uniformity coeffi cient (Cu) is 1.83, and the minimum and maximum void ratios are 0.60 and 0.92, respec tively.Experimental data for the SWCC, along with the best-fit curve given by Equatio (5), are presented in Figure 7. CSLs in the space of q-p' and  − are show in Figure 8a,b, respectively.The complete set of model parameter values used in simulat ing the behavior of Kurnel sand is presented in Table 2.      [41].The matric suction was maintained consta tests.Associated with the results shown in Figures 9-12 were suctions of 51, 10 kPa, respectively.In these tests, the radial net stress varied from approximatel proximately 100 kPa, and the initial void ratio varied between 0.658 and 0.780.
From Figures 9-12, it is evident that the present model accurately simula perimental deviator stress versus distortional strain and volumetric strain ver  the present model, it is less accurate at critical state.In addition, for the radial net stress of 50 kPa, the model predicted the peak deviatoric stress more accurately than the other two models.However, at the higher radial net stress of approximately 100 kPa, the Russell and Khalili [32] model provides a somewhat more accurate simulation.the present model, it is less accurate at critical state.In addition, for the radial net stress of 50 kPa, the model predicted the peak deviatoric stress more accurately than the other two models.However, at the higher radial net stress of approximately 100 kPa, the Russell and Khalili [32] model provides a somewhat more accurate simulation.distortional strain (experimental data from Russell and Khalili [32]).

Simulation of the Behavior of a Silty Sand
The present model is next used to simulate the axisymmetric triaxial behavior of silty sand.The experimental data for this soil have been reported by Patil et al. [29].The soil consists of 55% sand, 37% silt, and 8% non-plastic clay-size fraction.For this soil, Gs is 2.67.The maximum dry density and optimum water content of the soil are 1.87 g/cm 3 and 12.2%, respectively.Experimental data for the SWCC, along with the best-fit curve given by Equation (5), are presented in Figure 13.Critical state lines in q-p' and ( )  3. The triaxial compression tests were performed at four different suction values (50, 250, 500, and 750 kPa) and three different values of radial net stress (100, 200, and 300 kPa).In Figures 15-18, the simulations obtained using the present model are compared with experimental data and with the simulations reported by Patil et al. [28].Values of suction, radial net stress, and initial void ratio for each test are indicated on these figures.distortional strain (experimental data from Russell and Khalili [32]).
From Figures 9-12, it is evident that the present model accurately simulates the experimental deviator stress versus distortional strain and volumetric strain versus distortional strain results.Comparing these simulations to those obtained using the bounding surface models of [32,41], it is evident that they are in better agreement with the experimental results.In particular, simulations from the model of Morvan, Wong and Branque [41] give a "stiffer" stress-strain response.This model is also unable to reproduce the critical state data, in terms of either volumetric strain or shear strength, especially at the higher radial net stress of 100 kPa.The model of Russell and Khalili [32], on the other hand, simulates a "softer" behavior as compared to the experimental data.Compared to the present model, it is less accurate at critical state.In addition, for the radial net stress of 50 kPa, the model predicted the peak deviatoric stress more accurately than the other two models.However, at the higher radial net stress of approximately 100 kPa, the Russell and Khalili [32] model provides a somewhat more accurate simulation.

Simulation of the Behavior of a Silty Sand
The present model is next used to simulate the axisymmetric triaxial behavior of silty sand.The experimental data for this soil have been reported by Patil et al. [29].The soil consists of 55% sand, 37% silt, and 8% non-plastic clay-size fraction.For this soil, G s is 2.67.The maximum dry density and optimum water content of the soil are 1.87 g/cm 3 and 12.2%, respectively.Experimental data for the SWCC, along with the best-fit curve given by Equation (5), are presented in Figure 13.Critical state lines in q-p' and e − p ′ /p re f (1 + g(ξ)) ζ spaces are shown in Figure 14a,b, respectively.The complete set of model parameter values used in simulating the behavior of the silty sand is presented in Table 3.The data of Patil et al. [29] provided in  shows that the post-peak tening increases with the increasing amount of suction.In particular, the tests at the h est suction value of 750 kPa (Figure 18) show apparent anomalies in the experimenta viator stress and volumetric behavior.The post-peak softening for this test has led to s localization (i.e., the formation of a shear band).The displayed behavior preclude simulation of all features of the behavior using any existing model for unsaturated However, to be consistent, the results of the present model are presented for all the It can be observed that the present model is able to simulate the trend of the deviator s and volumetric behavior shown in Figures 15-17.Compared to the bounding su model of Patil et al. [28], the present model yields better agreement with the experim data.In particular, the model of Patil et al. [28] predicts a stiffer initial deviator st axial strain response, while the present model is smooth and close to the experim data.Also, at the critical state, the present model shows a better overall agreement the experimental data.The data of Patil et al. [29] provided in  shows that the post-peak s tening increases with the increasing amount of suction.In particular, the tests at the hi est suction value of 750 kPa (Figure 18) show apparent anomalies in the experimental viator stress and volumetric behavior.The post-peak softening for this test has led to str localization (i.e., the formation of a shear band).The displayed behavior precludes simulation of all features of the behavior using any existing model for unsaturated so However, to be consistent, the results of the present model are presented for all the da It can be observed that the present model is able to simulate the trend of the deviator str and volumetric behavior shown in Figures 15-17.Compared to the bounding surf model of Patil et al. [28], the present model yields better agreement with the experimen data.In particular, the model of Patil et al. [28] predicts a stiffer initial deviator stre axial strain response, while the present model is smooth and close to the experimen data.Also, at the critical state, the present model shows a better overall agreement w the experimental data.and q-p' (data from [29]).q-p' (data from [29]).

Conclusions
A state-dependent, critical state-compatible constitutive model that accounts for the effect of matric suction on the hydro-mechanical response of unsaturated granular soils was presented in this paper.The model is formulated in a rate-independent bounding surface plasticity framework and incorporates a newer generation definition of effective stress.The model is an enhanced extension of the hyperelastic model for fully saturated granular soils that was originally proposed by Lashkari, Karimi, Fakharian and Kaviani-Hamedani [43].Based on the observations of physical phenomena during laboratory testing of unsaturated granular soils, the stress-strain response of the unsaturated soil is assumed to be exclusively elastoplastic, with no purely elastic region.To satisfy the law of conservation of energy, the elastic components of the elastoplastic strains are calculated from a Gibbs free energy function in a hyperelastic formulation.Associated with the most general form of the present model are fourteen parameters whose values must be determined for a particular soil.Six additional parameters characterize the SWCC.The model was calibrated for two different sets of unsaturated experimental results, one for pure sand and one for silty sand.The predictive capabilities of the model were assessed by comparing simulations generated using the model with experimental data and with the simulations of the same data from three previously developed bounding surface constitutive models for unsaturated soils.The present model was shown to accurately simulate the behavior of unsaturated granular soils at various levels of initial matric suction, confining pressure, and void ratio.In addition, in many instances, the present model gave more accurate simulations than earlier bounding surface models for unsaturated soils.The data of Patil et al. [29] provided in Figures 15-17 shows that the post-peak softening increases with the increasing amount of suction.In particular, the tests at the highest suction value of 750 kPa (Figure 18) show apparent anomalies in the experimental deviator stress and volumetric behavior.The post-peak softening for this test has led to strain localization (i.e., the formation of a shear band).The displayed behavior precludes the simulation of all features of the behavior using any existing model for unsaturated soils.However, to be consistent, the results of the present model are presented for all the data.It can be observed that the present model is able to simulate the trend of the deviator stress and volumetric behavior shown in Figures 15-17.Compared to the bounding surface model of Patil et al. [28], the present model yields better agreement with the experimental data.In particular, the model of Patil et al. [28] predicts a stiffer initial deviator stress-axial strain response, while the present model is smooth and close to the experimental data.Also, at the critical state, the present model shows a better overall agreement with the experimental data.

Conclusions
A state-dependent, critical state-compatible constitutive model that accounts for the effect of matric suction on the hydro-mechanical response of unsaturated granular soils was presented in this paper.The model is formulated in a rate-independent bounding surface plasticity framework and incorporates a newer generation definition of effective stress.The model is an enhanced extension of the hyperelastic model for fully saturated granular soils that was originally proposed by Lashkari, Karimi, Fakharian and Kaviani-Hamedani [43].Based on the observations of physical phenomena during laboratory testing of unsaturated granular soils, the stress-strain response of the unsaturated soil is assumed to be exclusively elastoplastic, with no purely elastic region.To satisfy the law of conservation of energy, the elastic components of the elastoplastic strains are calculated from a Gibbs free energy function in a hyperelastic formulation.Associated with the most general form of the present model are fourteen parameters whose values must be determined for a particular soil.Six additional parameters characterize the SWCC.The model was calibrated for two different sets of unsaturated experimental results, one for pure sand and one for silty sand.The predictive capabilities of the model were assessed by comparing simulations generated using the model with experimental data and with the simulations of the same data from three previously developed bounding surface constitutive models for unsaturated soils.The present model was shown to accurately simulate the behavior of unsaturated granular soils at various levels of initial matric suction, confining pressure, and void ratio.In

Figure 3 .
Figure 3. (a) Closed loop in net stress space vs.(b) open loop in effective stress space (adopted from Sheng [13]).

Figure 4 .
Figure 4. Bounding, critical, and dilatancy surfaces and stress ratio  in p'-q space.

Figure 5 .
Figure 5. (a) Critical state data for an unsaturated soil at different suction values (s) and (b) normalization of critical state data using Equation (23) (data from Chiu and Ng [10]).

Figure 5 .
Figure 5. (a) Critical state data for an unsaturated soil at different suction values (s) and (b) normalization of critical state data using Equation (23) (data from Chiu and Ng [10]).
state effective confining pressures at saturated and unsaturated conditions, respectively.

Figure 6 .
Figure 6.Variation of the function f(s) with matric suction.
The seven parameters e 0 , λ, a, b, ζ, M c , and M e are associated with the definition of the critical state.Values of e 0 , λ, a, b, and ζ are determined by fitting a straight line to data plotted in e − p ′ p re f (1+g(ξ)) ζ space.The parameters a, b, and ζ are then obtained using Equations ( Figure 8a,b, respectively.The complete set of model parameter values used in simulating the behavior of Kurnel sand is presented in Table2.

Figure 7 .
Figure 7. Experimental and analytical soil-water characteristic curves for Kurnel sand.

Figure 7 .
Figure 7. Experimental and analytical soil-water characteristic curves for Kurnel sand.

Figure 8 .
Figure 8. Critical state line for unsaturated Kurnel sand in the space of (a)  −

Figures 9 -
Figures 9-12 compare the model simulations with experimental data and with the simulations obtained using the bounding surface models of Russell and Khalili [32] and Morvan, Wong and Branque [41].The matric suction was maintained constant in these tests.Associated with the results shown in Figures 9-12 were suctions of 51, 100, and 200 kPa, respectively.In these tests, the radial net stress varied from approximately 50 to approximately 100 kPa, and the initial void ratio varied between 0.658 and 0.780.

Figure 9 .
Figure 9. Simulations from multiple bounding surface models for unsaturated soils at constant matric suction of s = 51 kPa: (a) deviator stress vs. shear strain, and (b) volumetric strain vs. distortional strain (experimental data from Russell and Khalili [32]).

Figure 9 .
Figure 9. Simulations from multiple bounding surface models for unsaturated soils at constant matric suction of s = 51 kPa: (a) deviator stress vs. shear strain, and (b) volumetric strain vs. distortional strain (experimental data from Russell and Khalili [32]).

Figure 9 .
Figure 9. Simulations from multiple bounding surface models for unsaturated soils at constant matric suction of s = 51 kPa: (a) deviator stress vs. shear strain, and (b) volumetric strain vs. distortional strain (experimental data from Russell and Khalili [32]).
in Figure14a,b, respectively.The complete set of model parameter values used in simulating the behavior of the silty sand is presented in Table

Figure 14 .Figure 14 .
Figure 14.Critical state line for unsaturated silty sand in the space of (a) e − p ′ p re f (1+g(ξ))ζ and (b) q-p' (data from[29]).The triaxial compression tests were performed at four different suction values (50, 250, 500, and 750 kPa) and three different values of radial net stress (100, 200, and 300 kPa).In Figures15-18, the simulations obtained using the present model are compared with experimental data and with the simulations reported by Patil et al.[28].Values of suction, radial net stress, and initial void ratio for each test are indicated on these figures.

Table 1 .
Strain decomposition for unsaturated granular soils.

Table 1 .
Strain decomposition for unsaturated granular soils.

Table 2 .
Model parameters used to simulate the behavior of unsaturated Kurnel sand.

Table 2 .
Model parameters used to simulate the behavior of unsaturated Kurnel sand.Figures 9-12 compare the model simulations with experimental data an simulations obtained using the bounding surface models of Russell and Khal Morvan, Wong and Branque

Table 3 .
Model parameters used to simulate the behavior of unsaturated silty sand.

Table 3 .
Model parameters used to simulate the behavior of unsaturated silty sand.

Table 3 .
Model parameters used to simulate the behavior of unsaturated silty sand.