A Constitutive Model for Unsaturated soils based on a Compressibility Framework dependent on Suction and Degree of Saturation

The Modified Cam Clay model is extended to account for the behaviour of unsaturated soils using Bishop’s stress. To describe the Loading – Collapse behaviour, the model incorporates a compressibility framework with suction and degree of saturation dependent compression lines. For simplicity, the present paper describes the model in the triaxial stress space with characteristic simulations of constant suction compression and triaxial tests, as well as wetting tests. The model reproduces an evolving post yield compressibility under constant suction compression, and thus, can adequately describe a maximum of collapse.


Introduction
This paper presents a simple constitutive model for the numerical analyses of unsaturated soils.The proposed model is an extension of the Modified Cam Clay (MCC), in the unsaturated regime.It incorporates a Loading -Collapse (LC) Surface, derived on the bases of a compressibility framework proposed by the authors [1].
The utilized framework includes compression lines that are both dependent on suction (s) and on the effective (macro-structural) degree of saturation ( e r S ).This double dependence allows for a comprehensive description of an evolving compressibility during constant suction tests.

Compressibility framework
Sitarenios and Kavvadas [1] extended an intrinsic compressibility framework for saturated soils to account for the effects of partial saturation.The framework utilizes Bishop's average skeleton stress, assuming e r S   , while in general it is in favour of an increased yield stress and a decreased post yield compressibility with increasing suction (s).
It comprises an enhancement of the Barcelona Basic Model (BBM) to account for the slope of the unsaturated compression lines in the Bishop's stress domain, based on a relatively simple idea to scale the effect of suction by incorporating the effective degree of saturation: The term (1-S r e ) γ in fact implies that an increase in degree of saturation acts as a structure degradation mechanism which progressively increases the post yield compressibility of the unsaturated material.Μaterial constant γ is added (it was not included in [1]), to scale the effect of degree of saturation, allowing for increased versatility.For the   s  , the well-known BBM equation is used: and after some algebra we end up with: In equation 3, λ is the slope of the compression lines under saturated conditions and r, β, γ are material constants.Parameters r, β originate from the BBM, nevertheless they do not fully retain their original role.The proposed framework (see fig. 1) suggests that unsaturated soil states, lay on a compression curve which originates from the saturated compression line on a reference pressure c P , while its slope is given by eq. 3. Thus they plot to the right of the saturated compression line in the metastable domain, a common assumption, reflecting the beneficial effect of inter-particle forces (induced by partial saturation) on soil stiffness.It additionally accounts for the increased yield stress that unsaturated soils exhibit when compressed.The latter, in the present framework can be calculated using a similar to BBM concept to obtain: * ( , ) 0 0 ( , ) Equation 4 gives the evolution of the apparent preconsolidation pressure with partial saturation for a material with a saturated preconsolidation pressure equal to P 0 *.
The similarity of eq. 4 with its BBM counterpart is apparent, nevertheless in the proposed framework its additional dependence on the effective degree of saturation, ensures that saturated states will always lay on the compression line of the saturated material, irrespectively of the level of the applied suction.This is a prerequisite for frameworks that utilize Bishop's stress as the latter coincides with Terzaghi's effective stress upon saturation.Frameworks solely dependent on suction cannot capture such a behaviour unless some kind of discontinuity, associated with the air entry value of the material, is assumed.Moreover, the dependence on degree of saturation, can reproduce a constantly evolving compressibility for soils compressed under constant suction, as a natural outcome of the soil's increased capacity for water retention with decreasing void ratio.Such a behaviour is schematically depicted by the grey dotted line in figure 1.
In the last few years, many researches proposed that the slope (i.e., [2]) and/or the position (i.e., [3]) of the unsaturated compression lines in the Bishop's stress domain can be solely dependent on degree of saturation with suction only implicitly influencing the behaviour through the Water Retention Curve and Bishop stress.Such frameworks address the issue of unique compression lines for saturated states, and can also represent a constantly evolving compressibility under a constant given suction level.
Nevertheless, we believe that a double dependence is more representative of the fundamental mechanisms describing the shift of the compression lines in the metastable domain.As was mentioned before, this "structure" inducing mechanism is attributed to the bonding effect of the forming inter-particle forces.As explained in [1], similar to Gallipoli's proposal [4] we believe that the intensity of this bonding effect depends both on suction and on the quantity of water in the soil's pores.
The latter can be better understood in a relatively coarse-grained soil, where capillary phenomena prevail.In such a case, inter-particle forces are mainly attributed to the forming water menisci.The strength of the forming menisci depends mainly on suction, while how many menisci are present in a given soil mass depends on degree of saturation.As a soil progressively saturates under a constant suction level, some of the menisci are destroyed and thus the intensity of the inter-particle forces drops.
On the other hand, examining the theoretical case, where suction changes under a given degree of saturation, the strength of "constant" menisci increases or decreases, following suction changes.This example, although severely simplifying and fictitious to some extent, is illustrative of the assumed behaviour.

