Non-Associated Flow Rule-Based Elasto-Viscoplastic Model for Clay

We develop a non-associated flow rule (NAFR) based elasto-viscoplastic (EVP) model for isotropic clays. For the model formulation, we introduce the critical state soil mechanics theory (CSSMT), the bounding surface theory and Perzyna's overstress theory. The NAFR based EVP model comprises three surfaces: the potential surface, the reference surface and the loading surface. Additionally, in the model formulation, assuming the potential surface and the reference surface are identical, we obtain the associated flow rule-based EVP model. Both EVP models require seven parameters and five of them are identical to the Modified Cam Clay model. The other two parameters are the surface shape parameter and the secondary compression index. Moreover, we introduce the shape parameter in the model formulation to control the surface shape and to account for the overconsolidation state of clay. Additionally, we incorporate the secondary compression index to introduce the viscosity of clay. Also, we validate the EVP model performances for the Shanghai clay,the San Francisco Bay Mud (SFBM) clay and the Kaolin clay. Furthermore, we use the EVP models to predict the long-term field monitoring measurement of the Nerang Broadbeach roadway embankment in Australia. From the comparison of model predictions, we find that the non-associated flow rule EVP model captures well a wide range of experimental results and field monitoring embankment data. Furthermore, we also observe that the natural clay exhibits the flow rule effect more compared to the reconstituted clay.


Introduction
In a saturated clay medium, the liquid phase occupies the interparticle void spaces of the solid phase. When such a clay deposit experiences external loading, saturated soil either exhibits slow loading conditions or fast loading conditions. The first one is also known as the drained behavior, while the latter one represents undrained behavior. In most cases, the liquid removal of soft clays is not instantaneous due to its low permeability, compressibility characteristics and the viscous nature [1,2]. Therefore, depending on the loading state and the in situ condition, in many cases, the deformation of clay may continue for a long time [3]. For example, in August 1173, the Leaning Tower of Pisa in Italy was constructed on a highly compressible clay deposit, and it displaced horizontally to the magnitude of 4.7 m in 1990, which is also increasing 1.5 mm/year [4]. A similar situation also may happen in any geotechnical structure founded on soft clay (see also Brand and Brenner [2]). In this regard, the viscosity of clay most often contributes to the long-term time-dependent creep of clay, and subsequent damage to the structure, which requires billions of dollars in annual maintenance costs [5].
Additionally, only in the USA, the sum of onshore abandoned wells is about 3 million, and every year approximately 40,000 new deep boreholes are drilled [6]. The estimated cost using the bentonite clay-based plug of those abandoned wells is about $160,000 per shallow borehole, which also demonstrates a billion-dollar legacy (see also Islam et al. [6]). Furthermore, in an engineered barrier system for hazardous wastes (e.g., nuclear waste), the application of bentonite clay is also common and the clay type barrier requires 100's of years of monitoring due to the sensitive nature of the deposited materials [7]. Additionally, it is important to incorporate time-dependent viscous responses of clay in a constitutive model formulation to obtain realistic hydro-mechanical behavior [8]. Therefore, clay-based research is the active field of interest and the motivation of the present paper.
From the early 19 th century to the present time, to illustrate the clay behavior, a myriad of coupled constitutive models have been developed, including viscous-inclusive (e.g., Adachi and Okano [9]) and viscous-exclusive (e.g., Roscoe and Burland [10]) models. However, we limited our discussion only to the first group of models. In this regard, among others, Liingaard et al. [11] and Chaboche [12] provided literature reviews. Nevertheless, to avoid the mathematical formulation complexity, in most cases, the time-dependent constitutive models are limited to the associated flow rule (AFR). But, capturing the legitimate behavior of soft clay, the non-associated flow rule (NAFR) is imperative (Zienkiewicz et al. [13]). In the literature, there are a couple of NAFR-based EVP models, where the number of material parameters ranged between six to 44. In this regard, Islam et al. [14] also presented a summary of NAFR-based elasto-viscoplastic (EVP) models. It is worth mentioning that a constitutive model with too many model parameters may capture geomaterials' behavior very well, but most often, their practical applications are not convenient [15]. Additionally, for the engineering application of any EVP model, the model formulation simplicity, its finite element implementation and the objective determination of model parameters are essential.
In this paper, we develop a non-associated flow rule-based elasto-viscoplastic (EVP) model considering the Modified Cam Clay (MCC) model [10] framework, Perzyna's overstressed theory [16], the Borja and Kavazanjian [17] concept, the bounding surface theory and the mapping rule (see also Hashiguchi [18]). The EVP model here requires a total of seven parameters, and among them, five parameters are identical to the MCC model. The other parameters are the secondary compression index and the surface shape parameter. We also introduce a non-linear secondary compression index to account for the viscosity of clay. Additionally, the shape of surfaces in this paper is different than the original MCC model or EVP model with the MCC model equivalent surface (see also Islam et al. [14]). We consider a non-circular surface in the π-plane. Also, we introduce a composite boundary surface, and the shape parameter to control the bounding surfaces' shape.
Furthermore, we also discuss the importance of the non-associated flow rule, the non-circular shape surface and the composite bounding surface. For validation of the EVP model, we compare numerical results with a wide variety of laboratory observed test data considering the Shanghai clay, the San Francisco Bay Mud clay and the Kaolin Clay. Moreover, after validation of the developed non-associated flow rule EVP model, we also implement it in a coupled finite element solver named A Finite Element Numerical Algorithm (AFENA) [19]. Additionally, for a field application of the developed EVP model, we compare the predicted response with the long-term monitoring measured response of the Nerang Broadbeach Roadway (NBR) embankment in Australia [20].

