A viscoplastic constitutive model for unsaturated soils based on nonstationary flow surface theory

A viscoplastic constitutive model is developed to describe the viscoplastic behavior of unsaturated soils. The proposed model accounts for the stain‐rate and suction effects on yield by adopting an unsaturated isotach concept. The nonstationary flow surface theory (NSFS) is applied for modeling the viscoplastic behavior, with the yield surface which can evolve with the viscoplastic strain, viscoplastic strain rate and suction. Meanwhile, the progressively hardening concept is adopted for reproducing the viscoplastic behavior of soil at an over‐consolidated state. For the validation, a series of loading conditions are considered based on the data from the literature. Results show that the proposed model is able to reproduce the main viscoplastic behaviors of unsaturated highly compacted Gaomiaozi (GMZ) bentonite, Glenroy silt and Qianjiangping landslide (QL) soil, including CRS compression tests, rate‐dependent triaxial shear tests, and triaxial creep tests.


INTRODUCTION
Over the past decades, the time-dependent behavior of geomaterials has been widely investigated in geotechnical engineering. 1,2The time-dependent deformation of geomaterials is mainly induced by its viscosity, which has an adverse impact on the long-term safety and performance of subgrade, foundation, tunnel, and so forth.As a result, a large number of experimental and modeling works were conducted in this field.
4][5][6] Besides, constant strain tests (CRS) were performed on reconstituted and natural clays in oedometer, allowing a unique relationship among stress, strain and strain rate to be identified. 7,8The yield stress and shear strength were found to increase with increasing strain rate.Similarly, CRS and triaxial tests under different suctions and strain rates were carried out to understand the effect of unsaturated parameters (degree of saturation, suction, etc.) on the time-dependent behavior of unsaturated geomaterials. 6,9ased on the experimental results for saturated soils, numerous empirical relations and rheological models were proposed to predict the time-dependent behavior of soils, including creep, stress relaxation and rate effects. 2,10,11Some more advanced relationships between stress, strain and time (or strain rate) were constructed within the elasto-viscoplastic (EVP) frame, for instance, the over-stress theory 12 and the non-stationary flow surface (NSFS) theory. 11Based on the over-stress theory, many viscoplastic models using different over-stress functions were built to describe the strain rate effect, creep and stress relaxation behavior of saturated soils.In particular, 3D EVP models were derived from the generalized viscoplastic overstress theory (e.g., Qu et al. 13 ) Besides, using the overstress concept, Yin et al. 14 explored the anisotropic viscous behavior of natural soft soils.It is however worth noting that the overstress-based models cannot meet the consistency condition.In other models, the isotach approach proposed by Suklje 15 was generally adopted to describe the viscosity of soils, which was usually incorporated into the flow yield surface theory. 11,16The NSFS-based models can meet the requirement of consistency condition, and a higher convergence rate can be achieved in comparison with the overstress-based models. 17Several constitutive models have successfully applied the NSFS theory in describing the viscoplastic behavior of soils. 11,18,19However, most of these EVP models were developed for saturated soils.
To days, the few existing EVP models for unsaturated soils are based on either the overstress theory 20 or the ratedependent theory. 21,22The constitutive model built by Wu 20 is able to describe the time-dependent behavior of unsaturated Glenroy silt.However, the model still needs to be verified with more experimental data and with consideration of more stress paths.Mac et al. 22 employed the bounding surface plasticity theory to define the yield function by using the viscoplastic strain as one of the variables, while neglecting the influence of strain rate on the direction of viscoplastic strain.Based on the isotach concept, De Gennaro and Pereira 21 extended the Barcelona Basic Model (BBM) to the effect of loading rates for chalk at different suctions.Nevertheless, this model cannot describe the stress relaxation phenomenon and viscoplastic effect in low stress range due to the adoption of total volumetric strain rate as the viscoplastic variable and single yield function.As aforementioned, these unsaturated EVP models cannot reproduce some viscoplastic features (e.g., stress relaxation and creep in low stress range) of unsaturated soils due to the disadvantages in theory.In this case, it seems that the NSFS theory is regarded as a promising approach to model the viscoplastic behaviors of unsaturated soils.However, to the authors' knowledge, there has been no development of NSFS-based unsaturated viscoplastic model.
In this study, the NFSF theory is first extended to the unsaturated condition.Then, a viscoplastic model for unsaturated soils is developed based on the unsaturated NSFS theory.Within the consistency viscoplastic framework, the proposed viscoplastic model enables the description of the viscoplastic behavior of unsaturated soils.Besides, a gradual hardening law is incorporated into this model, allowing the description of smooth stress-strain curves and creep phenomenon at a very low stress.[25]

