Finite Element Simulations of an Elasto-Viscoplastic Model for Clay

In this paper, we develop an elasto-viscoplastic (EVP) model for clay using the non-associated flow rule. This is accomplished by using a modified form of the Perzyna's overstressed EVP theory, the critical state soil mechanics, and the multi-surface theory. The new model includes six parameters, five of which are identical to those in the critical state soil mechanics model. The other parameter is the generalized nonlinear secondary compression index. The EVP model was implemented in a nonlinear coupled consolidated code using a finite-element numerical algorithm (AFENA). We then tested the model for different clays, such as the Osaka clay, the San Francisco Bay Mud clay, the Kaolin clay, and the Hong Kong Marine Deposit clay. The numerical results show good agreement with the experimental data.


Introduction
Saito and Uezawa Saito and Uezawa [1] showed that the Takabayama landslide in Japan was mainly due to the time dependent deformation of the clay, resulting in the failure of the slope. In 1963, the failure of the Vayont reservoir in Italy, cost more than 3000 lives [2]; this was again due primarily to the creep of clay. These and many other examples in geotechnical engineering, reservoir geomechanics, wellbore plugging in petroleum and mining industry, etc., have made studying and modeling of clay a very important topic.
Clay is considered to be a multi-component system composed of solid particles, and fluid components (gas and/or liquid) [3]. For many contaminated soils, there might be additional components, such as bubbles, oil spill, etc. In a partially saturated clay, the void space is filled with gas and/or liquid, and, in a fully saturated clay, only a liquid phase is present. Removal of any fluid from this system depends on the magnitude of the external load and the drainage path, the imposed duration, and the properties of the clay. When an external load is applied, pressure is developed in the fluid, which may dissipate instantaneously or gradually [4]. On the other hand, for a solid medium with non-crushable fine materials, when the fluid pressure dissipation is negligible, the settlement of clay does not end. In this case, creep may continue for a long time under a constant pressure [5]. This process constitutes a complex hydro-mechanical phenomenon [6,7].
Many constitutive models have been developed for clay. These range from simple elastic models to very complicated elasto-viscoplastic models. An overview can be found in the works of Owen and Hinton [8] and Desai and Siriwardane [9]. Chaboche [10] also presented a review of constitutive theories for plastic and viscoplastic type models (see also Tatsuoka et al. [11]). Liingaard et al. [12] suggested that the time-dependent viscous type models can be separated into: (i) empirical models (Singh and Mitchell [13]); (ii) rheological models (Feda [14]); and (iii) models considering "stress-strain-time" (Adachi and Okano [15]). The first two categories are important for the basic understanding of the behavior of clay. In clay-based research and its finite element implementation, the third category is most popular. For the remainder of this paper, the discussion is limited to the elasto-viscoplastic (EVP) type models. In the formulation of the EVP models, the strain rate consists of an elastic component and an inelastic component. To obtain an expression for the latter, it is essential to include the time-dependent viscous property in the model formulation. Sekiguchi [16] subdivided EVP models into: (a) the overstressed type [15]; (b) the non-stationary flow surface type (Sekiguchi [17]); and (c) others, such as the Bounding surface model [18] and Borja model [19]. Details of EVP models and their strength, as well as their limitations, can be found in the work by Liingaard et al. [12].
In this paper, we develop a new EVP model with a non-associated flow rule, considering Perzyna's viscoplastic theory [20] and the MCC framework [21]. The time-dependent viscous aspect is incorporated using the Borja and Kavazanjian [19] approach. The new model requires six parameters. The current study is a generalization of the previous work by Islam and Gnanendran [22], who showed the strengths and the limitations of the existing methods when using the twosurface approach in the EVP models by incorporating the associated flow rule. Here, we use the nonassociated flow rule and the three-surface approach. We also look at the Osaka clay, the San Francisco Bay Mud clay, the Kaolin clay, and the Hong Kong Marine Deposit clay in a triaxial loading situation. The wide range of experiments included the drained and the undrained tests, the strain rate test, the relaxation test, and the over consolidation effect tests. The model was implemented in a finite-element code using the numerical algorithm (AFENA) [23]. The pertinent details of the finite element implementation and the development of the non-associated flow rule are discussed in this paper.