Importance of the Non-Associated Flow Rule
In this paper, we formulate the non-associated flow rule elasto-viscoplastic model considering the Modified Cam Clay (MCC) model framework, which was formulated considering the associated flow rule. Therefore, at first, we discuss the importance of the non-associated flow rule in the context of the MCC model as follows.
In the triaxial space, the incremental strain (݀ߝ) is divided into the volumetric (݀ߝ ௩ ) and the deviatoric (݀ߝ ) component, while each of them is also comprised of the elastic part and the inelastic part (e.g., plastic or viscoplastic) as follows (see also Yu [15]; Roscoe and Burland [10]): where, ݀ߝ ௩ and ݀ߝ are the total volumetric strain and the total deviatoric strain, respectively; while superscript ݁ and ‫‬ represent their corresponding elastic and plastic components. ݀ߝ ଵ and ݀ߝ ଷ are the incremental major strain and the incremental minor principal strain, respectively. Again, for the 1D consolidation by noting ߝ ଷ = 0 in Equation (3), the ratio of the incremental volumetric strain to the incremental deviatoric strain is obtained as: Additionally, by neglecting the relative magnitude of the elastic shear strain to the plastic shear strain, Yu [15] presented Equation (4) as: Also, Yu [15] presented the stress dilatancy relation as follows: In Equation (6), ݉, ݊ are the material constants. ‫ܯ‬ and ߟ́ are the critical state line (CSL) slope, and the stress ratio, respectively (see also Yu [15]). Furthermore, for the MCC model, ݉ = ݊ = 2. In Equation (5) (6), the stress ratio for the normal compression (ߟ́ , ) and the ‫ܭ‬ conditions become 0.40. However, ߟ́ , = ଷெ ିெ = 0.60 (see also Yu [15]). Additionally, McDowell and his coworkers also proposed that ߟ́ , = ‫ܯ6.0‬ = 0.60 (see McDowell and Hau [22]). Therefore, it implies that in the associated flow rule condition, the predicted ߟ́ , is too low.
Moreover, the yield surface of the MCC model is given by: where, ‫,݀‬ ‫ݍ݀‬ and ‫݀‬ are the differential quantities of ‫,‬ ‫ݍ‬ and ‫‬ respectively. It is to note that during the undrained triaxial tests, the magnitude of ‫݀‬ is a negative quantity for the normally consolidated clay. Therefore, in the 'wet' part when ߟ́< ‫ܯ‬ (see also Schofield and Wroth [21]): As a result, if the yield surface is not permitted to contract, as ‫݀‬ > 0 and ‫ݍ݀‬ > 0, the strainsoftening phenomena will not occur.
In addition, when the stress state reaches the yield condition, a material is subjected to the plastic deformation, which is known as the plastic flow (see also Hashiguchi [18]). The plastic potential function illustrates the post-yield and the failure behavior of soils, and the plastic strain increment is normal to it. If the plastic strain is calculated on any surface other than the potential surface as in the associated flow rule condition, the predicted plastic shear strain is too high (see also Yu [15]; McDowell and Hau [22]).
It is worth mentioning that in the associated flow rule condition, the adopted yield surface overestimates failure stresses on the dry side, and the bifurcation is not possible in the hardening regime (see also Yu [15]; Schofield and Wroth [21]). Additionally, in the associated flow rule (AFR), the yield surface and the potential surface are the same. From the evidence of the triaxial test results in the literature, it is observed that in the AFR condition, if the contraction of the yield surface with hardening is not allowed, the deviatoric strain is suppressed. Therefore, strain hardening behavior in the drained triaxial shearing entailed that the yield surface shrinkage with the hardening should not be permitted. Hence, in the associated flow rule condition, the strain hardening in the drained triaxial test and the strain softening in the undrained tests will not transpire. For this reason, to explain the plastic volume expansion and the softening behavior of plastically compressible materials, e.g., soft clay, a plastic potential surface rather than the yield surface is essential (see also Yu [15]).
To resolve the limitations mentioned above in the Modified Cam Clay (MCC) framework, we introduce the non-associated flow rule-based elasto-viscoplastic (EVP) model. Additionally, the MCC model is not capable of modeling the long-term viscous behavior of soft clay (see Islam and Gnanendran [23]), which is also a motivation of the EVP model formulation herein. Moreover, it is noteworthy that the MCC model yield surface (see also Equation (7)) assume the Von-Mises type circular surface in the π-plane. Thereby, except in the triaxial compression state, in any other loading state, the MCC model or equivalent EVP model with the circular type surface overestimate the stress (see also Islam and Gnanendran [23]; Islam et al. [14]). Furthermore, from the literature, comparing the single surface-based model with the composite surface-based model, it is observed that the first one underpredicts the soil state in the overconsolidated state compared to the latter one.
Moreover, from experimental evidence, we find two phenomena (see Schofield and Wroth [21]). The first one is for the overconsolidated clay, and it is evident that with increases in the overconsolidation ratio (OCR), the strength locus for the overconsolidated clay approach to the zerotension line. The second case is for the normally consolidated clay where the strength locus intersects the critical state line in the mean pressure-deviatoric (p-q) plane. In this paper, to resolve such shortcomings, we also revise the MCC model's single yield surface, and we discuss details of them in the model formulation section.