2
NONSTATIONARY FLOW SURFACE THEORY FOR UNSATURATED SOILS

Development of isotach approach for unsaturated soils
In order to model the viscoplastic behavior of saturated soils, the isotach approach proposed by Suklje 15 has been widely used to describe the relationship of stress-strain-strain rate.In this approach, it is assumed that there is a unique relationship between the effective stress, the strain and the strain rate. 7,8,15This relationship is characterized by the effect of strain rate on the yield stress and can be expressed by Equation (1) 7 : where  ′  is the yield stress and ε is the volumetric strain rate.Following this approach, soils may experience both compression under constant strain rate or delayed compression under constant stress.It is worth noting that each isotach line has an equivalent rate-dependent yield stress which is larger for a higher strain rate, as shown in Figure 1.By mathematically normalizing this relationship, it can be applied in constitutive modeling, as presented by Equation (2). 8The parameter  A represents the slope of the linear relationship of yield stress and strain rate in the log-log plane, indicating the viscosity between soil particle and thus affecting the creep behavior of soils.However, with Equation (2), stress relaxation cannot be reproduced because the total volumetric strain keeps constant, which means that strain rate is equal to zero, defining a constant stress state.To describe the stress relaxation phenomenon, Laloui et al. 26 proposed to replace the relationship of stress-strain-strain rate by the relationship of F I G U R E 1 Schematic illustration of stress-strain-strain-rate behavior.