Basic Equations
We consider clay to be a mixture of two constituents composed of solid particles and a liquid. We use the basic ideas and principles of continuum theories of mixtures, where electromagnetic and thermo-chemical effects are ignored (see Atkin and Craine [24] and Bowen [25]). The motion for each component is given through a mapping where designates the mapping, and refers to (solid component-clay) and (fluid component). For clay, the kinematical quantities of interests are the displacement vector , the velocity vector , the deformation gradient tensor , and the linearized strain : where superscript designates the transpose of a tensor. If is the volume of the particles before the deformation and is the volume after the deformation, then a quantity of interest is the volumetric strain, also called the dilatation Furthermore, in elastoplastic theories, for small deformations, we use ̇ instead of the symmetric part of the velocity gradient . That is, for small deformation, we assume ̇= , where For the fluid component, the kinematical quantities of interest are the velocity vector , the gradient of the velocity , and its symmetric part , which are given by, respectively, We assume that the density of the solid and fluid components in their current configuration are and , respectively, while, in their pure (reference) state, i.e., before mixing, they are given as and . The density and the velocity of the mixture are then defined as Note that Equation (12) is one of the many possibilities for defining a mixture velocity. We can also define a relative velocity such that The densities in the current and reference configuration are related through the kinematical field ( , ), called the porosity, such that where 0 ≤ ( , ) ≤ < 1. Thus, when = 1, we have a fluid with no pores and no particles. Note that this automatically assumes that the mixture is saturated. Otherwise, the porosities are constrained by + ≤ 1. In geomechanics-related problems, we sometimes use the void ratio , which is related to the porosity , through Assuming no inter-conversion of mass between clay and the liquid in the pores, the equations for the balance of mass for the two constituents are given by where (.) is the time derivative and is the divergence operator. The balance of linear momentum for the two components are given by where is the unit vector acting on the boundary. In accordance with the basic principles of mixture theory, advocated by Truesdell [26], if we add the two equations for the conservation of mass or the conservation of the linear momentum, we obtain the balance of mass and the linear momentum for the mixture (note that, in this case, drops out of the equations), due to application of the Newton's third law Finally, the balance of angular momentum indicates that, in the absence of angular momentum supply, the total stress is symmetric. That is [see Rajagopal [27] Since we consider a non-isothermal problem, we do not discuss the balance of energy or the Second Law of thermodynamics.

Constitutive Modeling
In this paper, the emphasis is on the modeling of the stress tensor of the solid particles (clay). However, for the closure of the governing equations, we also need constitutive relations for the stress tensor for the fluid, , and the interaction forces .

Fluid Component and the Interaction Forces
The following assumptions, as elaborated in Martins-Costa et al. [28] and Rajagopal [27], would lead to the classical Darcy's equation. The interaction forces designated by , within the context of mixture theory and many of the multiphase theories, are usually based on generalizing the interaction force for very special cases, such as the Stokes drag on a single spherical particle.
We assume that the frictional (viscous) forces within the fluid can be ignored and as a result the partial stress tensor for the fluid can be given by a Eulerian fluid model: where is, in general, a function of density, and is the identity tensor [note that compressive stresses are assumed to be negative in these theories, whereas in, geomechanics-related problems, the opposite convention is used]. We further assume, as is customary in geomechanics problems and basic flows through porous structures (see Oka and Kimoto [29], p. 34): The interaction force is assumed to be where is a coefficient that can depend on porosity, viscosity, permeability, etc. This is basically a generalization of the Stokes' drag on a single spherical particle. In general, the interaction forces could depend on other kinematical quantities such as the relative acceleration, velocity gradients, etc. (see Massoudi [30,31] for a review of this topic). To obtain the Darcy's equation, we ignore the inertial effects, i.e., we ignore the left-hand side of Equation (18) (when written for the fluid component), and, by using Equations (13) and (23), we obtain the Darcy's equation: where , as mentioned before, is the density of the fluid. Furthermore, by assuming (see Martins-Costa, et al. [28], Williams [32]) where represents the specific permeability, which for anisotropic materials is generally a secondorder tensor. Equation (26) can be re-written as In soil mechanics literature, Darcy's equation is sometimes expressed using the concept of hydraulic conductivity ( ) ( see Bear and Bachmat (1990, p.294)), defined by = (29) Here, represents the volumetric weight of the fluid, , which is assumed to act in the vertical direction ( = 0, , 0 ). Equation (28) is then re-written as: We use this form of the equation in our finite element simulation.