Constitutive equations
The aforementioned framework is incorporated in the MCC model, to describe the LC curve.The model uses Bishop's average skeleton stress as its First Constitutive Variable (FCV), while also suction (s) and effective degree of saturation ( e r S ) are incorporated as extra constitutive variables, the three of them comprising the ensemble of the model's external variables.No further modifications are implemented, in fact subscribing to the idea that the shear strength behaviour as well as the behaviour in the elastic domain are fully interpretable in terms of Bishop's average skeleton stress, when the macro-structure's degree of saturation is used.
At this point we shall mention that although degree of saturation is assumed an external variable, in fact it is not an independent variable, as it is an outcome of the soil's Water Retention Curve (WRC).
Regarding the set of internal variables, alike MCC, the proposed model utilizes the specific volume (v), needed for the porous elastic behaviour and for the hardening rule, while the saturated preconsolidation pressure * 0 P is used as the hardening variable.
In the following lines, the model formulation is presented in the triaxial effective stress space ( , e r ps p S   q ) with the primes dropped for simplicity, as net stress ( p ) does not appear in the formulation, while dot over a symbol denotes an infinitesimal increment of the corresponding quantity.

Yield function
The yield surface is described by the following yield function: , , , , ( , ) 0 ee rr f p q a S P q M p P s S p     (5)   in combination with equations ( 4) and (3).M is the slope of the Critical State Line (CSL) in the p-q.The behaviour inside the Yield Surface is assumed elastic, while plastic states always lay on the Yield Surface. Figure 2 presents the yield surface in the p-q-s space.

Elasticity
The standard porous elasticity is utilized to describe the behaviour inside the elastic domain.The corresponding elastic increment of Bishop's mean stress p and of the deviatoric stress q are calculated through: where K, G are the Bulk and Shear Modulus, respectively.K depends on void ratio and on the level of the applied stress through: where  is the slope of the swelling lines.
A pressure dependant Shear Modulus is also assumed depending on; a) the value of the Bulk Modulus and b) the Poisson's ratio value ( ): Figure 2. The model's yield surface in p-q-s space.

Flow rule
The flow rule describes the onset of plastic strains during plastic loading.The MCC assumes an associated flow rule, where the plastic potential function (g) coincides with the yield function ( gf  ) and thus, plastic strain increments are described by the following expressions:  is the so called plastic multiplier and is defined as: where H is the Plastic Modulus, derived on existing assumptions, through the consistency condition.
Having defined the evolution of plastic strain increments, by further utilizing the basic kinematic decomposition of total strains to an elastic and a plastic part ( ep    ), we may write that: 3 ( ) 3 ( ) In eq.12 we may substitute p and q with eq.13 and eq.14 respectively, to finally compute  in terms of the total strain increments, suction and effective degree of saturation increments.

Hardening rule
As already mentioned, in the proposed model, the saturated preconsolidation pressure is selected as the hardening variable.The MCC volumetric, hardening rule is used to describe its evolution and thus, we may write: where  the MCC compressibility parameter describing the slope of virgin compression lines under saturated conditions.

Consistency condition -Plastic modulus
Finally, the consistency condition is employed to calculate the Plastic Modulus H. Consistency condition ensures that plastic states remain on the Yield Surface and is mathematically described by: Using eq.12 in eq.15 we can write that: While further substituting * 0 P with eq.16, Plastic Modulus can be finally expressed as: The derivatives of the yield surface with respect to the external variables and the hardening variable, necessary in the previous calculations, are given in Appendix I.

Model predictions
This section evaluates the model's predictions through a series of simple material point simulations, based on a fictitious set of parameters.
To do so, the model is implemented in a single point algorithm, while an explicit integration scheme is used to solve the constitutive equations.The algorithm is capable of imposing common laboratory stress paths.
Regarding simulations under unsaturated conditions, suction is an input to the algorithm while degree of saturation is calculated through the Gallipoli's WRM [5].
The selected WRM can reproduce an evolving degree of saturation even during compression under constant suction due to its void ratio dependence.Equation 21 is used to calculate the increment of degree of saturation with suction and volumetric strains.
Finally, to scale degree of saturation for macrostructure the Alonso et.al. [6] power law is used.The effective degree of saturation and its corresponding increment, are calculated according to: Table 1 presents the material parameters selected for the analysis.The selected parameters, although fictitious, are coherent, and representative of the typical range of values for a silt.Four (4) extra parameters with respect to the MCC model are needed for the mechanical behaviour.Of course additional variables are required for the WRM and to scale degree of saturation for macro-structure.
Table 1.Material parameters used in the analysis.We may observe that the proposed framework can successfully describe a non-linear increase in the yield stress with increasing suction, as well as a constantly increasing post yield compressibility for compression under constant suction.

Constant suction isotropic compression
In the elastic domain, as expected, the results coincide when plotted in terms of Bishop stress.Nevertheless, in terms of net stress, there is a clear decrease in the specific volume under constant net stress, associated with the drying stage, while during the subsequent compression under constant suction different swelling lines apply to DOI: 10.1051/ , 9 E 2016 different suction levels, while their slope evolves with net stress, as has been explained in Jommi [7].
Regarding the post yield compressibility, we may observe that the higher the suction the higher the specific volume that the soil can retain under any given net stress level, reflecting the beneficial effects of partial saturation in stiffening the response.Moreover, the constantly evolving compressibility reproduces convex net stress compression lines, indicating a decreasing initial post yield compressibility with suction, which constantly increases with compression, while at some point becomes even greater than that of the same material under saturated conditions.The reproduced compression lines are very reminiscent of frameworks usually underlying constitutive models which try to reproduce a maximum of collapse, while in the proposed model it is a natural outcome of the dependence on both suction and degree of saturation.

Constant suction triaxial compression
The isotropic compression simulations presented in the previous section, are repeated, and for each suction level imposed.Compression is now stopped at three different mean net stress levels ( p  500, 1000 and 1500kPa), and the soil elements are subjected to triaxial compression under drained, constant suction conditions until critical state is reached.
All elements are sheared from normally consolidated condition (wet side), as the selected net stress levels are higher than the apparent preconsolidation pressure developed during drying.All elements exhibit a strainhardening response and a contractant volumetric behaviour, while ultimate conditions dependent on both suction and net stress.
Figure 5 presents the ultimate stress points reached in terms of net stress.We may observe that as expected, different CSLs apply for different suction levels in the net stress domain.The tensorial translation s P corresponding to each line has been included in figure 2, to compare with the e r sS  .Figure 6 plots the normalized ratio ( , , 0 / u s u s qq  ) along with M(s), representing the strength increase with suction, with respect to the shear strength of the same material under the same confining stress when saturated under zero suction.
As expected, in general it follows the e r sS  term, nevertheless due the fact that the points corresponding to the same suction level, do not correspond to a common e r S value the obtained s P values represent a weighted behaviour.
Another effect of the evolving r S is that the CSLs are not parallel to the saturated CSL but their slope depends on suction.Figure 6 plots the evolution of M with suction.We may observe that it slightly increases with suction until s=1000kPa and then practically stabilizes.We shall mention that various researches have experimentally defined such suction dependent CSLs with Wheeler and Sivakumar [8] being amongst the first.The model can represent a non-linear increase in the shear strength with suction, similar to what is usually observed in experimental results.A maximum is reached for s=1MPa (for the selected parameters), while shear strength drops for higher suction values.Moreover, the results indicate that as net stress increases, the effect of partial saturation in strengthening the soil behaviour decreases.
Qualitatively, the evolution of shear strength follows the evolution of the e r sS  , a common outcome of Bishop's stress incorporation in effective stress models, and also very helpful in terms of calibration.Scaling the degree of saturation for macrostructure though is very essential, for the e r sS  term not to obtain extreme values.

Wetting tests -volumetric collapse
Finally, once again based on the isotropic compression tests previously reported, wetting tests are simulated under different suction and net stress levels.Wetting tests are simulated by gradually reducing suction to zero, while maintaining the net mean stress constant.
Figure 7 presents indicative results for simulations regarding s=100kPa under seven different constant confining pressures given in the graph.We may observe that swelling behaviour is reproduced inside the elastic region while for stress levels higher than the preconsolidation pressure, the stress state lays on the LC curve and thus, collapse strains are reproduced.In more details, swelling or collapse behaviour is reproduced depending on the level of the applied net stress.Moreover, the calculated strains increase with increasing suction under any given net stress level, while when examined under a constant suction level, collapse strains increase with increasing net stress, until a maximum is reached, followed by a reduction in the amount of collapse.
For the parameters selected, a maximum of collapse appears at 1 p  MPa, for the elements with an initial suction up to 250 s  kPa, while for higher initial suction no maximum appears up to the net stress tested.Compression to higher net stress levels is required with increasing suction, due to the increased soil stiffness.

Conclusions
The mathematical formulation and typical predictions of a constitutive model for unsaturated soils were presented.The model is based on the MCC with the addition of a relatively simple framework to account for structure induced by partial saturation.
The simulation results indicate that at least in a qualitative way the model can represent typical behaviour of unsaturated soils, like increasing yield stress and decreasing post yield stiffness with increasing suction, as well as non-linear evolution of the shear strength with suction.Moreover, as far as a void ratio dependent WRM is used, it can also reproduce a maximum of collapse.
Further evaluation of the proposed framework against actual experimental data is still necessary.

Figure 3
Figure 3 plots the WRCs assumed for different specific volume values as well as the corresponding evolution of the e r sS  term.

Figure 4
Figure 4 presents the results of isotropic, constant suction compression tests under seven (7) different suction levels, ranging from 50 s  kPa to 5 s  MPa, thus covering an extended portion of the WRC.Initially over-consolidated soil samples with * 0 100 P  kPa are dried under constant net mean stress 50 p  kPa by raising suction to the desired level.To calculate the initial specific volume, 2.10 iso N  was assumed.Isotopic compression up to a men net stress of 5 p  MPa is finally imposed.

Figure 4 .
Figure 4. Constant suction isotropic compression tests under various suction levels.Behaviour on the ln , v p p  planes.

Figure 5 .
Figure 5.The Critical State Lines on the p net -q.

Figure 6 .
Figure 6.Evolution of shear strength with suction.

Figure 7 .
Figure 7. Wetting tests at various stress levels for 100 s kPa  .

Figure 8
Figure 8 summarizes the obtained volumetric response.Only the volumetric strains associated with the wetting phase are plotted.The results reveal that, the proposed model can reproduce typical collapse behaviour, usually observed in experimental results (i.e.,[9]).