(A) (B)
F I G U R E 2 Variation of yield stress with strain rate: (A) Lixhe chalk (data from Priol et al. 27 ); (B) GMZ bentonite (data from Qin et al. 9 ) .
stress-strain-viscoplastic strain rate, as shown in Equation (3).Qiao et al. 11 modified the original formulation by adopting Equation (4) to model the unique relationship for saturated soils and also pointed out that the determination of the related parameters is the same while neglecting the elastic strain rate effect.In Equation ( 4), the parameter  is replaced by the effective yield stress  ′ , ε , at a given reference state, enabling the strain rate dependency of yield stress to be considered.is the reference effective yield stress at the reference viscoplastic volumetric strain of ε , .For unsaturated soils, the relationship between the net yield stress and the strain rate is found to be similar to that of saturated soils-as the applied strain rate increases, the corresponding net yield stress also increases, as shown by Priol et al. 27 on Lixhe chalk and by Qin et al. 9 on GMZ bentonite.From Figure 2 it appears that the yield stress at each suction depends on the applied strain rate, showing a linear yield stress-strain rate relationship.Besides, the slopes of the linear relationship seem to depend on suction and the type of fluid.This observation suggests that the strain rate dependency is also suction-dependent.Thereby, the strain rate dependency of net yield stress can be expressed by replacing Effect of suction on  A (): (A) rockfill (data from Oldecop & Alonso 28 ); (B) GMZ bentonite (data fromYe et al. 29 & Qin et al. 9 ).
in Equation ( 4) with   ().By adopting a reference stress and a reference viscoplastic volumetric strain rate, the yield stresses at different suctions can be described with a unified equation (Equation 5).Hence, the unsaturated effect is also incorporated into the isotach approach.
( ε  , where  A () is the slope of the linear relationship and reflects the viscosity of unsaturated soils at a given suction , σ  () and ε , are the reference net yield stress and the corresponding reference viscoplastic volumetric strain rate at the same suction, respectively.
To describe the changes of  A () with suction, the slopes in Figure 2 are determined and presented in Figure 3.It appears that the slope  A () decreases with increasing suction, which reflects the effect of suction on the secondary consolidation.Besides, the collected experimental data indicate that there is a liner relationship between  A () and suction.Equation (6)  recommended by Qin et al. 9 is adopted to express the suction dependence of  A .Here, by adopting  +   as a variable, the infinity at saturated state ( = 0) is avoided.
where  and  are model parameters governing the evolution of   () with ,   is the atmospheric pressure taken equal to 0.1 MPa.
According to the unsaturated isotach approach, the coupled effects of strain rate and suction on the net yield stress can be schematically illustrated, as shown in Figure 4.At a given suction  1 , the compression line at the volumetric strain rate ε,1 is parallel to that at ε,2 .The initial net yield stress  ,1 at ε,1 is smaller than the  ,2 at larger volumetric strain rate ε,2 and also evolves with the accumulation of viscoplastic volumetric strain after reaching the normal consolidated state, showing the unique stress-strain-strain-rate relationship.When suction increases to  2 at a given strain rate, the corresponding initial net yield stress increases, indicating the suction hardening effect.The experimental results by Qin et al. 9 and Bagheri et al. 30 strongly support such rate-dependent relationship.Moreover, the suction hardening effect can also reduce the viscosity of soils, decreasing the creep behavior.This phenomenon can be represented by the decrease of  A () with suction, as presented in Figure 3.

Extending the nonstationary flow surface theory to unsaturated soils
Based on the isotach concept, the strain rate effect on the volume change behavior of saturated soils can be analyzed under one dimensional or isotropic compression conditions.For the modeling in a general stress space, Qiao et al. 11 adopted the NSFS theory to describe the strain rate effect on the shear behavior of saturated soils.In the NSFS theory, it is assumed that the yield surface can evolve with time. 11,31,32Thus, the general yield function can be written as follows: where  ′ is the effective stress tensor,   is the vector of non-hardening parameters,   is the vector of time-independent hardening parameters and  is time.
In elastoplastic constitutive models for soils, plastic volumetric strain    is generally taken as the hardening parameter.However, for viscoplastic models, viscoplastic volumetric strain    is considered as the hardening parameter.To reproduce the time-dependent deformation, Equation (7) adopts  as a variable.Based on the previous validation of isotach approach for both saturated and unsaturated soils, a unique stress-strain-viscoplastic-strain-rate relationship has been defined.Hence, by using the two variables: viscoplastic volumetric strain    and viscoplastic volumetric strain rate ε  , Equation ( 8) is obtained, as proposed by Qiao et al. 11,32 Because, at a given strain increment d   , the time increment d can be calculated by d   ε  , which means that the viscoplastic volumetric strain rate ε  is correlated with time  and can also be regarded as a time variable.In this case, the rate-dependent behavior can be presented more directly.
For unsaturated soils, the unsaturation effect on the mechanical behavior is similar to the strain rate effect.The results from CRS tests support such statement, 9,30 as shown in Figure 4.The initial net yield increases with the increase of suction.Moreover, the slope of a compression line at a larger suction is smaller than that at a lower suction, as shown by Alonso et al. 33 To describe such suction effect on yield stress, Equation ( 9) proposed by Alonso et al. 33 is adopted.Meanwhile, linking the changes in p0 to    (Equation 10), the evolution of yield stress with suction and plastic volumetric strain can be described.
where  is elastic parameter, () is the slope of compression line at suction ,   is a model parameter,    is the plastic voluteric strain,  0 is the initial yield stress at saturated state, p and p0 are the net yield stresses at unsaturated and saturated states, respectively.Combining Equations ( 9) and (10) with Equation ( 5), the following equation is obtained which enables the description of the change of net yield stress with suction, viscoplastic volumetric strain and viscoplastic volumetric strain rate.
Incorporating Equation ( 11) into the NSFS theory, a unified viscoplastic model in the general stress space is formulated for unsaturated soils, as shown in Equation (12).
where  is the net stress tensor.Figure 5. schematically illustrates the coupled effects of viscoplastic strain rate and suction on the yield surface in the p-q plane.At a given suction, the range of viscoplastic yield surface at the volumetric strain rate ε ,2 is larger than that at ε ,1 , due to a larger net yield stress.Meanwhile, the viscoplastic yield surfaces can also evolve with changes in viscoplastic volumetric strain as the common elastoplastic yield surfaces.Besides, the suction hardening not only increases the net yield stress, but also the cohesion of soils.This is also considered in the viscoplastic yield surface using the cohesion-related variable   , as presented in Figure 5. Thus, the yield surface can change with the variations of viscoplastic volumetric strain rate and suction, allowing modeling the suction-and rate-dependent properties of unsaturated soils.Meanwhile, the creep and relaxation behaviors of unsaturated soils can also be described.
Qu et al. 13 reported that there is a strain rate threshold below which the strain rate effect on the yield stress can be ignored.This is adopted in the viscoplastic model by Qiao et al. 11 for saturated soils, with consideration of a final stable state.A threshold value is also used in this study to represent the strain rate below which the rate-dependent effect disappears.For simplicity, the rate threshold εℎ  is assumed to be independent on suction.In this case, the variation of net yield stress with viscoplastic volumetric strain rate can be schematically represented in different strain rate ranges, as illustrated in Figure 6.Besides, as depicted in Figure 2, all lines can intersect at one point.In terms of net yield stress, there is a maximum value when the viscoplastic volumetric strain rate reaches a certain value.Thus, the intersection point is selected as the reference point to determine the reference stress and the reference viscoplastic volumetric strain rate in Equation ( 5), which was also recommended by De Gennaro and Pereira. 21In fact, the reference point will never be reached or exceeded due to the corresponding strain rate with a value up to 1 s −1 .At this point, the strain rate is so high that the mechanical problem is not static but dynamic.However, the viscoplastic modeling in this study belongs to statics.This point is also schematically represented in Figure 6.Combining this reference state with the final stable state, an upper bound yield surface and a lower bound yield surface can be defined.The yield surface evolves in between with viscoplastic volumetric strain rate, as illustrated in Figure 7.This is consistent with the variation of net yield stress, shown in Figure 6.Finally, combining the unsaturated isotach concept, the NFSF theory is extended to unsaturated condition with consideration of the upper and lower bound yield surfaces.F I G U R E 7 Schematic representation of the variation of yield surface with viscoplastic volumetric strain rate at a given suction  0 in p'-q plane.

Model description
In this study, the viscoplastic model is formulated in the framework of unsaturated effective stress.Thus, the coupling between water retention behavior and mechanical behavior can be better reproduced.Therefore, the Bishop effective stress is adopted for the stress expression to account for the combined effect of mechanical stress and suction (Bishop 34 ) as follows: where   is the degree of saturation, (  ) is the Bishop effective stress coefficient depending on   , and  is the identity tensor.
According to the generalized effective stress expression, the mean effective stress  ′ and the deviator stress  are defined in the triaxial stress space, as follows: where  ′ 1 and  ′ 3 are the axial and lateral effective stresses, respectively,  is the net mean stress, and  1 and  3 are the axial and lateral net stresses, respectively.Equation ( 16) is adopted to describe the variation of (  ) with   , as recommended by Alonso et al. 35 This approach can ensure a smooth transition of (  ) from low to high suctions and better reproduce the contribution of water in small pores to unsaturated effective stress under high suction.Van Genuchten (VG) model is adopted to link suction and degree of saturation, as expressed by Equation (17). 36 (  ) =    (16) where  ( > 1) is a material parameter controlling the evolution of (  ) with   ;   , n and m are model parameters.
The strain increment is partitioned into two parts: an elastic part and a viscoplastic part, as follows: where ,   and   are the total, elastic and viscoplastic strain increment tensors, respectively.
In the conventional triaxial space, the corresponding volumetric strain   and deviatoric strain   are defined as follows: With the Bishop effective stress, the elastic volumetric and deviatoric strain increments can be calculated by the following equations: where K and G are the elastic bulk and shear moduli, respectively;  is the specific volume, () is the suction-related elastic parameter, and  is the Poisson's ratio.Different experimental results have revealed the effect of suction on the elasticity of unsaturated soils. 37,38Thus, Equation ( 23) is developed to describe the variation of κ with suction.This equation was validated with the experimental data on highly compacted GMZ bentonite by Wang. 39 where (0) is the elastic parameter at saturated state,  1 and  1 are the model parameters for describing the changing rate and the smallest value of (), respectively.For the hardening law, the LC (loading collapse) curve from the BBM 33 is considered.The effective yield stress contains two parts: the net yield stress part and the suction part, as presented by Equation (24).This was also recommended by Sheng et al. 40 and Le Pense et al. 41 for the elastoplastic modeling of unsaturated soils.The evolution of the net yield stress with suction is expressed by Equation ( 25), similar to that in Alonso et al. 33 p′  = p +  (  ) * 0 =  0 exp Similar to that in Section 2.2, Equation ( 27) is obtained, enabling the description of the evolution of net yield stress with suction and strain rate: For the modeling of NFSF yield surface, the modified Cam-Clay (MCC) yield function (Equation 28) is adopted, as follows: where M is the slope of the critical state line and  ′ is the friction angle.
Inside the conventional elastic domain, creep and relaxation cannot occur. 21To make the creep and relaxation possible, a gradual hardening law is used in the proposed model (Figure 8), as in Hong et al. 42 and François and Laloui. 43The concept of mapping is adopted to define the progressive state.The image point is obtained as the projection of the current stress state on the normal yield surface, as shown in Figure 8.The stress ratio  between the current stress point and the image point is defined as follows: where  ′  and p′  are the actual yield stress and the reference yield stress, respectively,  ′ and  are the actual mean effective and deviator stresses, respectively, p′ and q are the reference mean effective and deviator stresses, respectively, From Figure 8, it can be found that by definition  corresponds to the inverse of the overconsolidation ratio (OCR). 42hereby, the introduction of  in the model allows the effect of overconsolidation on the viscoplastic behavior of soils to be described.Plastic volumetric strain is generally adopted to model the evolution of .However, Rouainia and Muir Wood 44 reported that the shear effect can result in the change of soil structure, thus inducing the soil hardening/softening.In order to account for the shear effect in progressive hardening, Hong et al. 42   where  0 is the initial value of , d   is the plastic deviatoric strain increment,   is a model parameter controlling the contribution of deviatoric loading to the gradual hardening, and  is a material parameter.
The time effect on the progressive hardening can be considered by replacing the plastic strains in Equations ( 31) and ( 32) with viscoplastic strains, as follows: where    is a generalized viscoplastic strain, and    is the viscoplastic shear strain.Incorporating the progressive hardening law in the NFSF yield function, the inner yield function and the normal yield function are obtained, as expressed by Equations ( 35) and (36), respectively.The corresponding two yield surfaces are schematically illustrated in Figure 9.
A non-associated flow rule is adopted, similar to that in the BBM. 33The inner viscoplastic potential function is as follows: where   is a parameter related to the stress state under zero lateral strain.Then, the corresponding dilatancy equation is determined with Equation (37), as expressed by Equation (39).

𝑑 = 𝑑𝜀
The consistency condition   = 0 is derived, as follows: It appears from Equation ( 40) that the effects of net stress, suction, viscoplastic volumetric strain and viscoplastic volumetric strain rate are all taken into account.Besides, by correlating the current stress state with the stress state projected on the normal yield surface based on the progressive hardening law, not only the gradual hardening but also the time-dependent behavior of soils inside the conventional elastic domain can be described.

Incremental stress-strain relation and loading criteria
For the determination of the incremental stress-strain relationship at each loading/unloading step, the elastic part can be calculated using Equations ( 21) and (22).The viscoplastic strains are determined using the dilatancy equation (Equation 39) and the consistency condition (Equation 40).The viscoplastic strain increment at each step can be computed as follows: The corresponding viscoplastic multiplier d and hardening modulus  are expressed as follows: Besides, according to the consistency condition, the viscoplastic strain depends on the net stress, suction and viscoplastic strain rate.Therefore, a more general form (Equation 44) can be deduced for the viscoplastic strain calculation, as follows: where  1 ,  2 and  3 are the partial multipliers corresponding to different loading factors.By combining the elastic and viscoplastic strains, the incremental stress-strain matrix can be deduced.In addition, the inner yield surface varies with net stress, suction and viscoplastic volumetric strain rate and viscoplastic volumetric strain, indicating that the loading criterion is also related to these factors.By referring to the plastic loading concept, the viscoplastic loading criteria are expressed as follows: Moreover, with respect to the yield surface in Figure 7, soil generally cannot reach the reference state due to the limitation of strain rate.Thus, the corresponding loading criterion is not considered.While for the final stable state, the viscoplastic volumetric strain rate ε  in the yield function is replaced by εℎ  .This means that below the threshold value εℎ  , the corresponding yield surface is rate-independent.Thus, the corresponding loading criterion is only dependent on d  .Abbreviations: QL, qianjiangping landslide; GMZ, Gaomiaozi.

Determination of model parameters
The parameters in the viscoplastic model involve the elastoplasticity, the hydraulics and the viscosity.Correspondingly, for their determination, it is divided into three parts: 1.For the elastic parameters (κ[0] and ν) and the general plastic parameters (M and λ[0]), they can be determined from triaxial and oedometer tests under different suctions.These tests also allow the determination of β 1 , r 1 , β, and r based on the relationships between suction and κ or λ, which can be identified from different unsaturated volumetric compression curves.k and A d can be obtained through sensitivity analysis.2. The water retention parameters (α v , m and n) can be determined by data fitting using VG model.p c and a can be identified based on the unsaturated shear strength of soils.3. The viscous parameters ( ε , ,   0 , ε ,ℎ ,  and b) can be obtained from CRS tests at different suctions and from creep tests.At each suction, the viscous parameter can be identified according to the evolution of yield stress with volumetric strain rate, as depicted in Figure 6.Then, the viscosity-related parameters can be determined (See Qiao et al. 11 for more details).

MODEL VALIDATION AND DISCUSSION
4][25] The calibrated model parameters for the three soils are presented in Table 1.

CRS compression behavior
Qin et al. 9 performed a step-wise saturated CRS test on GMZ bentonite at variable vertical strain rates of 10 −5 s −1 -10 −6 s −1 -10 −7 s −1 -10 −5 s −1 .Figure 10A shows comparisons between the experimental results and the model predictions in terms of void ratio variation with vertical stress.The predictions are in good agreement with the experimental result, indicating the performance of the proposed model.Moreover, the evolutions of γ are also plotted in Figure 10B.The initial values of γ are significantly smaller than 1, showing the initial over consolidation states of GMZ bentonite.As the vertical strain increases, γ increases, approaching 1.As aforementioned, the change in parameter γ represents the evolution of soil over consolidation.The application of parameter γ ensures the smooth transition from the conventional elastic region to the viscoplastic one, as shown in Figure 10A.Besides, in the previous models, the effect of strain rate on the compression behavior in the conventional elastic region cannot be reproduced, while this phenomenon can be satisfactorily described with the proposed model (although the difference is small in Figure 10A).This can be explained by the application of viscoplastic strain in Equation (34), enabling the consideration of the rate-dependent effect on the deformation in low stress range (i.e., conventional elastic region).

Triaxial shear behavior
Wu et al. 6 carried out a series of triaxial tests on reconstituted Glenroy silt under a net confining pressure of 400 kPa.Three different axial strain rates (0.1, 0.01 and 0.001 h −1 ) and three different initial suctions (s = 100, 200 and 300) were considered.The model parameters of Glenroy silt can be identified based on the water retention and triaxial test results, and are presented in Table 1.
Figure 12 shows the experimental and simulated results for three unsaturated triaxial shear tests conducted under the axial strain rates of 0.1, 0.01 and 0.001 h −1 with an initial suction of 100 kPa.Basically, the simulations are in good agreements with the experimental results in terms of stress-strain curves (Figure 12A).As the axial strain rate increases, the deviator stress increases at a given axial strain.Besides, the volumetric strains increase continuously (Figure 12B), resulting Comparisons between experiment and model for unsaturated CRS compression tests at different suctions at a constant vertical strain rate data from Qin et al. 9 ). in the increase of saturation degree under the constant water content condition.This induced the decrease of suction down to zero.Meanwhile, the rate dependency of the change in volumetric strain is also observed, which is attributed to the effect of soil viscous behavior.It appears that for the test under the axial strain rate of 0.1 h −1 , slight softening is observed but negligible, which may be related to the slight dilatancy in Figure 12B.Wu et al. 6 reported that in triaxial tests the suctions were determined using the measured air pressure and water pressure through a saturated ceramic diskette with a high air entry value from only the bottom of the sample, not representing the whole sample.Besides, when the strain rate is high enough (e.g., test under the axial strain rate of 0.1 h −1 in Figure 12C), the required equilibrium time for the water transfer in sample may conflict with the imposed strain rate, which may further result in the measured suction is not accurate.In spite of that, the rate dependency of suction change during shearing is reproduced by the proposed model in Figure 12C.
Figure 13 presents the simulations of three triaxial shear tests on Glenroy silt at an initial suction is 200 kPa with different axial strain rates.As the axial strain rate increases, a higher deviator stress, a smaller volumetric strain and a higher suction at a given axial strain.The comparisons in Figure 13A shows that the proposed model is able to well reproduce the stressstrain relationships at different axial rates without insignificant softening.In addition, neglecting the slight dilatancy at the strain rate of 0.1 h −1 in Figure 13B, the simulated evolutions of volumetric strain show a good rate dependency, in agreement with the experimental results.Compared to the test result at suction 100 kPa, the test result under the axial strain rate of 0.1 h −1 at suction 200 kPa shows more significant dilatancy and larger difference with the simulation.This is also attributed to the aforementioned relationship between equilibrium time and water transfer in the sample.Under the axial strain rate of 0.1 h −1 , the suction in experiments did not reach the real equilibrium state, as observed in terms of the measured suction in Figure 13C.In return, this induced an inaccurate volume change measured in Figure 13B.Nevertheless, the rate dependency and the suction decrease in the triaxial tests are also satisfactorily simulated by the proposed model, as shown in Figure 13C.
Figure 14 shows the experiments and simulations for three unsaturated triaxial shear tests performed under different axial strain rates and at a constant water content with an initial suction of 300 kPa.As the axial strain rate increases, the deviator stress increases, the volumetric strain decreases, and the suction also decreases for a given axial strain.The evolution trends of deviator stress, volumetric strain and suction were reproduced in simulations.The experimental stress-strain curve in Figure 14A show a significant softening phenomenon at the strain rate of 0.1 h −1 .As before, this may be related to the local transfer of pore water between the micropores and macropores at high strain rates, 6 which can further affect the evolutions of volumetric strain and suction at high strain rates, as shown by the experimental results in Figure 14B and C. Furthermore, comparing the triaxial experimental and simulated results at different suctions in Figures 12, 13 to 14, it appears that with increasing suction, the differences between experiments and simulations become larger.The reason is that the triaxial tests were carried out at a constant water content.In that case, as suction increased, the viscosity of soil increased, and the unsaturated water permeability decreased.As a result, the required equilibrium time for the water transfer under the same strain rated becomes longer, inducing more difficulties in accurately measuring the volume and suction changes.Apparently, this induces the more significant differences between experiments and simulations when increasing the initial suction.However, considering the triaxial test as an element test, the effect of water transfer cannot be reproduced in constitutive modeling.Nevertheless, on the whole, the proposed model can satisfactorily reproduce the rate-dependent shear behavior of Glenroy silt.

Triaxial creep behavior
The triaxial creep results of Qianjiangping landslide (QL) soil reported by Lai et al. 23,24 and Zou et al. 25 are considered in this section.The first set of tests by Lai et al. 24 were carried out under the same net confining pressures of 200 kPa and the same suction of 100 kPa.The samples were loaded to different deviator stress levels (i.e., 0.317, 0.489, 0.617 and 0.746 of maximum deviator stress q f = 350 kPa), respectively.For the second set of tests, the triaxial samples were loaded to the same deviator stress level (i.e., 0.55 of q f = 187 kPa, 256 and 370 kPa) under the same net confining pressures of 100 kPa and different suctions of 50 kPa, 150 and 300 kPa, respectively, as reported by Lai et al. 23  Figure 15 shows that comparisons of axial strain between experiments and simulations under different deviator stress levels at the same suction and confining pressure for the first set of triaxial creep tests.As the time increases, the axial strain increases, and a higher deviator stress level results in a larger axial strain.The proposed model can satisfactorily reproduce the evolution of axial strains with time.Compared to the axial strains in experiments, the simulated ones increase faster at the beginning, then match the experimental ones after a long term.This may be related to the loading procedure in the triaxial creep tests-the triaxial samples were fast loaded to a designed deviator stress.In this case, the higher the deviator stress, the larger the excess pore pressure.This induced a slower increase of axial strain with time at the beginning in the case of high deviator stress (e.g., ∕  = 0.617 and 0.746).Besides, as discussed in the simulation of CRS tests, the proposed viscoplastic model can reproduce the time dependent deformation in low stress range.This can explain the increase of axial strains with time under the low deviator stress levels (e.g., ∕  = 0.317), because the stress states were always located on the inner yield surface in the loading process, as presented in Figure 9.The normal yield surface evolves with time due to the change of viscoplastic volumetric strain rate, which in turn resulted in the time dependent deformation at the current stress state.
In Figure 16, the experimental and simulated results of axial strain evolutions at the same deviator stress level and confining pressure but different suctions are presented.As the suction decreases, the axial strain is larger, showing the suction-related time-dependent deformation.On the whole, the proposed model can well describe the unsaturated creep behavior.

F I G U R E 1 4
Comparisons between experiments and simulations for the triaxial tests under different axial strain rates at an initial suction of 300 kPa: (A) stress-strain curve; (B) volumetric strain change; (C) suction change (data from Wu et al. 6 ).

CONCLUSIONS
In this study, a viscoplastic constitutive model for unsaturated soils was proposed based on the nonstationary flow surface theory.The Bishop effective stress and the VG water retention model were incorporated to represent the effect of suction on the viscoplastic behavior of soils.The isotach approach for saturated soils was extended to suction effect by considering the effect of suction on the viscous parameter   which represents the variation of yield stress with the viscoplastic volumetric strain rate.An unsaturated nonstationary flow surface theory was developed, which can well describe the flow of yield surface with the changes in viscoplastic volumetric strain and viscoplastic volumetric strain rate accounting for the suction effect.Besides, the application of a gradual hardening parameter γ allowed the smooth transitions from the conventional elastic region to the viscoplastic one and the reproduction the time dependent deformation in low stress range.
The proposed unsaturated viscoplastic model was validated based the experimental results from unsaturated CRS tests under constant suctions, unsaturated triaxial tests under constant water content and unsaturated triaxial creep tests under constant suction.Comparisons between measurements and simulations showed the good performance of the model.More experimental data under other stress paths are needed for better validating the proposed model, for instance, the tests on the time-dependent behavior under the change of suction and the ones on the relaxation behavior under the constant suction for unsaturated soils.Besides, in order to better reproduce the shear behavior of unsaturated soils, the present F I G U R E 1 5 Comparisons of axial strain evolutions between experiments and simulations for the triaxial creep tests at a suction of 100 kPa (data from Lai et al. 24 ).

F I G U R E 1 6
Comparisons of axial strain evolutions between experiments and simulations for the triaxial creep tests under different suctions with the same q/q f of 0.55 (data from Zou et al. 25 ).
model can be implemented in a finite element method platform by adopting a permeability law.Thus, water transfer during shearing and its effect on softening and dilatancy can be reproduced through numerical modeling, but cannot reproduced through constitutive modeling.

D ATA AVA I L A B I L I T Y S TAT E M E N T
The data in this study are available from the corresponding author upon reasonable request.

F I G U R E 4
Schematic illustration of the isotach concept for unsaturated soils.

F I G U R E 5
Schematic illustration of yield surfaces in p-q plane.

F I G U R E 6
Schematic representation of the variation of net yield stress with viscoplastic volumetric strain rate after (after De Gennaro & Pereira21 ).

F I G U R E 8
Inner and normal yield surfaces.

F I G U R E 9
Schematic illustration of the yield surfaces.

F I G U R E 1 0
Model validation with saturated CRS compression tests at different constant strain rates: (A) comparisons between experiments and simulations (data from Qin et al. 9 ); (B) evolutions of γ.

Figure 11
Figure11presents the comparisons between measurements and simulations for the CRS compression tests at different suctions.As the suction increases, the yield stress increases for a given vertical strain.The simulated compression curves are in good agreement with the experimental ones, showing the performance of the model in describing the unsaturated volume change behavior.Similarly, the application of parameter γ contributes to the better simulations.In addition, it can be seen from Figure11that the experimental compression curves have different initial slopes at different suctions, showing the suction dependency of elastic parameter .This suction effect is considered through Equation (23).

F I G U R E 1 2
Comparisons between experiments and simulations for unsaturated triaxial tests at different axial strain rates with an initial suction of 100 kPa: (A) stress-strain curve; (B) volumetric strain change; (C) suction change (data from Wu et al.6 ).

and Zou et al. 25 F I G U R E 1 3
Comparisons between experiments and simulations for the triaxial tests at an initial suction of 200 kPa and with different axial strain rates: (A) stress-strain curve; (B) volumetric strain change; (C) suction change (data from Wu et al.6 ).
Model parameters.
TA B L E 1