A comparative study of different model families for the constitutive simulation of viscous clays

Correspondence Merita Tafili, Karlsruhe Institute of Technology, Institute of Soil Mechanics and Rock Mechanics, Engler-Bunte-Ring 14, D-76131 Karlsruhe, Germany. Email: merita.tafili@kit.edu Summary The simulation of the viscous behavior of some clays is of high importance in many geotechnical problems. The literature offers a vast amount of constitutive models able to simulate the rate dependence observed on these materials. Although most of these models are calibrated to very similar experimental observations and share similar definitions of material parameters, some discrepancies of their response have been detected, which are related to their mathematical formulations. In this work, the causes of these discrepancies are carefully studied. To that end, four different model families are analyzed, namely, nonstationary flow surface (NSFS) models, viscoplasticity with overstress function (OVP), viscoplasticity with Norton's power law (NVP), and visco-hypoplasticity (VHP). For the sake of a fair comparison, single constitutive models using the same set of material parameters, and following other requirements, are developed for each model family. Numerical implementations of the four resulting models are performed. Their response at different tests are carefully analyzed through simulation examples and direct examination of their constitutive equations. The set includes some basic tests at isotropic stress states and others as responses envelopes, undrained creep rupture, and an oedometer test with loading, unloading-reloading, creep, and relaxation. The article is concluded with some remarks about the observed discrepancies of these model families.

Many rate-dependent models for clays are actually based on very similar fundamentals: the basis of critical state soil mechanics, 34,35 accounting for the stress and void ratio dependence on stiffness and strength characteristics, and their viscous effects capabilities, adjusted to some empirical relations. In addition to this, their mathematical structure follows mostly from viscoplastic formulations, ie, the irreversible strain rate, and not, for example, the stiffness tensor, is responsible for viscous effects. With the aforementioned similarities, one may believe that two different viscous models whose formulations are adjusted to identical critical state relations, and to the same empirical relations for viscous effects, should deliver the same response. This is, however, not true: The literature offers different theoretical approaches to develop the model formulation, presenting important discrepancies between them, of which users should be aware of. One may distinguish, for example, viscous models based on nonstationary flow surface (NSFS), [36][37][38][39][40][41] viscoplastic models based on overstress functions (OVP), 13,[42][43][44][45][46][47][48][49][50][51] others based on potential laws similar to Norton's law (NVP), [52][53][54][55][56][57][58][59] and yet recently, the so-called visco-hypoplasticity (VHP), 60 among other formulations. 39,61,62 Model developers know that it is possible to formulate viscous models with identical material parameters but having different formulations. Works showing the capabilities of single models through predictions of experimental curves are often encountered in the literature. 13,60,[62][63][64][65] Nevertheless, very few works are devoted to compare different formulations of viscous models, 66,67 and to the authors' knowledge, a work showing their differences directly from simulation examples is not found in the literature, probably because it would demand an exhaustive work for their formulation under the same conditions and for their numerical implementation.
The present review article aims to show the differences emerging from different model formulations for viscous clays directly from simulation examples. To that end, four different models will be proposed and implemented according to the following model families: NSFS, OVP, NVP, and VHP. In order to provide a fair comparison, they will adjusted to the same empirical relations and calibrated with the same parameters. Most simulations presented herein do not aim to show the congruence of the models with a particular material but rather to provide simulation examples to analyze their differences due to their mathematical formulations. The selected set of model parameters follows from well-known modified Cam clay (MCC) and Buisman-type 27 relations. The article is structured as follows: first of all, some overall concepts often used in viscous constitutive formulations are introduced. Then, the constitutive equation for each model family is developed. These equations are adjusted to the same set of material parameters. Subsequently, the material parameters and the numerical implementations are briefly described. Several tests showing rate-effects are then simulated and carefully analyzed with the different models. The set of simulations include isotach curves, creep tests, stress relaxation tests, response envelopes, and undrained creep tests. This analysis is complemented with the simulation of an oedometric test of a Kaolin clay showing many rate effects. The article is closed with some final remarks.

NOTATION AND CONVENTION
The notation and convention of the present work is as follows: Italic fonts denote scalar magnitudes (eg, a, b), bold lowercase letters denote vectors (eg, a, b), bold capital letters denote second-rank tensors (eg, A, ), and special fonts are used for fourth-rank tensors (eg, E, L). Indicial notation can be used to represent components of tensors (eg, A ij ), and their operations follow the Einstein's summation convention. The Kronecker delta symbol is represented by ij , ie, ij = 1 when i = j and ij = 0 otherwise. The symbol 1 denotes the Kronecker delta tensor (1 ij = ij ). The unit fourth-rank tensor for symmetric tensors is denoted by I, where I i kl = 1 2 ( ik l + il k ) . Multiplication with two dummy indices (double contraction) is denoted with a colon "∶" (eg, A ∶ B = A ij B ij ). The symbol "⊗" represents the dyadic product (eg, A ⊗ B = A ij B kl ). The brackets || ⨆ || extract the Euclidean norm (eg,

CHARACTERISTIC ISOTACHS AND OVERCONSOLIDATION RATIO
In classical soil mechanics, isotropic compression yields to a unique e-p curve under normally consolidated states. 29 For viscous soils, a shift of this e-p line can be obtained depending on the strain rate velocity. Similar to that of other authors, 7,30,62,68 the term "isotach" is hereafter used to denote a characteristic e-p curve related to a particular volumetric strain rate The e − p space wherein isotachs exist is bounded according to some works. More specifically, some models consider the existence of a lower boundary at which isotachs of extremely low strain rates . v ≈ 0 are bounded, or equivalently, at which creep paths cease (eg, previous studies 43,44 ). This bounding curve is in this work referred to as "minimum" isotach. Other works consider a "maximum" isotach as the upper boundary of isotachs with very high strain rates . v ≈ ∞ (eg, previous studies 38,39,41 ). These and other characteristic isotachs are often used for the formulation of many existing viscoplastic models. They are summarized in the following lines: Reference isotach: isotach for . v = D r , whereby D r is a . For the sake of convenience, the reference isotach is usually set to coincide with an experimental curve. In this work, for the case of the OVP model, the reference isotach is selected to coincide with the minimum isotach. Maximum isotach: isotach for the limit value . v → ∞. It corresponds to the upper boundary of all possible isotachs in the e − p space. Minimum isotach: isotach for the limit value . v → 0. It corresponds to the lower boundary of all possible isotachs. It also bounds creep and relaxation paths under isotropic stress states (q = 0). OCR=1 isotach: isotach at which the overconsolidation ratio under isotropic states is equal to one, ie, OCR= 1. This coincides with the definition of the normal consolidation line.
Notice that while the reference isotach is thought to match an experimental curve, the maximum and minimum isotachs are rather mathematical boundaries in the e − p space. Several relations proposed in the literature are useful for their formulation, eg, previous studies. 34,35,[69][70][71][72] For instance, consider the MCC relation to describe the characteristic isotachs denoted collectively by e = e x (p): where is the compression index and e x0 = {e ref0 , e max0 , e min0 , e i0 } corresponds to the characteristic void ratio of the reference, maximum, minimum, and OCR= 1 isotachs at p = p r = 1 kPa respectively; see Figure 1A.
One may show that the studied models in the present work can be adjusted to a reference isotach for a compression curve with . v = D r . However, the incorporation of a maximum or a minimum isotach depends on the model type and is actually not mandatory. Particularly, the NSFS model and the VHP model incorporate only the maximum isotach in their formulations. This curve has been previously referred to as "instant normal compression line" by Yao et al 37 or "maximum void ratio isotach" by Fuentes et al. 60 OVP models incorporate a minimum isotach, referred to as "static yield surface" in previous studies. 13,73,74 The NVP model is the only one incorporating neither the maximum nor the minimum isotach. It will be shown that many salient discrepancies of their response emerge from the mentioned differences.
The overconsolidation ratio in one dimension is defined as whereby p ei corresponds to the Hvorslev pressure at the OCR= 1 line, solved from Equation (1). For the sake of generality, we introduce the ratio R x , being similar to the OCR, but considering other characteristic isotachs as reference curves at which R x = 1: Figure 1B illustrates the Hvorslev-type variables p e,x at the characteristic isotachs. The definition of the overconsolidation ratio for three-dimensions, denoted hereafter with OCR 3D , has been differently proposed by some works, or actually omitted in others. Most authors propose a function describing a surface in the stress space at which OCR 3D = 1, represented by the condition F OCR = 0. To keep consistency, this surface intersects the equivalent Hvorslev pressure p ei at the isotropic axis, ie, For the sake of simplicity, MCC-type yield surfaces will be used for the formulation of the OCR 3D surface of the OVP, NVP, and NSFS models. The following MCC function F OCR = F e is adopted for these models: where M is the critical state slope in the p − q space. To make our analysis simpler, we consider M a constant. However, it is strongly recommended that a Lode's angle dependence be incorporated as in other works. 60 The resulting OCR 3D definition being consistent with Equation (2) and (5) is , (MCC-type), whereby p + has been solved from Equation (5) for F e = 0. An analogous relation is also used for the ratio R x under three-dimensional states: OCR 3D functions based on different relations than MCC-type functions have been used in the literature as well. The VHP model in previous studies 60, 75 proposed a stress surface for the OCR with the following features: an open-wedge type surface that intersects the isotropic axis at p = p ei and the critical state surface at p = p ei ∕2. The function proposed by Fuentes et al 60 reads where M is the critical state slope in the p-q space, and f b is a scalar function defined as and f b0 is considered a material constant controlling the slope of the surface at p → 0 in the p − q space. Accordingly, a value of approximately f b0 ≈ 1.1 has been found to reproduce well many experimental observations. 60,75 Hence, we set this value as a constant for the simulations. Figure 2 illustrates the OCR 3D stress surfaces for MCC-type functions and the VHP model by Fuentes et al. 60

CREEP RATE
Under isotropic stress conditions (q = 0), the creep rate describes the rate of volumetric deformation . v upon a constant stress state . = 0. Among many relations proposed in the literature, 7,27,28,33 a Buisman-type relation 27 is selected for this whereby v,0 and t 0 are initial values of the void ratio and time, respectively, and C is analogous (but not equal) to the Buisman constant. 27 The differential form of Equation (10) is Considering that creep paths are performed under constant mean pressure . p = 0 while a change on the Hvorslev stress p ei is experienced, the difference − v,0 is then related to v − v,0 = − ln(p ei ∕p) = − ln(OCR). One can demonstrate that substitution of the latter relation in Equation (11), and assuming that creep volumetric strains are produced by viscoplastic strains . vis v , yields to the following relation: where D r = C ∕t 0 is a reference creep rate, already introduced in other works (eg, Niemunis 55 ), and I v is, by some authors, referred to as the viscosity exponent or also the Leinenkugel index, 76 which reads Note that D r = C ∕t 0 depends on the reference time t 0 = 1 s. This suggests that D r is a referential factor with units of 1∕s. It can be adjusted to a particular isotach, as for example, the reference isotach. The three-dimensional extension for Equation (12) incorporates the flow rule m, first referenced in Niemunis, 55 as follows: where the norm ||m|| = 1∕ √ 3 has been imposed to match the behavior for 1D case in terms of . p and . v .

DEVELOPMENT OF CONSTITUTIVE EQUATIONS
With the aim of performing a fair comparison, the viscous models are required to account for a set of characteristics, which are established in the following lines: 1. Isotropic compression under the reference strain rate, ie, Poisson ratio is used as parameter for this purpose. 5. Under isotropic stress conditions, a point lying at the reference isotach exhibits a creep rate equal to The viscosity index I v is introduced as a parameter for this purpose.
The aforementioned requirements are related to the usage of the set of parameters {e ref,0 , , D r , , M, , I v } in the constitutive formulations. In the following sections, the different formulations are presented considering the aforementioned requirements.

NSFS model
NSFS models are based on two main features: the formulation of a time-dependent yield surface responsible for the creep behavior and the consideration of the maximum isotach as the OCR= 1 isotach, ie, The model will be based on some MCC relations, as in other NSFS models in some previous studies. [77][78][79] The Hvorslev pressure p ei serves as hardening variable of the yield surface in Equation (5). It presents an isotropic hardening rule as in conventional MCC theory 35 with the following relation: , whereby e 0 is the initial void ratio, p ei0 = p e,max0 is the initial Hvorslev pressure at the maximum isotach, e i0 = e max0 is the void ratio of the maximum isotach at p = p r = 1 kPa, and vp v is the viscoplastic volumetric strain. An additional term f t = f t (t) is added in Equation (5) in order to account for the time dependence of the yield surface. Following this, the time-dependent yield surface function reads Combination of Equations (16) and (17) in conjunction with Equation (15) yields to the following relation: Differentiation of Equation (18) leads to the consistency condition . F e = 0: .
In order to find an expression for . t , the creep condition . = 0 is imposed at the consistency condition . F e = 0. By doing this, the following relation is obtained: Substitution of Equations (12) and (13) in Equation (20) gives The integrable form of the latter equation is t = ∫ According to the Kuhn-Tucker conditions, the model is elastic for F e < 0 and elasto-viscoplastic for F e = 0, ie, for F e < 0, where . vp is the viscoplastic strain rate tensor and E denotes the (hypo-)elastic (fourth rank) tensor, which reads where 1 is the Kronecker delta tensor, K is the Bulk modulus, and G is the shear modulus, both having MCC relations The viscoplastic strain rate tensor . vp follows from conventional plasticity relations: whereby . vp is the consistency parameter and m is the flow rule. The consistency parameter . vp is solved from the consistency condition The solution of . vp can be decomposed according to where . pl and . vp are the plastic and viscous consistency parameters, respectively, and read This decomposition allows to rewrite the constitutive Equation (22) with the following form: where . pl = . pl m and . vis = . vis m, such that . vp = . pl + . vis . For the sake of simplicity, the flow rule is selected to be associated with the yield surface: Notice that ||m|| = 1∕ √ 1∕3. Similar to that of other authors, 80,81 we distinguish between continuum Jacobian, defined as J = ( . )∕( . ), and algorithmic Jacobian, defined as J = ( Δ )∕( Δ ). While the first depends on the model formulation, the second depends on its numerical implementation algorithm. One may demonstrate that the continuum Jacobian J of the NSFS model is similar as in classical elastoplasticity 80 : where H is the hardening modulus and reads As an additional remark, notice that this model is able to simulate an inviscid clay by setting I v = 0 because . vis = 0; see Equation (27). For that case, the normal consolidation line would coincide with the maximum isotach.

OVP model
According to the OVP, the total strain rate . splits into an elastic component . e and a viscoplastic component . vp , 78,[82][83][84][85] that is, Then, the general constitutive equation reads: The formulation of the (hypo-)elastic stiffness E follows from MCC relations and is identical to Equations (23) and (24). The following lines are now devoted to deduce the viscoplastic strain rate . vp . Assume that the OCR = 1 isotach coincides with the minimum isotach, ie, OCR= R min and e i0 = e min0 . Considering this assumption, the yield surface function under isotropic states (q = 0) reads F e ≡ ln(p∕p e,min ) = 0 for q = 0.
We now demonstrate that the function exp (F e ∕I v − 1) is proportional to the viscoplastic strain rate 1∕I v under isotropic stress states (q = 0); see Equation (12). Evaluation of function exp (F e ∕I v − 1) yields to the following relation: For typical values of 1∕I v > 10, one may show that when p > p e,min , the approximation ln (p∕p e,min ) 1∕I v − 1 ≈ ln (p∕p e,min ) 1∕I v holds. Using the mentioned approximation in Equation (35), in conjunction with OCR=R min , gives The last equation proves that exp In light of this demonstration, OVP models adopt the following relation for . vp : where D r is a parameter described in Section 4 and is often referred to as the "functional of overstress" 44 or the "viscous nucleus". Once more, an associated flow rule has been adopted and reads Considering that the viscoplastic strain rate . vp does not depend (explicitly) on the strain rate . ilon, the continuum Jacobian J of the OVP model simplifies to

Visco-plastic model with Norton's power law
Viscoplastic models with Norton's power law are very often used in the literature. 7,13,52,55,59,62,86 The main difference of their formulation with respect to the others is the fact that they do not incorporate any boundary isotach. Its formulation is very simple and is explained in the following lines. The general constitutive equation reads whereby E follows the MCC relations according to Equations (23) and (24). The viscoplastic strain rate vp is proposed to be proportional to the Norton's power law, yet demonstrated in Equation (14) : The overconsolidation ratio is determined through the relation given in Equation (6). Note that for highly overconsolidated states OCR≫ 1, the viscoplastic strain rate is negligible . vp ≈ 0, and the model becomes practically (hypo-)elastic.
Once more, an associated flow rule m is adopted: Similar to the OVP model, the continuous modulus is equal to the elastic stiffness tensor, ie, J = E.

VHP model
VHP models assume a strain rate decomposition into three strain rate components as follows: whereby . e is the elastic strain rate, . hp is the hypoplastic strain rate, and . vis is the viscous strain rate. In contrast to other formulations, these three components are always active and are not subjected to any yield constraint. Accordingly, the constitutive equation reads vis is the viscous stress rate and E is the (hypo-)elastic stiffness, defined as The factors K and G are the bulk and shear modulus, and r = * ∕p is the stress ratio. The bulk modulus K and shear modulus G are the same presented by Fuentes et al, 87 adjusted to match the MCC relations: The strain rate components . hp and . vp are defined by the following relations: with Y being a scalar function termed the degree of nonlinearity, 55 and the overconsolidation ratio OCR 3D according to Equation (8). The degree of nonlinearity Y reads where r c = √ 2∕3M c − → r . One may notice that the maximum isotach with OCR= R max and therefore e i0 = e max0 is implemented. Finally, the visco-hypoplastic flow rule m reads Notice that one may simplify the constitutive model to different formulations according to the problem of interest. For example, for models specializing on monotonic loading without viscous effects, one may use the assumption I v = 0, which reduces the above model to a hypoplastic model for inviscid clay. The continuous modulus J is computed similar to conventional hypoplastic models:

ADJUSTMENT TO A UNIQUE REFERENCE ISOTACH
One of the requirements of the present study is that all models are adjusted with identical reference isotach. Let the reference isotach be the curve described by the equation Differentiation of Equation (55) in conjunction with the condition Let us now simplify the constitutive equations of the NSFS, OVP, and NVP models, for isotropic stress states q = 0 subjected to a reference strain rate compression . v = D r . According to Equations (22), (33), and (42), the following simplified forms of these models for q = 0 and The particular case of the Visco-hypoplastic model is treated apart considering the nonlinearity of the hypoplastic strain rate . hp with the strain rate, ie, . hp ∼ || . ||. Its simplified form under the aforementioned conditions is Substitution of vis v , K, and Y 0 with their respective definitions (see Section 5) yields to the following relations: The resulting equations depend on the OCR at the reference isotach. Notice that the value of OCR can be directly solved for the OVP and NVP models, from Equations (60) and (61), respectively, as follows: Hence, the normal consolidation line (OCR = 1) corresponds nearly to the reference isotach for the OVP and the NVP model. In contrast, simple numerical methods are required for the solution of OCR for the NSFS and VHP models, ie,   (59) and (62), respectively. Once the OCR is solved, the position of the OCR = 1 isotach of the different models is adjusted through parameter e i0 , which can be solved from the condition With this method, all models are adjusted to have the same reference isotach. Note that the normal consolidation line of the NSFS and the VHP model corresponds to their maximum isotache, as also depicted respectively in Sections 5.1 and 5.4. Hence, in order to back-fit an experiment result of OCR = 1 line (corresponding to a strain rate used in the experiment), different models require different reference strain rates, obeying Equations (59) to (62). Hence, instead of OCR, Equations (59) to (62) can be solved for D r and then the parameter e i0 corresponds to e i0 = e ref .

MATERIAL PARAMETERS
The four models were formulated with the same set of material parameters. In total, seven parameters are required for their calibration. They are subdivided in those related to the MCC theory and those describing the viscous behavior. The MCC parameters correspond to compression index , the swelling index , the void ratio at the reference isotach e ref0 , the critical state slope M in the p − q space, and the Poisson ratio . The additional two parameters controlling the viscous behavior correspond to the viscosity index I v and the reference strain rate D r . Detailed description of the procedure for their determination is encountered in Hadzibeti 88 and other works. However, in Appendix A, some hints are given for their calibration. Table 1 lists the required parameters, with their respective names, approximate range, and some suggested tests for their determination. In summary, two isotropic compression tests with different strain velocities and one triaxial compression test are required as minimum to calibrate all parameters. Obviously, a higher number of tests increases the calibration quality.

NUMERICAL IMPLEMENTATION
All constitutive models were implemented in a FORTRAN subroutine as "User Material" (UMAT), compatible with the software ABAQUS Standard. 89 A substepping scheme was adopted for all implementations to avoid numerical issues. Very small substepping sizes in order of ||Δ || = 10 −5 were selected to assure numerical convergence. In addition, some numerical recommendations for viscous models pointed in previous studies 55,62,90 were considered for the OVP, NVP, and VHP models. Specifically, while the viscous strain rate was semi-implicitly integrated, other components of the models were explicitly computed. This required the computation of an algorithmic Jacobian J alg , which is explained in the sequel. Numerical integration of the models under different paths assuming element test conditions (homogeneous field for stresses and strains) was solved with the use of a Newton-Raphson solution scheme, whereby tolerances and a maximum number of iterations were specified to minimize the numerical error. For the NSFS model, an elastic predictor was always computed to evaluate if an elastic or rather a viscoplastic step should be performed, as in Simo and Hughes. 91 At each viscoplastic step, a small correction was applied to the stress increment Δ to assure the condition F e = 0 at the end of the increment. Other models were directly implemented without distinguishing between elastic and viscoplastic steps due to the lack of a yield surface. All implementations showed rapid convergence, and no numerical issues were detected. All the subroutine UMAT are now available for their use in finite element computations. The algorithmic Jacobians of the different models were found as follows.
The NSFS model has been explicitly integrated with a substepping technique. Hence, the algorithmic Jacobian J alg has been set equal to the continuum Jacobian J (see Equation 30), ie, J alg = J. For the OVP and NVP models, a semi-implicit integration for the viscous stress increment Δ vis has been implemented. The stress increment Δ is therefore computed according to the following mathematical procedure: where the explicit stress increment Δ exp is computed as The required derivatives ( . vp )∕( ) and ( . vis )∕( . ) are given in Appendix B. A similar integration scheme has been also performed with the VHP model. The only difference with Equation (66) is the computation of the explicit stress increment Δ exp , which, according to Equation (46), reads as

ELEMENT TEST SIMULATIONS
In this section, a set of element test simulations are presented to analyze the performance of the different model families.
These simulations aim to show the distinct responses of the models under different conditions. Firstly, a set of illustrative element tests are analyzed, ie, they are not calibrated to any particular material. The selected parameters are typical on some clay-like soils and are listed in Table 2. The mathematical procedure in Section 6 was followed to adjust all models to the same reference isotach, being characterized by the set of parameters e ref0 , , and D r . According to Section 6, the selection of a unique reference isotach for all models implies that they present different OCR = 1 isotachs. Therefore, the resulting characteristic void ratios e i0 , describing the OCR = 1 isotach, were computed for each particular model according to the relations given in Section 6. The values of e i0 are listed in Table 3. The illustrative simulations correspond to some tests showing strain rate dependency. They include isotropic compression tests with different strain rates, creep tests, relaxation tests, oedometric tests, and undrained creep tests. In addition, responses envelopes analysis 92 were included as well to analyze the models performance on different strain/stress rate directions. Response envelopes will allow to detect directions at which viscoplastic effects act the most. At the end, simulations of two experiments are included, the first consisting of an undrained triaxial test with tertiary creep and the second of an oedometric test with different strain rate, creep, and relaxation stages.

Isotropic compression with constant strain rate
Simulations of isotropic compression tests with different volumetric strain rates . v were performed. These simulations cataloged as basic, allow us to obtain an overview of the isotach boundaries acting at each model. Volumetric strain rates ranging between . v = {10 −1 , 10 −10 } 1∕s were systematically applied to generate different compression paths. Figure 3 presents the simulations of the NSFS, OVP, NVP, and VHP models. They show that under reference strain rate . v = D r , yield to the reference isotach e = 1.2−0.1 ln(p), as a result of considering the procedure described in Section 6. The figures show that the compression curves of the NSFS and VHP models are bounded by maximum isotachs, which are different for each model. The isotachs of the OVP and NVP models are bounded by their elastic response at very high velocities, with slope equal to . This upper bound actually depends on the initial conditions of e and p and has been previously referred to as "instant time line" by other authors. 37,66 Simulations with the OVP model show that a lower bound exists because of the consideration of a minimum isotach. Below this curve, the existence of isotachs for this model is not possible. But this does not have to be necessarily a negative point, because the lower bound can be shifted by choosing a very low reference creep rate D r .
One of the reasons to use a minimum isotache is to limit the creep deformation at an infinitely long time, even though this is not proven by experimental evidence yet, but it is not refuted either. Its definition requires an additional parameter, such as R min . The problem is that an infinitely long time test is needed for its calibration. Otherwise, special procedures can be developed in order to determine this parameter like the construction of prototype models for the considered geotechnical problem whereby the final creep rate can be derived as a function of time. Note that in the literature, it is controversial whether the creep reaches a final value or not. In the current paper, the minimum isotache for OVP simply corresponds to the reference isotache line, which is not always the case. Some applications, such as problems dealing with very long time, may be disadvantageous because the creep is limited. Yet, for issues with aging effects, the introduction of this isotache may be rendered as convenient.
The maximum isotach has been introduced into some models, because it is inferred that above the normal consolidated state OCR = 1, there must be a compression line for which the creep time is "0." 37 The line then describes the normally and instantaneously compressed behavior of clays, thus the name "instant normal compression line." For geotechnical analysis, engineers may be used to have a bound of the void ratio, which implies a bound of the stress ratio, as otherwise macro pores and their impact on the constitutive modeling should be treated. Because these effects are not subject of the models described here, we find it appropriate to introduce a maximum bound of isotachs. The attentive reader may believe that the introduction of a maximum isotache may not be sufficient for installations (shaking) or especially very fast earthquakes. Yet, both the NSFS and the VHP model reach the maximum isotach for . → ∞, hence every strain rate (even very fast loading) is considered. Figure 4 presents the resulting mean stress p at a given void ratio e = const for different strain rate velocities. In particular, Figure 4A presents the results for a void ratio of e = 1.15 and shows that while the NSFS and the VHP model have yet not reached the maximum isotach at high velocities, the OVP model has already reached the minimum isotach at low velocities. The NSFS and VHP reach their maximum isotach at a void ratio of e ≈ 1.0, as shown in Figure 4B.
The stress ratio, namely, the overconsolidation ratio OCR reached at the different isotachs, is depicted in Figure 5. The figure shows that while the NSFS and VHP models show values of p ei ∕p ≥ 1, the OVP model shows always p ei ∕p ≤ 1. The NVP model does not show any bound for the stress ratio p ei ∕p. For geotechnical analysis, engineers may be used to have a bound of the void ratio, which implies a bound of the stress ratio, as otherwise macro pores and their impact on the constitutive modeling should be treated. However, these effects are not subject of the models described here. Worthy to note is also the fact that the models normal consolidation line (OCR=1) corresponds to different strain rates, hence isotachs, for each model. While for the OVP and NSFS model the OCR = 1 line corresponds to the reference isotache, for the VHP and NSFS models, the maximum isotache corresponds to the normal consolidation line OCR=1; see also Figure 5.
Examination of the constitutive equations allow us to determine the relations describing the isotach formulation. In the following lines, relations depending on the void ratio e, mean pressure p, and strain rate . v are found to describe the From Equation (69), one can show that from the constant strain rate condition̈v = 0 results the following relation: where the subindices a and b denote two different isotachs. Substitution of the corresponding definitions of . vp v for the NVP and OVP models in Equation (70) leads to the relations . Similarly, the isotach condition for the VHP model is Substitution of the corresponding definitions of K and . hp v in Equation (73) results in the following relation: Hence, two isotachs are related with the following relation for the VHP model: .
The inequality x a ≠ x b is true considering that x = x(p, e) is not a constant but rather a function of the mean stress p and void ratio e. Finally, we inspect the NSFS model. Considering the NSFS relations x Substitution of One can show that for the NSFS model, the relation x = ( ∕ −1+ . pl ) holds. Equations (71) It should be emphasized that most existing one-dimensional isotache models tend to overestimate the creep strains measured after surcharge loading due to their viscous formulation = 0.434 C ln( Yuan et al 93 proposed a two-parameter modification of this relation, which has resulted in a better match with laboratory data. Because of the potential dependence of the viscous strain rate on the overconsolidation ratio (see Equation 14), the creep settlements after surcharge loading may be predicted well by these models, as will be shown in Section 9.2. Furthermore, the viscous behavior of highly overconsolidated clays, ie, OCR ≥ 4, is likely to be dominated by secondary swelling instead of creep. But certainly, the creep may resume after for example a decade-long wait.

Creep test
Creep tests under isotropic states (q = 0) are now performed to evaluate the models response. The creep condition is simulated by restraining the effective stress for a large time, ie, . p = 0 and q = 0. Following Section 4, a similar response to the one described by the Buisman-type Equation 10 is expected from the simulations. For each model, two creep tests after isotropic compression were performed: Their initial conditions correspond to e 0 = 1.2, p = 1 kPa. An isotropic compression is applied but with distinct strain rates . v = {D r , 5D r } for each test. In order to achieve the same creep response, D r is chosen according to the values listed in Table 2, thus D r = 10 −10 for OVP and D r = 10 −6 for NSFS, NVP, and VHP. At the end, a creep path of 10 10 seconds is performed. Figures 6 and 7 present the simulations of all models. The response of the models can be summarized as follows: A similar slope of the creep test can be observed for all models in the e vs ln(t) space, see Figure 7, except by the OVP model. The slope is equal to C in the e vs ln(t) space, see The accurate simulation of creep is indispensable for predicting the flow pressure on dowelling elements of creeping slopes. The reduction of creep with OCR proves that the models can capture well the surcharge-induced stress history of the material. Hence, they may provide a realistic creep prediction and are able to control long-term settlements using surcharge loading.
Another, very important geotechnical issue in engineering practice is the prediction of long-term consolidation settlement, especially the process of residual settlement. [94][95][96] In construction site, the consolidation settlement is generally estimated based on the e − ln(p ei ) line obtained from a 24-hour incremental loading oedometer tests, which often corresponds to a strain rate of . = 10 −7 /s. In order to achieve a more accurate prediction, Watabe and Leroueil 96 derived a relation for the creep strain using the compression index C c at the consolidation pressure, the lower limit of the preconsolidation pressure, and the preconsolidation pressure at the respective strain rate. However, in situ strain rates are several orders smaller than those used in experiments. Considering Equation (11) in one dimension, which presents the basis for the multiaxial creep rate incorporated into the models, it can be concluded that one-dimensional compression creep can be captured well by the models as a function of in situ strain rate, C , C c , and void ratio e.

FIGURE 6
Creep tests after isotropic compression, e vs p space. D r corresponds to the values listed in Table 2, thus D r = 10 −10 for OVP and D r = 10 −6 for NSFS, NVP, and VHP

Stress relaxation test
Stress relaxation paths are now simulated under isotropic states (q = 0). The stress relaxation condition is achieved by restraining deformations under isotropic states, ie, . = 0 with q = 0. Similar to the creep test, the stress relaxation is started after an isotropic compression path. For each model, two relaxation paths were simulated after isotropic compressions with strain rates of In order to achieve a similar relaxation response, D r is chosen according to the values listed in Table 2, thus D r = 10 −10 for OVP and D r = 10 −6 for NSFS, NVP, and VHP. Both relaxation paths start at different isotachs but with same void ratio e = 0.8. They simulate a period of t = 10 10 seconds. The results are plotted in Figures 8,  9, and 10. Figure 9 shows that for a given void ratio, both simulations of each model yield asymptotically to the same curve in the p vs ln(t) space. The analytical solution of this asymptote is deduced later for all models except for the OVP model and is plotted with dashed lines in Figure 9A,C,D. Notice that the OVP model path is stopped by the minimum isotach, which acts as lower boundary. Of course, by choosing , the relaxation behavior of OVP is more pronounced and similar results to NVP are obtained. In Figure 10, the curves are gathered together for comparison purposes.
We now proceed to deduce the analytical solution of the stress relaxation paths. Consider the NVP model constitutive Equation (42). The relaxation condition . = 0, together with Equation (42) yields to which is an ordinary differential equation with solution p = p(t). Let c nvp = c nvp (e) denote the function of the void ratio ∕ . Substitution of c nvp in Equation (78) leads to the simplified form where L 1 is an integration constant, for instance, found from the initial relaxation stress condition with p(t = 0) = p 0 : Substituting L 1 into Equation (79) yields to the relation A similar mathematical procedure can be followed to find the analytical solutions for the NSFS and VHP models. Consider the constitutive Equations (22) and (46) under the conditions q = 0 and . = 0: One may show that the same differential Equation (79) Table 2, thus D r = 10 −10 for OVP and D r = 10 −6 for NSFS, NVP, and VHP VHP model. Following the same procedure, one can show that the simplified differential form describing the relaxation process (q = 0, ∕ is a void ratio function. The latter equation lacks of analytical solution. Figure 9 includes the analytical solutions of the NSFS, NVP, and VHP models.

Response envelopes
Response envelopes were originally proposed by Gudehus 97 to analyze the directional response of constitutive models. The shape of response envelopes in the p-q space results from joining the end points of different p − q paths, all of them beginning with the same initial conditions (e = e 0 , = 0 ) and subjected to identical strain increment magnitude ||Δ || = 0.2 % on different strain rate directions − → . . The analysis is usually performed under triaxial conditions. For each model, four different response envelopes are constructed under triaxial conditions. Two initial conditions were considered: at isotropic stress states, with p = 100 kPa, q = 0, and at anisotropic stress states at p = 100 kPa and q = 75  Table 2, thus D r = 10 −10 for OVP and D r = 10 −6 for NSFS, NVP, and VHP kPa. For each initial state, two different strain velocities were applied, || . || = 10 −5 1/s and || . || = 10 −7 1/s. The results are plotted in Figure 11. From the results, the following observations are commented: • Response envelopes showed larger stress increments Δ for the lowest strain rate case (|| . || = 10 −7 1/s). They showed also that the effect of the strain rate is only evident in a particular zone of the response envelope. The zone showing no strain rate dependency is governed by the elastic response of the model. • The NSFS and VHP models showed lower strain-rate dependency than other models. This is related to the fact that both models incorporate the maximum isotach and the selected initial conditions are close to the maximum isotach. Therefore, their viscous effects are in someway bounded. • In contrast to other models, response envelopes of the NSFS model present corners at points dividing the elastic and viscoplastic response. Other models showed rounded shape. Corners on responses envelopes are typical on elastoplastic formulations and coincide with the division of the elastic and (visco-)plastic response. 55 To the date, these corners have not been experimentally found. • The response envelope of the VHP model with initial conditions of p = 100 kPa and q = 75 kPa shows a "rotation" of its main axis with slope similar to the critical state. Other models do not show this "rotation." This rotation is attributed

Simulation of undrained triaxial creep tests
Undrained creep has been observed in many experimental works with soft soils, 100-104 and its simulation has been addressed in other works. 38,39,44,50,59,105 The usual procedure is to perform an undrained triaxial test until a small value of axial strain is reached. Thereafter, the axial stress 1 and confining pressure 3 are kept constant. In terms of Roscoe invariants q and p, the deviator stress is kept constant . q = 0, while the mean (effective) stress p is allowed for variations. These conditions result in a creep process leading to the development of pore pressures p w and axial strains 1 under undrained conditions ( . v = 0). Details of undrained creep can be found in other works. [106][107][108][109] Under the mentioned conditions, an acceleration of the creep process is sometimes evident when the deviator stress is close to the maximum deviator stress q f , ie, q ≈ q f . This creep acceleration has been refereed to as "tertiary creep" by some authors. 107,110,111 Eventually, a pronounced acceleration of the axial strains may lead to the failure of the sample and loss of controllability. This case is known as "undrained creep rupture" and has been also studied in other works. In the following lines, we examine the capabilities of the different models with respect to these observations. Undrained creep is simulated with the four models. Initial conditions under isotropic states q = 0 were assigned to match a void ratio e = 0.739 and mean stress of p = 100 kPa. These initial conditions coincide with a point at the reference isotach, ie, R ref = 1. An undrained shearing under triaxial conditions is then performed with a strain velocity of || . || = 10D r . The undrained shearing is stopped at different axial strain levels 1 . Subsequently, undrained creep (q = 0) is simulated for a period of 10 6 seconds. Figures 12 and 13 present the simulations in the q vs 1 space and p vs q space. A monotonic undrained test with the same strain velocity has been also included for comparison purposes. Four different creep paths were simulated, each at different q. Simulations show that all models are able to simulate undrained creep, which is obvious because of their formulations. The resulting strain velocity under the undrained creep is examined in Figure 14. The only models evidencing simulations of tertiary creep, ie, a strain acceleration, correspond to the NSFS and VHP models; see Figure 14A,D. Figure 14D shows that these special cases ended with creep rupture because the condition . s → ∞ is reached. Other models did not show evidence of this capability. According to Vermeer and Neher, 59 a NVP-type model can also reproduce tertiary creep, when the viscous strain is formulated as follows: Note that is a function of the stress and tends to → ∞ for states very close to critical states. Some simulation examples with this modification are shown in Figure 13 in p−q and q− 1 plane, and in Figure 14(E), the strain acceleration leading to creep rupture is presented.
We now examine the creep rupture condition under undrained triaxial tests for the models NSFS and VHP. Consider triaxial conditions, such that one may simplify the constitutive Equations (46) and (28) to where . p = −( Computation of the scalar function K 22 is not an easy task. It depends on the mean stress p, the deviator stress q, and the void ratio e. For the case of the NSFS model, it also depends on the state variable vp v , which makes its solution path dependent. Hence, we find a solution only for the VHP model. For a given void ratio e, one may numerically solve a set of points p − q solving the condition K 22 = 0. For this purpose, a special script in the software WOLFRAM MATHEMAT-ICA has been implemented. Some details of the programming lines are given in Appendix C. Figure 14(F) shows some envelopes of the creep rupture condition K 22 = 0 for the VHP model. Three different solutions are found for different void ratios e = {0.7, 0.75, 0.8}. The surface described by the condition OCR 3D = 1, see Equation 8, is also included for comparison purposes. Notice that for a given void ratio, the rupture curve is found to the right side of the CSL line at high deviator stresses q. This is in accordance with simulations and also experimental observations. Now, we will compare the simulated response of undrained triaxial creep tests (UTCTs) of Haney clay with test data provided by Vaid and Campanella. 112 The material parameters are listed in Figure 2 Figure 15, it is also evident that all models, except for OVP, reach a vertical tangent. Hence, NSFS, NVP, and the VHP models are able to simulate tertiary creep in accordance with experimental results and the aforedescribed qualitative simulations. At this point, it is worthy to note that the flow rule term , introduced in Equation (85), could also be applied to an overstress-type model, making OVP models capable of capturing tertiary creep.
The accurate simulation of these effects are of great importance for the assessment of the failure probability of a creeping slope due to accelerated (tertiary) creep. Creep in general is indispensable for predicting the flow pressure on dowelling elements of these slopes.

Simulation of an oedometric test with different strain rates
The assessment of the models is now examined with simulations of an experimental study. The tested material corresponds to a Kaolin clay reported by Niemunis and Krieg. 7 The Kaolin is a medium plasticity clay with plasticity index equal to PI = 12, 2%. Reconstituted samples of height and diameter equal to 50 mm were used for this purpose. The experiment consists of a volumetric compression under oedometric conditions considering strain rate effects and unloading-reloading    Table 2, and the characteristic void ratios e i0 describing the OCR= 1 isotach in Table 3. In Figure 16, dashed lines represent the experimental curve, while continuous lines represent the simulations.
The simulations showed some relevant aspects of the model capabilities: All models, except by the NSFS model, were able to reproduce creep after an unloading path. The sequence loading-unloading-creep forced the NSFS model to remain under elastic response and thus to present at that instance a complete inability to simulate creep. This implies that, for the NSFS, an infinitesimal unloading disables the creep simulation. This is of course against experimental observations and may be considered as one of the main disadvantages of classical NSFS models. Other strain-rate dependencies were satisfactory captured by all models.

DISCUSSION AND FINAL REMARKS
The present work was devoted to review four different models for viscous clays belonging to different model families, namely, NSFS, NVP, OVP, and VHP. In order to perform a fair comparison, the developed formulations accounted for the same material parameter definitions, same reference isotach, and other rules established in Section 5. Numerical investigations allowed to detect a number of discrepancies related to their formulations. Some of them are in the following lines mentioned, while others were directly within the text indicated: Discrepancies observed at simulations of creep and stress relaxation are mostly related to the consideration of a maximum and/or minimum isotach in their formulations. While the OVP model incorporates a minimum isotach, at which no more rate effects are experienced for creep and stress relaxation, the NSFS and VHP model present a maximum isotach, which bounds instant time effects. The NVP model does not incorporate any of these bounding isotachs. Response envelopes showed that viscous effects are exhibited only at certain stress rate directions by all models, and a loss of convexity in conjunction with the appearance of corners were only noticed on those from the NSFS model. The undrained creep analysis showed that models having a strain rate with three components, such as the NSFS and VHP, are prone to include the simulation of undrained creep rupture in their capabilities. The simulation of a real oedometer test, including strain-rate dependencies, showed that the NSFS model can not simulate creep after an unloading path. Finally, it seems that the consideration of a maximum isotach allows the formulation to enable the simulation of viscous and nonviscous clays. The latter point was not illustrated through simulations but shown directly through their formulations. Table 4 provides a summary of these observations. As a final note, it is recalled that the observations listed in Table 4 hold only for the formulations inspected herein, which were proposed according to their original works for each model family. The consideration of a Buismann-type relation for the description of creep processes leads to an accurate estimation of the viscous behavior such as for example the prediction of creep settlements after surcharge loading. Creep in general is indispensable for predicting the flow pressure on dowelling elements of creeping slopes. All models can be used for this purpose. The accurate simulation of tertiary creep is of great importance for the assessment of the failure probability of these slopes due to accelerated (tertiary) creep. The NSFS, VHP, and NVP (with the extension proposed in Vermeer and Neher 59 ) models have shown the ability to describe these effects well. Furthermore, it can be concluded that the simulation of tertiary creep can be achieved by different assumptions ie, using three components models (elastic+plastic+viscoplastic) or a modified flow rule (with → ∞ for states close to critical states).
Of course, special extensions of the studied models may provide different features and therefore different conclusions but were not studied to keep simplicity.
• The characteristic void ratio e ref0 is found by adjusting the reference isotach, described by e = e ref0 − ln(p), to the compression curve with . v = D r . • The viscosity index I v can be found through two isotachs under isotropic compression (q = 0). To that end, consider the overconsolidation ratios of OCR a and OCR b and volumetric strain rates of . v,a and . v,b for the two isotachs, respectively. One may use the relation for the NVP model (see Equation 71): ) . (A1) • The poisson ratio can be adjusted to match the stiffness of an undrained triaxial test. At the q vs. 1 space, the initial stiffness of an undrained triaxial test is equal to 3G, where G is related to through Equations 24. A simple trial and error procedure may be also useful to find a which renders accurate simulations for undrained tests.

APPENDIX B: REQUIRED DERIVATIVES
In the following lines, the required derivatives for the semi-implicit integration of the viscoplastic strain rate . vp for the OS and VP models and viscous strain rate . vis for the VHP model are given.

B.1 OVP and NVP models
For these models, the derivative ( . vp )∕( ) is found as follows: with Note that the third stress invariant has not been considered.

B.2 VHP model
For the VHP model, the derivative ( . vp )∕( ) is found as follows: with . vis The derivative ( . vis )∕( . ) is found as follows: . with

APPENDIX C: PROGRAMMING THE CREEP RUPTURE CONDITION FOR THE VHP MODEL
For triaxial conditions, we may simplify the VHP model to the following matricial form: whereby

33
For the same conditions, the elastic stiffness matrix [E] may be written as follows: [E] = K [ 1 1 1 and the stress ratio is computed as The hypoplastic flow rule is computed as

APPENDIX D: MATHEMATICA SCRIPT FOR THE CALIBRATION OF E I0
In order to compare the VHP and UH models with the VP model, we calibrated them to have the same reference isotach.