Numerical Modeling
We assume that the porous media: (i) comprises of two phases (the solid phase and the liquid phase (see also Figure 1), (ii) fully saturated, (iii) obeys the small deformation, (iv) fluids follow Darcy's law, (v) supports the isotropic state, the static equilibrium and the isothermal equilibrium conditions, (vi) individual soil grains and the liquid phase are incompressible. Additionally, we ignore the geochemical effect, interaction forces and dynamic actions. Thereby, we also assume each phase density is constant. Furthermore, we consider the soil mechanics' basic principles (see Terzaghi [1]). First, deformations of solid-phase originate due to the liquid-phase removal and rearrangement of the solid's grain. Second, the total stress is the summation of the load carried by the soil (effective stress) and the fluid (pore pressure) (see also Schofield and Wroth [21]). Introducing Terzaghi's effective stress concept [1], we couple the solid phase and the liquid phase.

Governing Equations
Considering the assumptions above, we present governing equations as follows (see also Bear and Bachmat [24]).
(9) Mass balance equation: where ‫'ݒ݅݀'‬ is the divergence operator. is the total stress tensor. ߩ, ߩ and ߩ ௦ are the total density of the porous medium, the liquid phase density and the solid phase density. is the gravitational acceleration acting along the z-axis (see also Figure 1). ߮ is the porosity of the porous medium.
is the time derivative. ௦ and are the solid phase and the liquid phase velocity, respectively.

Constitutive Assumptions
In Equation (9), can be defined as (see Bear and Bachmat [24]): Here, is the identity tensor. ᇱ is the effective stress tensor and ‫‬ is the liquid pressure. Additionally, in Equations (10) and (11), we assume ߩ and ߩ ௦ are constant. Hence, after the summation of Equations (10) and (11), we obtain: For the two-phases porous media, Bear and Bachmat [24] presented the relation between and ௦ using the relative velocity ( ) term and the relative specific discharge ( ) as follows: where, ∇ is the gradient operator. = ఘ ఓ is the hydraulic conductivity which depends on the fluid phase or fluidity ቀ ఘ ఓ ቁ and the specific permeability tensor of soil (). Additionally, ߩ is the liquid phase volumetric weight (ߛ ) , which also express as = ሾ0,0, ߛ ሿ ் , where 'T' represents transpose. ߤ is the liquid phase viscosity and is the gravitational acceleration. Again, in Equation (13), ‫ݒ݅݀'‬ ௦ ' term is given by (see Bear and Bachmat [24]): Here, ߝ ௩ is the volumetric strain and reads: Assuming the small deformation, we also obtain the strain-displacement relation as follows: where u is the displacement component. Substituting, Equations (16) to (18) into Equation (13), then rearranging, we obtain: It is worth mentioning that in the non-associated flow rule EVP model formulation, we assume that ߝሶ (see Equations (17) to (19)) consists of the elastic component and the viscoplastic component as follows.

Strain Rate Tensor of the EVP Model
The total strain rate ൫ߝሶ ൯ in the non-associated flow rule-based EVP model which is given by (see also Lubliner [25]): where, ߝሶ and ߝሶ ௩ are the elastic strain rate tensor and the viscoplastic strain rate tensor, respectively. We obtain ߝሶ as follows (see also Lubliner [25]): where, ߪሶ ᇱ is the effective stress tensor. ܵ is the fourth-order compliance tensor and written as: Here, ߜ is the Kronecker delta. ‫ܧ‬ is the modulus of elasticity and ߥ is the Poisson's ratio. Again, assuming Perzyna's overstressed theory [16], we obtain ߝሶ ௩ in Equation (21) as: where ߔ is the rate sensitivity function; ⟨ ⟩ is the Macaulay's bracket and ‫ܨ‬ is the overstress function. Additionally, ݂ , ݂ and ݂ are the potential surface, the loading surface and the reference surface, which we discuss in the next section. It is worth noting that if ݂ < ݂ , the geomaterials behave elastically, while ݂ > ݂ , similar material will experience the viscoplastic strain. Additionally, we present the derivation of ߔ in Appendix A. Moreover, we also discuss details derivation of ߝሶ ௩ for the finite element implementation in Appendix B.

Bounding Surfaces of the EVP Model
Constitutive models that adopt the classical plasticity theory, such as the MCC model, generally consider a single yield surface (SYF). The limitations of the SYF can be summarized as follows [15,18]: (i) The SYF separates the elastic domain from the plastic state and forms an elastic state boundary within the yield surface. From the comparison of the experimental data and the model predictions, it is evident that the predicted elastic domain is larger than the observed one.
(ii) For the SYF models, the observed transition from the elastic state to the plastic state is in contrast to the experimentally observed gradual changes in the stiffness.
(iii) The SYF provides limited scopes to exemplify the plastic modulus in the loading direction.
(iv) The SYF model is usually incapable of capturing the proportional loadings. During the last couple of decades, limitations of a single-surface model have opened up a more comprehensive research area. There are several methods to overcome these shortcomings. However, the two most popular theories are the multi-surface plasticity and the bounding surface model (see also Hashiguchi [18]; Yu [15]).
In this paper, for any loading history, we consider three bounding surfaces in the EVP model formulation (see also Equations (25) to (27) and Figure 2). Each surface has two ellipses: ellipse 1 (see also Kaliakin and Dafalias [26]) and ellipse 2 (see also Kutter and Sathialingam [27]). In the composite surface, two ellipses of each surface meet at common tangents and allow control of the shape of surfaces. In Figure 3, we illustrate surfaces where ߰ is the slope of the surface at any point on the potential surface. Additionally, as the deviatoric stress ‫)ݍ(‬ decreases with increases in the mean effective pressure ‫)(‬ , the magnitude of ߰ is negative (see also Figure 3). To predict the overconsolidation effect of clay, the surface shape in the 'wet side' is higher than the 'dry side'. In this paper, assuming the ellipse shape parameter (R) is equal to two, we obtain the extended MCC model equivalent ellipse shape (see also Islam et al. [14]).
Loading surface, In Equations (25) to (27), the suffixes ‫,‬ ‫ݎ‬ and ݈ represent the potential surface, the reference surface and the loading surface, respectively. ‫‬ , ‫‬ and ‫‬ are the intersection of the corresponding surface with the positive mean pressure axis (see also Figure 2). Additionally, ‫‬ = ఙ́ೖ ೖ ଷ and ‫ݍ‬ = ቂ ଷ ଶ (ߪ́ௗ) (ߪ́ௗ) ቃ భ మ are the mean effective normal stress and the deviatoric stress, respectively (see also Schofield and Wroth [21]). It is worth mentioning that we considered ‫‬ ௧ = ఙ ೖೖ ଷ as the total mean stress. Also, ‫ܯ‬ represents the slope of the critical state line and is defined as (see also Prashant and Penumadu [28]): where ߶ is the maximum internal friction angle for any specific stress path. b (0 ≤ ܾ ≤ 1) represents b-value (see also Islam and Gnanendran [23]; Islam et al. [14]). ܾ = 1 and ܾ = 0 designate the triaxial compression and the triaxial extension test, respectively. By changing the b-value, we obtain any stress path in the present model formulation. In Equations (25) to (27), we revise the expression of ‫ܯ‬ to obtain the realistic surface in the π-plane. Additionally, in Figure 4, we present a comparison of the MCC surface and proposed modification in the surfaces (see also Figure 2) with the true triaxial test results on Kaolin clay (see also Prashant and Penumadu [28]). Moreover, Lade [29] presented the relation between b-value and Lode angle as follows (see also Figure 5): In Equation (29), ߠ and ܾ represent, the Lode angle and the b-value ቀ= ఙ మ ିఙ య ఙ భ ିఙ య = ఙ́మିఙ́య ఙ́భିఙ́య ቁ respectively, while ߪ ଵ , ߪ ଶ and ߪ ଷ are the major principal stress, the intermediate principal stress and the minor principal stress, respectively, while ߪ́ represents their corresponding effective stress (see also Prashant and Penumadu [28]).

Image Parameters of the EVP Model
Using the triaxial test, we find the reference state for any clay and at any time. Additionally, the current loading stress state ‫(‬ , ‫ݍ‬ ) is associated with its image stress on the reference surface ‫(‬ , ‫ݍ‬ ) and the potential surface ൫‫‬ , ‫ݍ‬ ൯ through the 'radial mapping rule' (see also Hashiguchi [18]). Moreover, similar to the Modified Cam Clay (MCC) model, we also assume that the projection center is in the origin of the stress space (see also Figures 2 and 3). In this regard, Kutter and Sathialingam [27] reported that separation of the projection center and the stress state origin results in an elastic nucleus. Also, such a nucleus requires additional model parameters that split up the elastic domain from the inelastic domain (see also, Hashiguchi [18], Chapter 7). Therefore, ignoring the elastic nucleus, we obtain the image stress of the loading surface on the reference surface (see also Hashiguchi [18]) as follows: Now substituting Equation (30), into Equation (26) for Ellipse 1 (ߟ < ‫)ܯ‬ (see also Figure 3): Similarly, for ߟ > ‫ܯ‬ and Ellipse 2, we find: Again, for the potential surface, we also find:    In Appendix A, we discuss the derivation of ‫‬ and ‫‬ (see also Figure   2).

Couple Finite Element Formulation
We assume that a soil mass occupies in a domain Ω and its surface is ߁, which is subdivided into a subdomain and segment of the surface during the discretization of space, as ߲Ω and ߲߁, respectively. In the following sequences, to obtain a coupled solution for the governing equations (see Equations (9) to (11), (20)), we use the weak form solution and the Galerkin weighted residual method (see also Zienkiewicz et al. [31]). We use the effective stress relation (see Equation (12)) and the Darcy's law (see Equation (16)) to obtain the coupled hydro-mechanical relationship for the elasto-viscoplastic model. In our model formulation, we have a total of 16 equations. One equation for the equilibrium relation (see Equation (9)). Three equations for the mass balance relations (see Equations (10) and (11)). Moreover, we present six equations for the stress-strain relationship (see Equations (22) and (24)) and six equations for the strain displacement relation (see Equation (19)). Additionally, we have 16 unknowns (six for the effective stress, six for the strain, three for the displacement and one for the liquid pressure). Thereby, our coupled solutions represent a "wellposed" problem definition (see also Bear and Bachmat [24]) Additionally, the element matrices for two phases porous media are given by (see also Zienkiewicz et al. [31]): We also find the global matrix by a sum over the number of elements in the element matrix. Also, in a simplified form, we re-write Equation (35) as (see also Owen and Hinton [32]): = ሾ( ) ିଵ + ሿ ିଵ , = ሾ1,1,1,0,0,0ሿ ் , = ሾܰ ଵ . . . ܰ ሿ (48) In the above equations, ா is the tangential stiffness matrix. ா is the coupling matrix. ா is the load vector. ா is the flux matrix and ா is the fluid conduction matrix. T and T represent transpose and the traction force, respectively. is the elastic constitutive matrix related to (see also Equations (22) and (23)). ௨ = ख ௨ is the strain-displacement matrix, where ख is the tangential operator. ௨ and are the displacement shape function and the pore pressure shape function, respectively. ‫ݐ‬ represents time. ߠ ௗ demonstrates the integration parameter (see also Segerlind [33]). is a mapping vector. is the gradient matrix (see also Owen and Hinton [32]).
It is worth mentioning that for the isoparametric element, the shape function and the interpolation function are identical (see Zienkiewicz et al. [31]; Potts and Zdravkovic [34]). This simplification allows flexibility to consider any arbitrary shape of elements. Additionally, the shape functions are defined in terms of the local coordinates (ߦ, ߟ, ߞ) and the strain interpolation matrix requires global derivatives, with respect to the global coordinates ‫,ݔ(‬ ‫,ݕ‬ ‫.)ݖ‬ To map both coordinates, the chain rule can be applied as follows (see also Zienkiewicz et al. [31]): In Equation (49), ሾሿ is the Jacobian matrix. ܰ ഥ is the shape function for nodal values, where ݅ represents nodal points of elements. In the literature, there are two types of assumptions to account for the shape function for the displacement and the pore pressure. In the first case, it is assumed that both shape functions are identical, while in the second case, different shape functions are adopted for variables groups. For example, the displacement and the pore pressure variables vary linearly in the linear triangular element (e.g., three nodes) and the bilinear rectangular element (e.g., four nodes) (Zienkiewicz et al. [31]; Potts and Zdravkovic [34]). However, in a six nodes triangular element and an eight nodes rectangular element, the displacement and the pore pressure field change quadratically. Additionally, when the displacement varies quadratically, the effective stress change linearly (see also Zienkiewicz et al. [31]). Thereby, there is a variation between the pore pressure and effective stress. To achieve the same order of variation in the primary variables, for the eight nodes rectangular element, the degrees of freedom (DOF) for the pore pressure can be obtained at four corners of the rectangular element. In contrast, for the triangular element, a similar DOF for the pore pressure is calculated only from the apex of the triangle (see Potts and Zdravkovic [34]). Hence, the shape function for the displacement and the pore pressure can be separated. Moreover, to account for the large deformation for the homogeneous porous media (e.g., triaxial creep or relaxation test) and heterogeneous porous media (e.g., embankment with multiple layers), special attention is essential to consider the element type (see also Zienkiewicz et al. [31]). For example, to model the consolidation behavior of porous media comprising sand deposit overlying clay deposit, the sand is modeled assuming no pore pressure DOF at the nodes, which behaves like drained media or nonconsolidating elements. In contrast, clay deposit is modeled as a consolidating element considering fluid pressure degrees of freedom at the nodes (see Potts and Zdravkovic [34]). Additionally, Zienkiewicz et al. [31] demonstrated that when porous media is approaching the undrained limit state, to satisfy the Babuska-Brezzi convergence condition (see Zienkiewicz et al. [31]; Potts and Zdravkovic [34]), the shape function for the nodal displacement and the pore pressure need to be separated. In such a case, the choice of element types is limited, and details can be found in Zienkiewicz et al. [31]. In any other state than the undrained limit state, element type selection is extensive in the finite element simulation practice.
In the Results and Discussion Section, we discuss both the drained and the undrained triaxial tests considering the short term loading and the long-term loading. Additionally, in this paper, for validation of the triaxial test, we use the first-order three nodes triangular element. Moreover, for the long-term prediction (e.g., embankment performance estimation), we consider the second-order six nodes triangular element (see also Zienkiewicz et al. [31]).

Time Integration
For finite element modeling (FEM) of time-dependent porous media (e.g., clay), there are several approaches to discretize the time domain (see also Segerlind [33])). However, for the FEM solutions, the ߠ-method is the simplest method (see Potts and Zdravkovic [34]) and in this paper we also use this method. The value of ߠ ranged in between 0 to 1 while ߠ = 0, 0.5, 1 represent the fully explicit time integration (also known as the forward difference method), the implicit trapezoidal time integration (also known as the Crank-Nikolson rule) and the fully implicit time integration (also known as the backward difference method), respectively (see also Zienkiewicz et al. [31]; Owen and Hinton [32]; Segerlind [33]; Potts and Zdravkovic [34]).
Again, nodal values of the interpolation polynomial (ܰ ) for variables (e.g., displacement, pore pressure) at nodal points varies with respect to individual element's nodal coordinates. For example, for the linear triangular element ܰ has three nodal points. Segerlind [33] presented, Equation (36) for the interpolation polynomial of the nodal values as follows: ܰ ሶ + ܰ = .
For ߠ = 0, 0.5 and 1, we also find a simplified form of Equation (55) (see also Segerlind [33]). Additionally, to demonstrate the effect of ‫ݐ∆‬ and ߠ in solutions, for simplicity, we assume an ordinary differential equation as follows: where ܶ and ܰ are the total time and time step incremental number, respectively. The exact solution of Equation (56) is given by: ‫)ݐ(ݕ‬ = ݁ ି௧ . (57) Moreover, we also find the solution of Equation (57) assuming the backward Euler or implicit Euler and the forward Euler or the explicit Euler as follows, respectively:  (58) and (59), we find that for any value ‫ݐ∆‬ > 1, there will be notable oscillation in the forward Euler or the explicit Euler solution. In Figure 7, we present the effect of the selection of ‫ݐ∆‬ and ߠ in solutions.
Again, the solution of the coupled Equation (55) for the non-associated flow rule-based EVP model herein is more complicated compared to Equation (57). However, the solution scheme of Equation (55) holds a similar degree of convergence challenge to select ‫ݐ∆‬ and ߠ as that evident in Figure 7. Thereby, for the solution of the hydro-mechanical coupled equation, careful consideration is essential to choose ‫ݐ∆‬ for the corresponding ߠ-method (see also Owen and Hinton [32]). For our solution, we assume ߠ = 0.5 and we obtain critical values of ‫ݐ∆‬ following Potts and Zdravkovic [34].

Initial and Boundary Conditions
For two-phase coupled porous mediums, we find initial conditions of primary variables as follows:
Additionally, the boundary conditions are given by: where, and ‫‬ are the initial displacement and the initial porewater pressure. represents the unit normal vector. denotes the traction force. We have two types of boundary conditions. They are the Neumann boundary conditions (e.g., ߁ ௧ ) and the Dirichlet boundary condition (e.g., ߁ ௨ ).
In this paper, for the validation of the non-associated flow rule-based elasto-viscoplastic model, we use the conventional triaxial test (e.g., cylindrical specimen). We discuss the details of the validation in the Results and Discussion Section. For validation, we consider different clay samples. Thereby, the initial conditions of those triaxial samples are also different, which we obtain from the published literature. Depending on the clay sample, the initial stress state (e.g., ߪ ), the initial pore pressure ൫‫‬ ൯, the loading rate (e.g., stress-controlled or strain-controlled test) and the confining pressure are different. We assume that in the initial state, all clay samples are fully saturated. Additionally, we consider an axisymmetric section of the triaxial sample for the EVP model validation. In Figure 8, we present a representative illustration of the initial and the boundary condition of the conventional triaxial sample. Again, both the drained and the undrained tests consist of two stages. They are the isotropic consolidation stage and the shearing stage (see also Figure 8d). Additionally, for both tests, the isotropic consolidation procedure is identical. During the first stage of the triaxial test, on the bottom, we assume that the horizontal and the vertical displacement components are fixed. Additionally, due to the symmetry of the axis, we also restrict the horizontal displacement along the radial or x-axis in the left side boundary. We also apply the confining pressure along the right boundary. In addition, depending on the stress-controlled or the strain-controlled test, we also introduce the corresponding value at the top. Moreover, in the initial state, we allow the drainage at the top boundary. Also, we invoke the "geostatic" stress condition considering the body force of the sample. The "geostatic" condition ensures the equilibrium state of the sample. Furthermore, we consider the stress ratio along the horizontal direction to the vertical direction is 1. Then, we run the axisymmetric triaxial sample for one incremental step, which we assume as the initial state for the shearing stage. The initial state makes sure that the triaxial sample satisfies the initial yield surface of the elasto-viscoplastic model.
Additionally, in the shearing stage, there are two essential criteria for the drained and the undrained tests. The first one is the boundary condition, and the latter one is the loading rate. For the drained test, the top and the bottom boundaries are permeable, while these boundaries are impermeable for the undrained test. Additionally, for the drained test, the loading rate is slow to avoid the development of the excess pore water pressure, while in the undrained test, the loading rate is fast. It is noteworthy that for the drained condition, selection of the time increment is crucial, and it needs to be small (see also Figure 7). After the selection of the optimal time step for the drained test, we plot the excess pore water (EPW) pressure with respect to time, confirming that the EPW is approximately zero. Moreover, in the undrained test, immediately after the incremental load application, the excess pore water pressure builds up. Therefore, initially, we consider a small-time step for the undrained test, and then we increase the time step incrementally.
In Figure 8d, we discuss two stress paths for the cylindrical sample. They are the triaxial compression (TC) and the triaxial extension (TE). After the isotropic consolidation state, we apply ߪ ଵ , ߪ ଶ and ߪ ଷ for corresponding TC ad TE tests. It is worth mentioning that the EVP model presented in this paper is not limited to any specific geometry shape or the specific stress paths, which are in most cases, the limitations of many elasto-viscoplastic models in the literature. In the present paper, by changing the ܾ-value, we will able to obtain any stress path in the stress space. Additionally, we will able to introduce a specialized stress path too. For example, the constant mean pressure test or the constant ߪ ଷ test (see also Figures 4,5 and Lade [29]). For such a special loading condition, we need to revise the corresponding ߪ ଵ , ߪ ଶ and ߪ ଷ with respect to the mean pressure, the deviatoric pressure and the ܾ-value. From three equations, we obtain three unknowns, ߪ ଵ , ߪ ଶ and ߪ ଷ .

Model Parameters
In this paper, the non-associated flow rule-based elasto-viscoplastic model requires a total of seven parameters. They are divided into (i) consolidation parameters, (ii) strength parameter, (iii) elastic property, (iv) state parameter, (v) creep parameter and (vi) surface shape parameter. The consolidation properties are the normal consolidation line gradient ቀߣ = ଶ.ଷ ቁ and the swelling line gradient ቀߢ = ೞ ଶ.ଷ ቁ, where ‫ܥ‬ and ‫ܥ‬ ௦ are the compression index and the swelling index, respectively (see also Figure A1 in Appendix A). The slope of the critical state line ‫)ܯ(‬ is considered as the strength parameter and related to the angle of internal friction at failure (see Figures 2,5 and Equation (28)). We consider the Poisson's ratio (ߥ) as the elastic property. Additionally, at any reference time and unit mean pressure, the void ratio (݁ ே ) is the state parameter (see also Figure A1 in Appendix A). Moreover, to account for the time-dependent behavior of clay, we introduce the secondary compression index ‫ܥ(‬ ఈ ) as the creep parameter. Finally, we also introduce a surface shape parameter (ܴ) to control the shape of the bounding surface (see Figures 2 and 3 and Equations (25) to (27)). In Table 1, we present a summary of the non-associated flow rule-based elasto-viscoplastic (EVP) model parameters in this paper and their determination method.  It is noteworthy that among seven parameters, five parameters are identical to the Modified Cam Clay (MCC) model. The other two parameters are ‫ܥ‬ ఈ and ܴ. Additionally, the MCC model parameters determination processes are well documented in many soil mechanics textbooks (see also Schofield and Wroth [21]; Roscoe and Burland [10]). Moreover, there are three approaches available in the literature to obtain ‫ܥ‬ ఈ . The first one is the empirical relation with ‫ܥ‬ (see also Mesri and Castro [35]) which is assumed constant over time. The second one is a direct calculation from the oedometer test or the triaxial test. This approach is also divided into the void ratio-based method ቀ‫ܥ‬ ఈ = ∆ ∆ ௧ ቁ and the strain-based method ቀ‫ܥ‬ ఈఌ = ∆ఌ ∆ ௧ = ഀ ଵା బ ቁ (see also, Liingaard et al. [11]). Additionally, using the cone-penetration test field data, we also may obtain ‫ܥ‬ ఈ (see also Tonni and Simonini [36]). Moreover, in the literature, there are two concepts regarding the secondary compression index: constant or linear function and non-linear function, while in both cases, strong opinions are available. In this paper, we use the void ratio-based non-linear ‫ܥ‬ ఈ which is neither tied to any specific clay nor requires any fittings parameters (see also Islam and Gnanendran [23]; Islam et al. [14]) and reads: where ݅ and (݅ − 1) represent the present step and the previous time step, respectively. Additionally, to obtain the surface shape function (ܴ), there are several methods. For example, Dafalias and Herrmann [37] proposed a fitting method using the undrained triaxial stress path as follows: It is worth mentioning that to predict the pore pressure for the undrained test, we also need the permeability (see Equation (16)), which is the material parameter of clay and not relevant to the EVP model parameters. Additionally, in Table 2, we summarize EVP model parameters for different clays.

Shanghai Clay
From Huang et al. [38], we obtain the Shanghai soft clay properties and model parameters. It is an undisturbed soft sensitive clay (sensitivity = 4.86). The water content, the liquid limit and the plastic limit of this clay sample are 51.80%, 44.17% and 22.40%, respectively. Additionally, the specific gravity of the Shanghai clay is 2.74. Moreover, the clay fraction (< 5 ߤ݉), the silt fraction and the sand fraction of this clay are 26.60%, 63.40% and 10.00%, respectively. Furthermore, Huang et al. (2011) demonstrated that diameter and length of triaxial samples were 39.10 mm and 80.00 mm, respectively, which we use for the preparation of the numerical model geometry.
In Figures 9, we compare numerical results with laboratory-measured stress controlled, isotropic and consolidated undrained triaxial tests (see also Huang et al. [38]). Additionally, we also compare the EVP model predictions with the Kutter and Sathialingam [27] proposed EVP model. From the comparison, we find that in the small strain zone and after 14% axial strain, both the EVP models over predict the deviatoric stress. It is to note that in the EVP model formulation, we ignore the hysteretic response of the clay (see also Whittle and Kavvadas [41]), which results in the overprediction in the early stage. Also, in the present model formulation, we did not include the destructured behavior of clay (see also Liu and Carter [42]). Therefore, we observe a slight overprediction at the higher axial strain. Inevitably, in the present EVP model formulation, the incorporation of the hysteretic response and the destructured response require extra model parameters.
In Figure 9a, we compare the associated flow rule and the non-associated flow rule EVP models predicted stress-strain responses with the Shanghai clay observed results. We find that for 150 kPa mean pressure up to 1% axial strain; both models capture experimental responses well. Then, before the measured peak deviatoric stress, we also observe overprediction. Afterward, under-prediction in the associated flow rule-based EVP model continues. We also observe similar results for 200 kPa. Comparing both the flow rule EVP model and the Kutter and Sathialingam [27] model, we also find that the non-associated flow rule model well captures the experimental results.
Furthermore, we present the stress path responses in Figure 9b. We observe that after attaining the peak deviatoric stress, the EVP model prediction gradually follows a 'narrow region' (see also Islam et al. [14])) and such a phenomenon is common in the natural clay than the reconstituted clay (see also Islam and Gnanendran [23]). In general, the MCC model (Roscoe and Burland [10]) and the associated flow rule EVP models based on the MCC framework (see also Kutter and Sathialingam [27]) are incapable of capturing the 'narrow region'. To predict such a behavior, incorporation of the additional model parameters in the MCC framework are also frequent (see also Liu and Carter [42]). However, in this paper, we obtained such a 'narrow region' of the natural clay without any extra model parameters.

San Francisco Bay Mud Clay
From Lacerda [39]; Kaliakin and Dafalias [26], we obtain EVP model parameters for the San Franciso Bay Mud (SFBM) natural clay. In this section, we compare the EVP models' predictions with the experimentally observed undrained triaxial test results for the relaxation test. Additionally, we demonstrate herein SR-I-5 test data (see also Lacerda [39]). The water content, the liquid limit, the plastic limit, and the plasticity index of this clay sample are 88-93%, 88.4-90%, 35-44% and 45-55 %, respectively. Moreover, the specific gravity of the SFBM clay is 2.66-2.75. Also, the isotropic consolidation pressure of the clay sample is 78.4. Additionally, the initial void ratio and its equivalent mean pressure are 1.30 and 156.9 kPa, respectively (see also Kaliakin and Dafalias [26]). In Table 3, we present the shear and the relaxation phase, axial strain, strain rate and the duration of the test. In Figure 10, we present a comparison of the measured data, and EVP models predicted responses. For the EVP models, we also compare the flow rule effect considering the associated flow rule (AFR) and the non-associated flow rule (NAFR). We observe that the NAFR based EVP model capture well the experimental response. It is worth mentioning that Islam et al. [14] also find an identical flow rule effects for the extended Modified Cam Clay equivalent surface.  In Figure 11, we present a comparison of the EVP models' predicted response with the observed results of the drained triaxial compression test on the San Francisco Bay Mud clay. The triaxial sample was isotropically consolidated to the confining pressure 156.9 kPa. Additionally, the axial strain rate was 0.0031%/minute. From Figure 11, we observe that the non-associated flow rule EVP model well captured the experimentally observed results.

Kaolin Clay
In this section, we present the over consolidation ratio effect of the Kaolin clay (see also Herrmann et al. [40]) for the undrained triaxial compression and extension tests. It is a reconstituted clay and comprised of Snow-Cal 50 Kaolin and 5% bentonite mixture. The specific gravity of this clay is 2.64. The liquid limit and the plasticity index of this clay are 47.0% and 20.0%, respectively (see also Herrmann et al. [40]).
The initial isotropic consolidation pressure of the triaxial Kaolin clay sample was 392.2, and corresponding the initial void ratio for the unit over consolidation ratio (OCR) was 0.613 (see also Islam and Gnanendran [23]; Islam et al. [14]). We calculate the confining pressure for different OCR'S considering the initial state condition for both the triaxial compression and extension tests.
From Figure 12 and for OCR = 1, we observe that the non-associated flow rule-based EVP model prediction of the stress-strain response is satisfactory before attaining the peak. Then, after 14.0% axial strain, the under-prediction in the non-associated flow rule-based EVP model is 1.05%. Additionally, such a magnitude comparing with experimental results for the associated flow rule are 3.5% and 1.40%. We also observe a similar prediction for the undrained triaxial compression with OCR = 2, 4 and 6 for the stress-strain responses. For the triaxial extension test, we present a similar comparison for OCR = 1 and 2. From such a contrast of EVP models with the measured experimental results illustrate the effect of flow rule in the EVP models predictions.
Additionally, for OCR = 1 and 2, in Figures 12 b and 12 c, we present the pore pressure response and the stress path response, respectively considering the triaxial compression and the triaxial extension test. We also find that the non-associated flow rule prediction is close to the experimental responses.

Application of the EVP Models
We use the EVP models in this paper to capture the long-term behavior of the Nerang-Broadbeach Roadway (NBR) embankment in Australia. Islam et al. [20] presented details of geotechnical properties, subsurface and geology, the geometry of the NBR embankment sections, measured instrumentation data for the settlement plates and piezometers, methodologies to determine the model parameters and laboratory, as well as field measured values. In this paper, we compare 590 days measured data of a surcharged-preloading section with the EVP models' predicted response considering the associated and the non-associated flow rule to illustrate the flow rule effect. Additionally, we also compare the Modified Cam Clay (MCC) predicted results with the EVP models' response to demonstrate the viscous response of the clay.
We present the geometry of the finite element section for the NBR embankment in Figure 13. Additionally, to reduce the boundary effect in numerical simulations, we extend the width and depth of the embankment. Moreover, to avoid the instabilities in the finite element simulations, we implement surcharge preloading incrementally. Furthermore, 3.0 m preloading was applied in 370 days. Then, 1.0 m additional surcharged load added and monitored for 220 days. We demonstrate the model parameter in Table 4. Moreover, in Figure 14, we present a comparison of the observed and FE predicted responses. We observe that the non-associated flow rule model well captured the measured settlement plate response compared to the associated flow rule model. Additionally, from the comparison of the MCC model with the EVP models, it is evident that the MCC model underpredicts the long-term behavior of the embankment. Such underprediction in the MCC model developed due to the exclusion of the viscous behavior of clays. Additionally, we also observe that for the long-term monitoring of EVP models, the flow rule effect is significant. We also noticed a similar response in the Results and Discussion Section for natural clays.

Conclusions
In this paper, we presented a non-associated flow rule (NAFR)-based three surfaces elastoviscoplastic model. Assuming the potential surface and the reference surface are identical, the NAFR Settlement (mm) model is also reduced to the two-surfaces associated flow rule mode. Both EVP models' surfaces comprises two ellipses, while the surface shape parameter (ܴ) controls each surface shape. Again, considering ܴ = 2 , we also find the extended Modified Cam Clay equivalent surface shape. Additionally, to obtain a realistic non-circular surface in the ߨ-plane, we also included the ܾ-value to define the slope of the critical state line. Thereby, by changing the ܾ-value, any stress path can also be attained.
It is worth mentioning that EVP models in this paper require only seven parameters. Also, all model parameters are well known in textbooks, and all of them can be deduced from simple laboratory tests. Additionally, in the present models' formulations, we did not consider any fittings parameters. It is to note that EVP models in this paper are formulated for the isotropic clay, and any other state additional model parameters with associated modification are essential. Moreover, EVP models herein are not capable of capturing the structured behavior of highly sensitive clay, which require further changes with extra model parameters. Also, the present model may not be the best fit for the dynamic performance of clay. Furthermore, the EVP model herein is limited to the Darcy type and immiscible fluid flow.
For validation of EVP models, we considered a wide range of triaxial tests for natural and reconstituted clays. From comparisons of observed and predicted responses, we found that the nonassociated flow rule EVP model well captured experimental results. Also, we found that the flow rule effect is noticeable for the natural clays than the reconstituted clay. Additionally, we implemented EVP models in a coupled consolidated finite element solver. Then, we used EVP models and the Modified Cam Clay (MCC) model to predict long-term monitoring data of the Nerang Broadbeach Roadway embankment. We also found that the non-associated flow rule EVP model well captured the settlement plate measured response compared to the associated flow rule EVP model and the MCC model. From Figure A1, we also find (see also Roscoe and Burland [10]): Now, combining Equations (A1) to (A3), we find: In the triaxial stress space, ߝሶ ௩ ௩ also can be written as: By comparing Equations (A4) and (A5), then separating ߔ for calculation of the viscoplastic strain component results in the flow rule independent elasto-viscoplastic model, which violates the flow rule theory (see also Islam and Gnanendran [23]). Therefore, for the non-associated flow rule, to obtain ߔ, we find ப డ considering the one-dimensional compression test criterion for relevant flow rule as follows: Again, from Figure A1, we also find expression for ‫‬ and ‫‬ (see also Islam et al. [14]):