Solid Component
The partial stress for the solid component can be defined as (Lewis and Schrefler [33]): where is the effective stress tensor and is the pore (fluid) pressure. We can relate the total stress tensor of the mixture [see Equation (19)] to the effective stress tensor [by adding Equations (19), (24) and (31)], namely We assume that the strain in clay can be decomposed into an elastic part and a viscoplastic part. In plasticity theory, we find it more convenient to assume this decomposition applies to for the strain rates [see Davis and Selvadurai [34], p.97] We also assume that the elastic part of the strain can be represented by the "small-strain" or the linearized theory of elasticity, where, as customary in soil mechanics, the strain is assumed to depend on the effective stress (Terzaghi [4], pp.11-15), Schofield and Wroth ( [35], p. 9). For an isotropic material, using the index notation, the elastic strain is given by (Matsuoka and Sun [36],p. 37) Where in accordance with the critical state theory, we assume = ( )( ) is the (modified form of the) Young's modulus, is the Poisson's ratio, is the slope of the unloading and reloading path (see [Matsuoka and Sun [36] (p.35, Figure 2.8); Desai and Sriwardane [9](p. 289, Figure  11.7)), = is the initial mean pressure, is the initial void ratio, and is the Kronecker delta ( = 1( = ), 0( ≠ )). In a more compact form (using the index notation), Equation (34) can be written as where is the fourth-order compliance tensor, related to the stiffness tensor. Since we have assumed that the material is isotropic, in short hand notation [recalling Hooke's law ̇ = ̇ ], In Equation (35) ̇ can be generalized to the case of an elasto-viscoplastic case (see Desai and Sriwardane [9], p. 294, Equation 11.32), i.e., where is a fourth-order tensor similar to the elastic moduli. For the viscoplastic modeling of the strain rate, we start with Perzyna, who used the associated flow rule and assumed the plastic potential function coincides with the loading function [20]. To apply this theory to geomaterials, the main challenge is to define the static loading function [12]. By definition, , represents the stress state where the strain rate is assumed to be zero. Here, we assume that the viscoplastic part of the strain rate in Equation (33) is based on Perzyna's approach, where In the above equations, is the rate sensitivity function, and its functional form can be obtained either experimentally or theoretically; ⟨ ⟩ is the Macaulay's bracket (in Equation (38), the Macaulay's bracket ensures that the function inside the bracket only has a value when it is positive, otherwise its magnitude is zero); is the over-stress function; and is a new term, representing the new potential surface (surfaces of the proposed EVP model are obtained by extending the Modified Cam Clay surface [Roscoe and Burland [21]; in our EVP model, we require a total of three surfaces (see Figure 1): the loading surface , the reference surface , and the potential surface ) given by Potential surface: and and are given by the same expression as those in the Perzyna model. They are the dynamic loading function (the potential surface) and the static loading functions, respectively, given by Reference surface: Loading surface: where = = , = √ = ( ) ( ) . Similar expressions for and (also known as the deviatoric stress) can also be found in the work by Borja and Kavazanjian [19]. In the principal stress space, and are the octahedral normal stress and the octahedral shear stress, respectively (see Matsuoka and Sun [36], pp. 29-30). In the above equations, the suffixes , , and represent the reference surface, the loading surface, and the potential surface, respectively. The meaning of these surfaces is shown in Figure 1. At any given time, the reference stress state and the reference surface are known from the laboratory test. The current stress state and the potential stress state are related to the reference stress state through the radial mapping rule, which was proposed by Phillips and Sierakowski [37]. In this paper, we used two image parameters and the details are given in Appendix A. It is worth mentioning that Islam and Gnanendran [22] demonstrated the strengths and the limitations of the existing methods using the two-surface approach in the EVP models. They used the associated flow rule for their EVP model, while, in the present paper, we use the non-associated flow rule and the three-surface approach. In the previous paper, the surface shapes are two ellipses, while, in the present study, we assume all the surfaces are given by a single ellipse, which is close to the MCC surface. To formulate the new EVP model, we assume that the "projection center" is in the origin of the stress space, which is identical to the MCC model. However, to define the surfaces, the expression of the slope of the critical state line ( ) is changed with respect is the internal angle of friction at the failure for each -value test, ranging from 0 (triaxial compression) to 1 (triaxial extension).
We have introduced a new parameter in order to obtain a more realistic surface shape in the -plane (see Islam and Gnanendran [22]). It is observed that, in any stress state (0 < − ≤ 1) other than the triaxial compression state ( − = 0), the MCC equivalent surface overestimates the stress. For the sake of completeness, in Figure 2, we compare the EVP model based on the MCC surface with the new modified surface presented in this paper. It is observed that the newly extended MCC surface captures and compares well with the experimental results.
In Figure 3, we can see that, at any arbitrary reference time ̅ , the soil state is at "A", where the corresponding void ratio is ̅ .With time changing from ̅ to t, due to creep, the soil moves from "A" to "B" where the corresponding void ratio is . Then, the following expressions are obtained from Borja and Kavazanjian [19]: where and are the slope of the normal consolidation line (the -line) and the unloadingreloading line (the -line), respectively, and is the void ratio corresponding to the -line when = 1 kPa at ̅ . It should be mentioned that , and are the necessary parameters in the EVP model; these can be obtained either from the oedometer test or the triaxial test. The meaning of these parameters is the same as those in the MCC model. As time increases, the -line changes to the ̇line and is transformed to ̇ . The ̇line will generate the new bounding surface. The -line and the ̇line are parallel, as shown by Bjerrum theory [39]. Using Borja and Kavazanjian's [19] concept and the multisurface theory, ̅ is the arbitrary time, representing the state of stress prior to the surface evolving. In Equations (40) and (41), and are also known as the creep exclusive preconsolidation pressure and the creep inclusive preconsolidation pressure, respectively (see Islam and Gnanendran [22]). The expression for is similar to the one used in the MCC model: After rearranging Equation (44), the expression for can be obtained From Figure 3 for the definitions of and , we can obtain an expression for : The detailed derivation of ⟨ ( )⟩ and ̇ are presented in Appendices B and C, respectively. In closing this section, we need to mention that the overstressed EVP models are usually based on the Perzyna's theory [20] in combination with the critical state soil mechanics theory, e.g., the Modified Cam Clay (MCC) model [21]. In these approaches, the viscous nature of the EVP model is introduced in the theory using a secondary compression index, a creep function and a relaxation function. However, in most cases, to reduce the complexity of the model and to minimize the number of parameters, the EVP models are developed considering the associated flow rule, where it is assumed that the yield surface is identical to the potential surface. To capture the behavior of geomaterials, using the non-associated flow rule is essential [40]. Depending on the application of the Critical States Soil Model (CSSM) in the EVP models using the non-associated flow rule, there are two approaches we can consider: (i) those with critical state [40]; and (ii) those without critical state [41]. The required parameters for the EVP models, satisfying the non-associated flow rule, ranges from 7 parameters [40] to 44 parameters [41]. A summary of the EVP models with the non-associated flow rule for different geomaterials is presented in Table 1.

Model Parameters
There are six parameters in our model that need to be determined: the consolidation parameters can also be found in many soil mechanics textbooks.

Summary of the Basic Equations and the Assumptions Used in the Code
The equations used in the code and the finite element formulation do not have the same exact forms as those given in Sections 3.1-3.2. In this subsection, we present a summary of the derivation of the equations used in the code and we show how these equations can be obtained from the equations in the earlier sections, if proper assumptions are made. For example, the momentum (equilibrium) equation for the mixture, as used in finite element code, can be obtained if we write Equation (18) for the two components, add them up, assume that the inertial terms can be ignored, use the definition of the total stress tensor [Equation (19)], and the definition for the total density [Equation (11)], then we have where is the acceleration due to the gravity. If we substitute Equation (32) into Equation (48), and use the convention that a compressive pressure is positive, we have Similarly, if the conservation of mass equations, given by Equation (17), are re-written for the two components, for the cases of constant densities, i.e., when = constant and = constant, they become Adding these two equations, we obtain In this equation, the term = ( − ) is known as the "specific discharge" and is related to the relative velocity . We also notice that [recall that = was defined as the dilatation (see Equation (6) ]. Now, using Equations (53) and (30) in Equation (52), we obtain This expression is also known as the "storage equation" which is the basic equation for the consolidation theory ( see Bear and Bachmat [48], Equation 4.1.46).

Finite Element Solutions
The set of equations that need to be solved is: (i) the mass balance equation; (ii) the equilibrium equation; (iii) stress-strain relations; and (iv) strain-displacement relations. The unknown variables are the stresses, the strains, the displacements and the pore pressure. A summary of the equations and the unknown variables is presented in Table 2. Using the non-associated flow rule EVP model presented in this paper, we can see that there are enough equations to solve for the unknown variables. For any deformable porous media, similar expressions can be found in the work of Bear and Bachmat ( [48], Chapter 4).
To solve these equations using finite element, we can use either "the strong form" or "the weak form". The details of these two techniques can be found in many finite element text books (Zienkiewicz et al. [49]). Even though both approaches provide similar results, considering the "continuity requirements" and the "symmetry of the stiffness matrix", the weak form solution is better compared to the strong form [49]. In this paper, we use the weak form approach. To obtain the weak form, the Galerkin weighted residual method [50] is introduced. Equations (48) and (54) are used to solve the problem. In the following section, a summary of the weak formulation is described.

Weak Form Formulation
Equilibrium equation Darcy's equation Continuity equation Effective stress relation In Equations (55)- (58), and ′are the total stress and the effective stress, respectively, is the body force, is the differential operator, ∇ is the divergence operator, ∇ is the gradient operator, is the pore pressure, is the superficial velocity, is the coefficient of permeability, is the unit weight of the fluid, is the body force vector for the fluid, and ̇ is the volumetric strain rate of soil. The general expression of , ∇, , , and in the above Equations are  Assuming that the soil is subjected to the traction , the following conditions can be written. Initial conditions: Assuming that soil is isotropic, the stress tensor , the stain tensor and the displacement vector can be expresses as To obtain the weak formulation for the equilibrium equation, basic necessary steps are: (i) replacing Equation (58) in Equation (55) and multiplying the resulting equation with an arbitrary function, which removes the essential boundary condition or the Dirichlet boundary condition [51]; (ii) integrating the system over the domain; and (iii) applying the Green-Gauss theorem [50] and integrating by parts .
Now, ′̇ (see Equation (37)) can be written as where, in Equation (37)  In the time interval Δ = − , the viscoplastic strain rate increment Δ can be defined as [7,8,52] where = 0, denotes the "fully explicit" Euler time integration scheme (or the forward difference method) and = 1 indicates the "fully implicit" Euler time integration scheme (or the backward difference method). For ≥ 0.5, the integration scheme is unconditionally stable and, when < 0.5 , the integration process is conditionally stable. The integration scheme for = 0.5 is the "implicit trapezoidal" scheme or the Crank-Nicolson rule.
The viscoplastic strain rate increment ̇ at time step (n +1) in Equation (106) can also be written as follows [8] where Δ is the stress change at any time interval Δ = − and (given in Appendix D) is the derivative of the viscoplastic strain-rate with respect to time at the nth time step.

Model Verification and Discussion
We investigated the performance of our EVP model in a variety of experimental applications. These applications included the consolidated undrained triaxial compression test and the extension test, the consolidated drained compression test, the over-consolidation ratio test, the confining pressure effect, the strain rate test, the creep test, and the relaxation test. In this respect, we considered both natural clay and laboratory obtained reconstituted clay. We specifically investigated the Osaka clay [53], the San Francisco Bay Mud (SFBM) clay [54], the Kaolin clay [55], and the Hong Kong Marine Deposit (HKMD) clay [56]. Clay properties are presented in Table 3 and the triaxial loading configuration is shown in Figure 4.

Simulation of the Undrained Triaxial Test on the Osaka Clay
Adachi et al. [53] presented undrained triaxial compression test for the natural Osaka clay. This is a moderately sensitive clay (sensitivity, = 14.50). The water content, the liquid limit and the plastic limit of this clay are 65.00-72.00%, 69.20-75.10% and 24.50-27.30%, respectively. Its specific gravity is 2.62-2.70. The particle size distribution of this clay consists of the clay fraction of 44.00%, the silt fraction of 49.00% and the sand fraction of 7.00%. The strain rate during the test was 1 × 10 %/minute. From the experimental results, it is evident that the Osaka clay shows strain softening behavior, which was captured using the new EVP model presented in this paper. In Figure 5, the experimental results for the Osaka clay for two different mean pressures (176.4 kPa and 235.2 kPa) are compared with the non-associated flow rule type EVP model, and the EVP model developed by Kutter and Sathialingam [42]. In Figure 5-a, we can see that, for the mean pressure of 235.2 kPa, the new EVP model and the Kutter and Sathialingam [42] EVP model predictions were identical up to an axial strain of 1.87%. After this point, their model over-predicted the observed experimental results. After an axial strain of 14%, the over-prediction in their model was 16%, whereas the new EVP model captured the strain softening. On the other hand, at an axial Deviatoric stress,q (kPa) strain of 14%, the under-prediction in the new EVP model was 5.20%. For small strain cases (axial strains of 1-3.75%), an over-prediction was observed in both EVP models. Jiang et al. [57] reported that, at small strains, the over-prediction of the nonlinear responses of clay are expected. Using the hysteretic response equation [58], such issues can be resolved; however, this requires additional model parameters. To make the present EVP model formulation simple, and to limit the number of model parameters, the hysteretic response equation concept was not introduced here. A similar trend was also observed for the mean pressure of 176.4 kPa.
In Figure 5-b, the stress paths obtained from the associated flow rule, the non-associated flow rule EVP models, and the Kutter and Sathialingam model are compared with the data for the Osaka clay. It is evident that, after the deviatoric stress reached the peak, the attributed stress paths gradually followed the "narrow region". This phenomenon indicates that the critical state concept is applicable to the natural soft clay, even at large strains [53]. It was observed that the new nonassociated flow rule EVP model could capture the "narrow region", whereas there was marginal under-prediction in the associated flow rule EVP model.

Simulation of the Undrained Triaxial Stress Relaxation Test on the San Francisco Bay Mud Clay
Lacerda [54] presented the results of undrained triaxial stress relaxation tests for the undisturbed San Francisco Bay Mud clay. In this paper, only SR-I-5 test data are presented and compared with the EVP model predictions. The sample was isotropically consolidated to a pressure of 78.4 kPa. Basic properties of this clay include the specific gravity = 2.66-2.75; the liquid limit = 88.4-90%; the plastic limit = 35-44 %; the plasticity index = 45-55%; and the moisture content = 88-93%. Kaliakin and Dafalias [59] reported that the initial void ratio ( ) and its corresponding mean pressure ( ) for the San Francisco Bay Mud are 1.30 and 156.9 kPa, respectively. The details of the undrained triaxial stress relaxation tests on the SFBM clay are presented in Table 4.  Figure 6, the measured data for the San Francisco Bay Mud clay are compared with the numerical simulations for the associated flow rule and the non-associated flow rule EVP models. In addition, our results are also compared with the EVP models of Kaliakin and Dafalias [59] and Kutter and Sathialingam [42]. It is evident that the non-associated flow rule prediction was close to the laboratory results.

Simulation of the Consolidated Undrained Triaxial Tests on the Kaolin Clay
Herrmann et al. [55] presented the effect of over-consolidation ratio (OCR) on the Kaolin clay for the undrained triaxial compression test and the extension test. This clay is a mixture of the Kaolin (Snow-Cal 50) and a 5.0% Bentonite by weight, which was prepared in the laboratory. The Kaolin clay properties are: the specific gravity = 2.64, the liquid limit = 47.0%, and the plastic limit = 20.0%. The sample was isotropically consolidated to a pressure of 392.2 kPa; to achieve OCR = 1, 2, 4 and 6, the confining pressure was changed accordingly. For the OCR = 1, the corresponding initial void ratio ( ) is 0.613. The deviatoric stress versus the axial strain predictions for the consolidated undrained triaxial compression test (OCR = 1, 2, 4 and 6) and the extension test (OCR = 1 and 2) were compared with the experimental results, as shown in Figure 7. For the normally consolidated clay (OCR=1), it was observed that the new EVP model captured the stress-strain responses well before the deviatoric stress reached its peak, but then slightly under-predicted them (at 14% strain, the under-prediction was only 1.05 %). In Figure 7a, it is evident that, when OCR = 2, 4 and 6, before reaching the peak deviatoric stress, the model over-predicted but then exhibited small amount of under-prediction. In the consolidated undrained extension tests, a similar trend was observed. When OCR = 1 and 2, we compared the experimental results with the predicted responses of the pore-water pressure for the consolidated undrained triaxial compression test and the consolidated undrained triaxial extension tests, as shown in Figure 7b. For the normally consolidated clay, predictions for both the compression test and the extension test were satisfactory. For the OCR = 2, in the triaxial compression test, the new model slightly over-predicted before the maximum pore pressure was reached and, in the triaxial extension, the negative pore pressure was well captured with a small degree of under-prediction. A similar pattern was observed for the stress path, as shown in Figure 7c.   Compression

Simulation of the Consolidated Drained Tests on the Hong Kong Marine Deposit Clay
Kaolin clay  The Hong Kong Marine Deposit clay [see Yin and his co-workers [56]] was used to study the isotropically consolidated drained test. This is a reconstituted soft to very soft illitic silty clay, which has been extensively studied.
The properties of the Hong Kong Marine Deposit (HKMD) clay are: the specific gravity = 2.66, the liquid limit = 60%, the plastic limit = 28%, the plasticity index = 32%, and the moisture content = 51.7%. The predictions of the new model for the normally consolidated drained triaxial test on the HKMD clay are shown in Figure 8a for the mean pressures of 300 kPa and 400 kPa. It is evident that the present model captured the deviatoric stress versus axial strain response very well (Figure 8). The volumetric response and the stress paths are presented in Figure 8b, c, which were also satisfactory.

Perturbation Analysis
We performed a perturbed analysis for each model parameter for the Osaka clay and the Kaolin clay. The first one is a natural clay, while the second one is a reconstituted clay. The objectives of such analyses were to quantify the sensitivity of each perturbed parameter to the calibrated reference parameters. The mean value of the error ( ) was calculated for each parameter as  where N is the number of error calculations at different strain values; and and are the predicted values of the deviatoric stress at strain using a reference parameter and a perturbed model parameter, respectively. The material parameters of the EVP models were divided into two groups: the Modified Cam Clay (MCC) parameters and the viscous parameters. Here, the parameters were perturbed by the same magnitude so that the influence of the perturbation of the parameter on model's prediction could be observed.
The parameters in the Osaka clay model were perturbed by 3% of their reference values, which were from experimental data. Then, for each parameter, (%) due to the perturbation was calculated for the associated flow rule and the non-associated flow rule EVP models, as shown in Equation (111). The perturbed results for the Osaka clay are shown in Figure 9. It was observed that parameters and were most sensitive, while M had medium sensitivity and and were less sensitive. To compare the perturbation influence on the flow rule, identical perturbation was also conducted for the Kaolin clay, as shown in Figure 10. In the perturbation analyses with the same magnitude, the effect of the flow rule was negligible.

Concluding Remarks
To describe the time-dependent behavior of clay, a new EVP model is formulated in this paper. The model requires six independent parameters, which are defined in the literature and can be extracted from the simple oedometer tests and the triaxial tests. A generalized nonlinear creep function is also used. The model captured a wide range of results when compared with experimental observations obtained from the drained and the undrained conditions for the triaxial compression test and the extension test, the creep test, the relaxation test, the over consolidation ratio effect test, and the strain rate test. Both undisturbed clay and reconstituted clay were considered in our study. From the comparison of the model predictions, it is evident that reasonably good agreement was obtained.

Appendix B. Derivation of the Rate Sensitivity Function ( )
Islam and Gnanendran [22] presented the derivation of the rate sensitivity function, ( ) for the associated flow rule and also discussed the limitations of these approaches. For the sake of completeness, we provide a summary of the derivation for the non-associated flow rule. From Islam and Gnanendran [22], the viscoplastic strain rate ̇ is given as If we decompose ̇ into the volumetric ̇ and the deviatoric ̇ parts, the expression of the viscoplastic strain rate ̇ in the conventional triaxial stress space can be shown to be Implementing Equation (A12) into Equation (A11), we obtain the complete form of the gradient matrix for two-dimensional case: