Numerical modeling of unsaturated soil behavior considering different constitutive models

Recent advances, not only in fluid flow but also in soil mechanics, have allowed the understanding and forecasting of common engineering problems such as slope stability, soil shrinkage and soil collapse. However, owing to limited access to data or more sophisticated numerical tools, the modeling of soil behavior is usually carried out considering simpler constitutive models which cannot predict some important features of unsaturated soils. This study is focused on the numerical modeling of unsaturated soils, adopting four constitutive models based on theories of elasticity and plasticity. For each model, a numerical simulation of a circular footing resting over a soil that is subject to drying and wetting processes is analyzed. Through the comparison of results, it is possible to highlight the use of more sophisticated constitutive models for unsaturated soil behavior, particularly forecasting the phenomenon of pore collapse during wetting processes.


Introduction
Several developments in soil mechanics have been progressively generalized in order to incorporate the effects of new phenomena and new variables to fit real soil behavior. Through Terzaghi's effective stress principle [1], it was recognized that not only is the soil itself problematic, but its complexity also arises from its interaction with water. The theory of poroelasticity [2] extended Terzaghi's principle and highlighted the rigorous coupling between mechanical and hydraulic variables. Since then, different non-linear models based on concepts such as yield or limit state surfaces have been adopted to represent the behavior of saturated soils. The soil mechanics theory was extended by the consideration of unsaturated soils in which the pore space is occupied by water and air. Bishop [3] was one of the first to establish an equivalent effective stress to that of Terzaghi for application in unsaturated soils. However, such stress, also known as Bishop stress, was not applied in further studies owing to the lack of a physical definition and experimental difficulties for its assessment. Later on, Fredlund and Morgenstern [4] proposed that, for a proper understanding of the shear strength of unsaturated soils, it was necessary to consider two independent stress variables: the suction and the net stress. Taking into account those variables, Alonso et al [5] introduced the Barcelona Basic Model (BBM) to represent both the strength and the deformability of unsaturated soils. That model was based on the critical state theory and developed through laboratory tests performed with suction control. Since then, several researchers have applied the BBM [6][7] [8] and extended its applications to more comprehensive problems [9] [10]. However, the use of the suction and net stress does not allow a natural recovery of the effective stress concept for saturated soils. Therefore, several studies have proposed the use of a "generalized effective stress" to model the soil behavior under saturated/unsaturated conditions [11][12] [13]. That concept has allowed the modeling of some important behavior of unsaturated soils such as their collapse under complex wetting processes [14][15] [16]. Furthermore, if this concept is adopted in other constitutive models, it can enable a fair comparison among them in terms of formulation and result analyses.
In this study, the formulation of some mechanical constitutive models is reviewed regarding their capabilities for a proper soil behavior modeling. For that purpose, a generalized effective stress that includes the suction effect on the soil mechanical behavior and allows the recovery of the effective stress concept in saturated conditions is adopted. Three classical constitutive models are considered: linear elastic (LE), non-linear elastic (NLE) and Modified Cam-Clay (MCC). In addition, a more comprehensive model that includes the generalized effective stress into the formulation of the BBM is used [15]. This model is referred in this paper as Barcelona Modified Model (BMM). To represent the fluid flow under unsaturated conditions, the van Genuchten model is applied. The presented numerical tests highlight the key differences between relatively simple and more comprehensive models to represent the behavior of soil foundations under drying and wetting processes.

Constitutive relations
In this section, the constitutive relations adopted to represent the mechanical and hydraulic behavior of unsaturated/saturated soils are introduced. The signs of the stress, strains and pressures follow the convention of soil mechanics, with positive stresses in compression.

Mechanical constitutive model
In this study, the concept of generalized effective stress [11][12] is given as follows in which is the total stress tensor; ‫‬ and ‫‬ ௪ are the air and water pressures, respectively; is the vector that introduces the influence of pore pressure on the effective stress; and ܵ ௪ is the effective degree of saturation given as: where ܵ ௪ is the degree of water saturation, ܵ ௪ is the residual value at very high suction values, and ܵ ௪ ௦ is the maximum value (normally equal to 1).
Observe that Eq. (1) defines only a constitutive variable that includes the combined effect of the net stress ( − ‫‬ ) and the suction (ܵ ௪ ሺ‫‬ − ‫‬ ௪ ሻ). This concept is different from the classical formulation of the BBM that considers both the net stress and the suction as constitutive variables. Once Eq. (1) is rewritten in an incremental form and any change in air pressure is disregarded, as it is usual in most geotechnical applications, the following equation is obtained For saturated soils, changes in the effective stress are related with strain changes (∆) through the stiffness matrix (). However, for unsaturated soils, an additional component related to the suction ‫ݏ(‬ = ‫−‬ ௪ ) must be considered to obtain changes in the generalized effective stress: where is a vector that includes the effect of suction on the development of plastic strains. Both and can be obtained through [11]: in which is the elastic stiffness matrix; ி = ‫߲/ܨ߲‬ ᇱ is the gradient of the yield function ‫;)ܨ(‬ = ߲ܲ/߲ ᇱ is the gradient of the plastic potential function (ܲ); ‫ܥ‬ = ‫ݏ߲/ܨ߲‬ is the derivative of the yield function with respect suction; and ‫ܪ‬ is the hardening/softening parameter.
Observe that for elastic models, = and = .
For the definition of the elastic stiffness matrix in the linear elastic model (LE), two constant parameters must be given: the bulk ‫ܭ(‬ ா ) and the shear ‫ܩ(‬ ா ) moduli. On the other hand, in the non-linear elastic model (NLE), the variable bulk ‫ܭ(‬ ோ ) and shear ‫ܩ(‬ ோ ) moduli are stress dependent and given as where N is the swell index, Q is the Poisson ratio, ‫′‬ is the mean generalized effective stress, and X is the specific volume.
Eqs. (7) and (8) where ‫ݍ‬ is the von Mises or deviatoric stress, ‫ܯ‬ is the slope of the critical state line in ‫ݍ-′‬ space, and ‫‬ is the preconsolidation pressure that controls isotropic hardening/softening. The preconsolidation pressure is related to the plastic volumetric strain (ߝ ௩ ) by: where O is the slope of the normal compression line for saturated soil, while the soil specific volume is updated according to (11) in which ܰ is known as the specific volume of normal compression line at unit pressure and it is a characteristic property of the type of soil.
In Eq. (9), the slope of the critical state line ‫ܯ‬ provides the shape of the failure surface in the deviatoric plane. Such parameter is expressed as a function of the Lode angle (θ) and of the friction angle of the soil at critical states (φ) as follows: Observe that suction is not included in the yield function of the MCC. Then, term ‫ܥ‬ in Eq. (6) is null. Therefore, the suction does not affect the elasto-plastic stiffness matrix in that model. For the Barcelona Modified model (BMM), the same elastic parameters are described in Eqs. (7) and (8). Furthermore, the following yield function is adopted where ‫‬ * is the preconsolidation pressure given as in which ‫‬ is a reference mean stress and O * is the slope of the normal compression line as follows where ‫ݎ‬ and ߚ represent model parameters. Observe that under unsaturated conditions, the yield function of the BMM depends on suction. Therefore, if the soil is over the yield surface, term ‫ܥ‬ in Eq. (6) is not null. Consequently, suction can trigger further plastic deformations. On the other hand, under saturated conditions with ‫ݏ‬ = 0 and O * ൌ O, the yield function of the BMM becomes that of the MCC model. Fig. 1 shows the yield function of the BMM considering the concept of generalized effective stress. On the top left, a projection over the plane ‫′‬ vs ‫ݍ‬ and the Critical State Line (CSL) with an inclination given by the parameter M are shown. Note that, owing to the use of the generalized effective stress, the yield surface is located between ‫‬ ᇱ ൌ 0 and ‫‬ ᇱ 0. Such feature is different from the BBM, in the case that the horizontal axis is given by the mean net stress that can be negative. On the top right of Fig. 1, the projection of the yield surface over the plane ‫′‬ vs ‫ݏ‬ through the Loading-Collapse (LC) and the curve that define the occurrence of irreversible strains after suction increase (SI) are shown. Observe that, the curve SI is given by a straight line defined by the maximum past suction ever experienced by the soil ‫ݏ(‬ ሻ. At the bottom of Fig. 1, it is seen a tridimensional view of the yield surface. It can be noticed that the size of the yield surface increases with suction while for the saturated condition ‫ݏ(‬ ൌ 0), the yield surface has the size of the MCC model. Finally, to complete the formulation of both MCC and BMM, it is necessary to define plastic potential functions (ܲ). For the sake of simplicity, this study assumes an associated flow rule in both models.

Fluid-flow constitutive model
A water retention curve defined in terms of the water content and suction is usually adopted to model fluid-flow of unsaturated soils. However, for hydromechanical modeling, it is more convenient to express such curve in terms of the degree of saturation. Different water retention relations can be found in the literature [17][18] [19] [20]. In this study, we adopt the retention curve proposed by van Genuchten [18] given by where ܽ, ݊ and ݉ represent van Genuchten parameters. Those parameters can be determined through field or laboratory testing. Observe that for saturated conditions, Eq. (16) becomes ܵ ௪ ൌ ܵ ௪ ௦ . An example of a retention curve is shown in Fig. 2. Such curve was built considering the following parameters: ܽ=1.0m -1 , ݊=0.5, ݉=1.0, ܵ ௪ = 0.2 and ܵ ௪ ௦ =1.0. It can be observed that for those parameters, suction values higher than 1000 kPa would be necessary to reach the minimum degree of saturation. The other characteristic relation adopted for a proper representation of fluid-flow of unsaturated soils is the permeability function shown in Fig. 3. Such function is used to obtain the relative permeability of unsaturated soils (݇ ௪ ) from known suctions. According to the van Genuchten model, the permeability function is given as Finally, the permeability matrix of the soil is obtained as ൌ ݇ ௪ , where is the permeability matrix for saturated conditions. Observe that for saturated conditions, ܵ ௪ ൌ 1 and ݇ ௪ ൌ 1. Therefore, ൌ .

Governing equations 3.1 Stress equilibrium equation
Under quasi-static conditions, the stress equilibrium equation can be written as where is the gravity vector and ߩ is the bulk density of the soil that under unsaturated conditions is given by in which I = ሺX − 1ሻ/X is the soil porosity, ߩ ௦ and ߩ ௪ are the density of solid grains and water, respectively. Following the definition of generalized effective stress and using Eqs. (3) and (4), Eq. (18) can be rewritten in a rate form as where ሶ , ‫‬ሶ ௪ and ߩሶ represent the rate of the strain tensor, water pore pressure and bulk density, respectively.

Mass conservation
The mass balance of water-flow through an unsaturated soil is given by where ߚ ௪ is the water compressibility and ௪ is the macroscopic Darcy velocity given by in which ߤ ௪ is the water viscosity.

Solution
The stress equilibrium and mass conservation equations are solved through GEOFLUX3D [15], an in-house finite element simulator for the analysis of problems with fully implicit hydromechanical coupling. GEOFLUX3D has been successfully employed in the modeling of geotechnical applications such as earth dams [21], groundwater flow [22], slope stability [23] and soil consolidation problems [24].

Numerical tests
In the numerical tests, we consider the problem of a flexible circular footing of radius R=0.5m, resting on a homogeneous soil layer. Such a problem is a classic test for the verification of an unsaturated soil model, owing to the singularity at the axisymmetric axis and the strong rotation of the principal stresses during loading. Similar problems have been analyzed in previous studies [14][10] [25]. Fig. 4 shows the geometry and the finite element mesh used in the numerical tests. The quadrilateral elements adopt quadratic and linear interpolation for displacements and pore pressures, respectively. In total, the mesh has 576 elements with 4-point Gauss integration and 1825 nodes. The material properties are shown in Table 1. The mechanical behavior of the soil layer is represented through four constitutive models: LE, NLE, MCC and BMM. Therefore, four analyses are performed. The fluidflow is represented through the van Genuchten model considering the retention and permeability curves given in Fig. 2 and Fig. 3, respectively. For the sake of simplicity, any hysteresis of soil-water characteristic curves is disregarded. Initially, it is assumed that the soil is saturated and the pressure is hydrostatic with the water level position on the surface. In order to facilitate the comparison among the models, the whole soil layer is under an isotropic effective stress of 10 kPa. Observe that for such condition, the initial bulk and shear moduli in the soil layer are 1000 and 461.54 kPa, respectively. For application of the MCC and BMM, an initial preconsolidation pressure ‫(‬ ) of 58.8 kPa and a maximum past suction ‫ݏ(‬ ) of 100 kPa are adopted. The soil presents an isotropic intrinsic permeability of 10 -13 m 2 .
Regarding the boundary conditions, the horizontal displacements are restricted at the model lateral boundaries while both horizontal and vertical displacements are restricted at the model bottom. In addition, the bottom and the lateral boundaries of the model are closed for flow. During the analyses, three phases are considered. In the first phase, the soil is dried through the application of a uniform suction at the top of the model. Such suction increases linearly along the time up to reach 40 kPa in 40 days. Then, in the second phase, the circular footing is incrementally loaded to 100 kPa during the next 40 days. As in this phase, the footing is loaded by a uniform vertical pressure; therefore, it is perfectly flexible. Finally, in the third phase, the top of the model, between points B and C, is wetted to zero suction in further 40 days. For result comparison, the vertical displacements at three observation points A, B and C (shown in Fig. 4) are considered. Fig. 5 shows the settlement history at point A along the whole simulation. At the end of the drying phase, it is noticed a maximum settlement of 7.3 cm obtained with the LE model while in the other models, the settlement is the same and reaches values of 4.5 cm. The similar behavior in the analyses with NLE, MCC and BMM are explained for the occurrence of pure elastic deformation during this drying phase. Observe that the applied suction during the drying phase was lower than ‫ݏ‬ . In comparison with the LE model, the results are consistent since the LE model considers constant elastic parameters during the whole simulation. On the other hand, the other models consider dependence on the evolution of the mean effective stress which are increased as the soil is dried. Consequently, the stiffness of the soil is increased, triggering lower settlements than those predicted with the LE model.
During the loading phase, it can be noticed that the settlement continues to raise with more intensity in the LE model, reaching at the end a maximum value of 16.6 cm. In turn, the other three models follow the same path up to 53 days and then, the MCC model triggers more intense settlements, reaching at the end of the loading phase, the highest settlement (18 cm). Such behavior in the MCC model is explained by the occurrence of plastic deformation for loading higher than 37 kPa. On the other hand, the settlement with the NLE model and the BMM is the same until 67 days when appear plastic deformation in the BMM. These late deformations, in comparison with those seen with the MCC model, are explained by the effect of the suction in the yield envelope of the BMM. As the suction was increased during the drying phase, the yield surface of the BMM grew up. Consequently, plastic deformations appear later for loading higher than 69 kPa. At the end of the loading phase, the BMM and the NLE analyses triggered 9.4 and 7.6 cm, respectively.
Finally, during the wetting phase, we can observe that the settlement in the analyses with the LE, NLE and MCC models is reduced owing to increments in the soil volume. However, in the BMM, the settlement continues raising. It can be explained by the stress state that is located over the yield surface at the end of the loading phase. Under such condition, wetting causes plastic yielding and the soil undergoes a decrease in volume. That behavior is in practice responsible for soil collapse and can be only captured by the BMM model.   6 shows the settlement history at point B. A similar behavior to those observed at point A is noticed, but lower settlements are triggered during the loading and wetting phases. Moreover, during the loading phase, the path of the BMM is practically the same as the NLE. However, at the end of that phase, the higher settlements predicted with the BMM indicate the occurrence of plastic deformations. Then, during the wetting phase, the settlements are increased with the BMM, showing a similar behavior of that observed at point A.   remain with elastic behavior. Consequently, during the loading phase, the settlement in this region increases but with lower intensity than that observed at points A and B. Then, during the wetting phase, the settlement at point C is reduced owing to soil swelling. In order to compare regions with plastic deformation, it is adopted an equivalent plastic strain (ߝ ) given as where ߝ ଵ , ߝ ଵ and ߝ ଵ are the principal plastic strains. Fig. 8 shows the equivalent plastic strains obtained with the MCC model at the end of the loading phase (on the left) and at the end of the wetting phase (on the right). It is observed that the plastic deformations are concentrated around the footing. Furthermore, from the comparison of both contour plots, we can notice that the plastic region at the end of the loading phase is the same as that obtained at the end of the wetting phase.    9 shows similar plots as those of Fig. 8, but considering the results obtained with the BMM. In comparison with results of the MCC model, we can notice a narrow plastic region at the end of the loading phase. However, at the end of the wetting phase, we observe a significant increase of the plastic region owing to the collapse of the soil foundation.

Conclusions
In this study, a comparison of different constitutive models to represent the behavior of unsaturated soils is presented. Considering the concept of generalized effective stress in the constitutive modeling, an equivalence between the formulation of four constitutive models, from the simplest LE model to the most comprehensive BMM, is introduced. In addition to those models, the NLE and the MCC models are considered for numerical simulations of a circular footing resting over a shallow foundation submitted to drying, loading and wetting. Through the comparison of the predicted settlements with each model, different soil responses are observed. The effects of suction are particularly evident during the loading and wetting phases of unsaturated soils. Those results highlight the importance of an appropriate choice of constitutive model to analyze the behavior of unsaturated soils, particularly in the forecast of potential geotechnical problems such as soil collapse.