A simple hypoplastic model with loading surface accounting for viscous and fabric effects of clays

This paper presents a simple hypoplastic model capturing mostly all salient features of clays: rate dependency, time dependency and inherent and induced anisotropy without being restricted to only viscoplastic clays. Therefore, due to the strain rate decomposition into three parts, nonviscous clays, that is, rate‐independent clays, can also be simulated. The incorporation of a loading surface allows to capture the behaviour of normal and overconsolidated clays. The model requires eight material parameters, which are simple to calibrate from standard laboratory tests. In total, 77 simulations of five different clayey‐like soils are compared with experimental data. The simulations contain one oedometer test with loading–unloading–reloading cycles, creep and relaxation stages, both undrained and drained triaxial tests in compression and extension, as well as eight incremental response envelopes capturing also the directional response of Beaucaire Marl clay. Some limitations of the model such as the description of temperature effects on the behaviour of clays are also pointed out.


INTRODUCTION
There are a number of important geotechnical engineering problems-notably the analysis of soil-structure interaction in deep excavations, shallow foundations and tunneling-where large differences typically exist between stress paths experienced in different regions of the surrounding soil, both in terms of magnitude and direction [1][2][3][4] and in terms of time for most fine-grained soils. 5,6 Furthermore, the question arises if the fine-grained soil is normally consolidated or overconsolidated. In these cases, the quality of the engineering prediction crucially depends on the ability of the adopted constitutive model to correctly describe the behaviour of the soil along a wide range of loading paths and rates. In addition, the important challenges are the anisotropy and bonding/destructuration, which influence the relation between strain and stress. Also, the in situ conditions in terms of natural and reconstituted samples or drainage of the soil effect remarkably the behaviour and stiffness of fine-grained soils. [7][8][9][10][11][12][13] Advanced models taking into account these different features are expected to simulate the behaviour of clays accurately. Both the number of parameters and the number of equations necessary for the model description are in most of these cases high and complex.
Up to now, a relatively thorough understanding was implemented in constitutive models for time-independent fine-grained soils. [14][15][16][17][18][19][20][21][22][23][24][25][26][27][28] Some of these works are restricted to either normally consolidated 22 or overconsolidated 21 fine-grained soils. Some models need various parameter sets for different overconsolidation ratios, which for the users is very confusing. In order to capture the strain-rate dependency typical for plastic clays, different models [27][28][29][30][31][32][33][34][35][36][37] have been developed. The majority of recent works in this field have opted for the viscoplastic strains to replace the original plastic strains in order to describe the time and strain rate effects observed in experiments with soft clays. Therefore, for the simulation of a nonviscous clay, another constitutive model should be used instead and the user may lose the connection between the models' mathematical formulations. Apart from this point, most of the models dealing with rate-dependent clays do not focus on the incorporation of the effects of soil structure and fabric anisotropy. Recently, a very few number of models for clays accounting for fabric have been developed, 32,35,[38][39][40][41] but none of them use directly the simple and pioneer theory of transversal isotropy for clays of Graham and Houlsby. 42 Whereas the extended Creep-SCLAY1 model 32 includes the effects of structure introducing four additional parameters, the isotach elastoplastic model presented in Yang et al. 35 introduces six additional material constants to capture the fabric of the material. On the other side, the model introduced by Masin and Rott 3 and Masin 43 introduces three additional material parameters while using transversal isotropic elasticity. Using the theory of Houlsby, only one additional material parameter is added to the proposed model in this work for the description of main effects resulting from fabrics. An evolution of this parameter with the hypoplastic strain is proposed.
Therefore, this paper focuses on the development of a simple hypoplastic model accounting for mostly all the salient features clays exhibit. In contrast to most conventional models, the time dependency is simulated by an additional strain mechanism, namely, the viscous strain rate, whose intensity is directly related to the secondary compression coefficient C (Buisman constant 44 ). The pure hypoplastic strain rate is responsible for the description of the intrinsic (time-independent) behaviour of the material. Both the viscous and the hypoplastic strain rates are influenced by the overconsolidation ratio accounting for the behaviour of both normally consolidated and overconsolidated clays. The structure of the material is introduced by transforming the elastic stiffness incorporating an additional material parameter according to the theory of Graham and Houlsby. 42 First, the mechanical formulation of the proposed model in the triaxial space is presented (please note that in this section, the inherent anisotropy is not considered). Then, the generalization of the triaxial formulation to the multiaxial stress space is described systematically for each equation considered, introducing herein also the fabric. Subsequently, an inspection into the performance of the model has been done. In total, 77 simulations of five different clayey-like soils are compared with the experimental data. The simulations contain one oedometer test with loading-unloading-reloading cycles, creep and relaxation stages, both undrained and drained triaxial tests in compression and extension, as well as eight incremental response envelopes capturing also the directional response of Beaucaire Marl clay. The model requires the calibration of eight material parameters, which are simple to calibrate from standard laboratory tests. Some limitations of the model such as the description of the dilatancy for normally consolidated clays and the incorporation of temperature effects are also pointed out. Finally, a finite element calculation of a shallow foundation is carried out with Abaqus to evaluate the dependence of the undrained shear strength on the displacement velocity and on the anisotropic coefficient. The major part of the article is also presented in the PhD thesis of the first author. 45 The notation and convention of the present work is as follows: italic fonts denote scalar magnitudes (e.g., a,b), bold lowercase letters denote vectors (e.g., a,b), bold capital or greek letters denote second-rank tensors (e.g., A, ) and special fonts are used for fourth-rank tensors (e.g., E,L). Indicial notation can be used to represent components of tensors (e.g., A ij ), and their operations follow the Einstein summation convention. The Kronecker delta symbol is represented by ij , that is, 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 . Multiplication with two dummy indices (double contraction) is denoted with a colon ":" (e.g., A : B = A ij B ij ). The symbol "⊗" represents the dyadic product (e.g., A ⊗ B = A ij B kl ). The brackets || ⊓ || extract the Euclidean norm (e.g., ||A|| = √ A i A i ). Normalized tensors are denoted by ⃗ ⊓ = ⊓ ||⊓|| , or in general as ⊔ → . The superscript ⊓ * extracts the deviatoric part of a tensor (e.g., A * = A− 1 3 (trA)1). Components of the effective stress tensor or strain tensor in compression are negative.
= 2∕3 | 1 − 3 |; and the Lode angle cos (3 ) = − √ 6 tr ( * · * · * ) ∕( * ∶ * ) 3∕2 triax. = −3 where the subscripts 1 and 2 and subscript 3 denote the axial and the horizontal directions, respectively. The stress ratio is defined as = q/p. The deviator stress tensor is defined as * = +p 1 and the stress-ratio tensor with r = * ∕p = The relation between the stress rates and the strain rates is proposed to The total strain rates consist of the elastic strain rate with the viscosity index I v , whereby the well-known Buisman constant 44 C = I v represents the intensity of creep, that is, of the viscous strain rate. The overconsolidation ratio is denoted as OCR and the degree of nonlinearity with Y and the volumetric and deviatoric part of the flow rule with m v and m q , respectively. The bulk modulus K is adjusted to conform the behaviour of clays at normally consolidated states. Therefore, we consider the model formulation for isotropic compression, neglecting, for instance, the viscous strain rate: The isotropic compression . v > 0, e = e i0 − ln(p∕p re ) with the maximum void ratio e i0 at the reference mean stress Combining Equation (4) with Equation (3), we obtain Along the same lines of thoughts, the (minimum) isotropic value of the degree of nonlinearity Y 0,max is formulated considering isotropic unloading with . v < 0, e = e i0 − ln(p∕p re ). Using the evolution equation of the void ratio and the constitutive equation considering isotropic unloadinġ we obtain for Y 0,max The material parameters and represent the compression and swelling index in the e − ln(p) space, respectively. The elastic shear modulus is expressed in terms of K using the Poisson ratio as a material parameter: Note that using these relations, any hypoplastic model for sand can be reformulated to account for the behaviour of soft soils (without viscosity if the viscous strain rate is dropped I v = 0). Of course, the definition of void ratios need to be adjusted as well.

Flow rule m
For the formulation of the flow rule m, which describes the direction of both plastic and viscous strain rate in this model, first the bounding surface (BS) in conjunction with the critical state surface (CSS) needs to be defined. The CSS characterizes the vanishing stress rates . p = . q = 0. We introduce the relation F c (p, q) = 2 − M 2 = 0 in triaxial space for the CSS. A comparison between F c , the CSS defined by Matsuoka-Nakai (M-N) and the one of Modified Cam Clay (MCC) is given in the principal stress space in Figure 1A,B.
The bounding surface F b (p, q) = 2 − M 2 b = 0 with its triaxial slope M b = f b M c intersects the isotropic axis q = 0 at the maximum value of the void ratio e = e i (p); hence, at the Hvorslev pressure, 47 p ! = p ei = exp((e i0 − e)∕ ). For this purpose, the following formulation has been adopted for f b 37,48 : The exponent n f can be adopted to fulfil the requirement that the bonding surface intersects the CSS at e = e c : The BS and CSS of the proposed model are illustrated in the principal stress space in Figure 1C,D and will be described in detail in Section 3.2.
Analogous to the two stress surfaces (critical and bounding), the model introduces two characteristic void ratios presented in Figure 1E, which are pressure dependent. Associated to the CSS is the critical void ratio: whereby e c0 is the critical void ratio e c at p ref = 1 kPa. According to the conventional MCC theory, 49 the bounding surface intercepts the critical state at p = p ei ∕2 = exp((e i0 − e)∕ )∕2. To fulfil this requirement, the reference critical void ratio e c0 at the reference pressure of p ref = 1 kPa has been adopted to The second characteristic void ratio, namely, the maximum void ratio at isotropic compression (normal consolidation) is associated to the BS and reads e i = e i0 − ln(p∕p re ), p re = 1 kPa, whereby the material parameter e i0 is the maximum void ratio of an isotropic virgin compression path at p ref = 1 kPa. The postulate of a maximum void ratio is reasonable as the soil cannot take arbitrarily large void ratios without considering the effect of macropores, which is not subject of the present work. For isotropic states lying on the BS, e = e i holds, and for those lying on the CSS ,e = e c is implied, illustrated also in Figure 1F. In order to simulate the behaviour of sand, the formulations of e c and e i can be adjusted similarly. Equation (29) allows the developer to define the normalized flow rule ⃗ m describing the direction of plastic flow by having defined a CSS. Experiments have shown that the direction of plastic flow at the critical state corresponds to the deviatoric direction and at isotropic states to isotropic direction. Thus, an associated flow rule m = F c ( ) ∕ can be very easily deduced. Considering the critical state surface in the p−q space F c (p, q)= −M c g with = q/p and g = 1 or g = c = M e /M c for triaxial compression or extension, respectively, the associated flow rule reads The disadvantage of an associated flow rule for hypoplastic models is documented by other authors as well. 17,22,30,50 In fact, an associated hypoplastic flow rule yields to a high overestimation of the K 0 value so that asymptotic states cannot be accurately reproduced; see the detailed description provided by Niemunis 30 and Wu and Niemunis. 51 Furthermore, as documented by Masin, 17 models with nonrotated state boundary surface (as is the case with the present model) cannot describe the behaviour of anisotropic clays accurately. A significant improvement in this direction can be achieved by a nonassociated flow rule (even though the BS remains nonrotated). Therefore, a reformulation of the flow rule according to the aforementioned experimental attributes ( with the Lode angle function g( ) described in Section 3.2, Equation (27).

Degree of nonlinearity Y and overconsolidation ratio OCR
The degree of nonlinearity herein is proposed to take into account the isotropic value Y 0,max , which corresponds to the minimum value of Y and is deduced in Section 2 to capture the slope of the isotropic unloading line in the e − ln(p) space starting from a normally consolidated state. If the unloading starts from an overconsolidated state OCR > 1, Y 0,max needs to be extended to Y 0 = Y 0,max (1/OCR) 2 , a slightly different but more consistent formulation than the former one proposed in Fuentes et al. 37 (Y 0 = Y 0,max (p/p ei ) 2 ), because it incorporates not only the OCR describing one-dimensional states but the general formulation of OCR considering also stress states with q ≠ 0. In order to fulfil the requirement of constant stress at the CSS . p = 0, . q = 0, the degree of nonlinearity should take its maximum value Y = 1 at the bounding surface. An appropriate function for this purpose is the following linear interpolation: whereby for isotropic states, = 0 ⇒Y = Y 0 and at the BS, For the determination of the viscous strain rate, we need a precise definition of the OCR considering the following requirements for isotropic states (q = 0): • at p = p ei ⇒ e = e i , the stress state is normally consolidated (OCR=1) and lies at the bounding surface; compare Figure 1E and Figure 1F. Thus, the maximum isotach is reached at OCR=1, which corresponds to . → ∞.
• p = p ei /2 ⇒ e = e c is the intersection point between the BS and CSS; see Section 3.2 and Figure 1E,F.
• the material parameter e i0 is the maximum void ratio at a reference mean pressure of p ref = 1 kPa. It can be calibrated from an isotropic compression test or an oedometer test with constant rate of strain as shown in Appendix B1.
According to these assumptions, the relation OCR ≥ 1 holds for the model. It is debatable whether OCR < 1 is realistic for soils. However, from the authors' point of view, the memory of the soil would rearrange in every occasion its overconsolidation ratio so OCR ≥ 1 seems to be more realistic.
In the isotropic axis, the simple relation OCR=p ei /p is adopted. For deviatoric stress states, a similar schema to the one proposed in previous works, 28 with the simplification A = 1 − ( M c g b0 ) 2 , the following expression is obtained: Hence, the OCR can now be calculated through In order to keep the mathematical formulation of the model as simple as possible, OCR can be formulated analogously to the degree of nonlinearity Y. Comparing the pure (hypo)plastic strain rate and the viscous one, then Y represents the counterpart of (1∕OCR) 1∕I v . Having defined the isotropic value of OCR=p ei /p, for general cases, the formulation is straightforward regarding the aforementioned requirements and the relation given in Equation (16): whereby an additional material parameter n OCR has to be introduced, which allows the user more flexibility for the calibration process, but on the other side, it complicates the numerical implementation. Figure 2A,B illustrates the OCR determination after Equation (19) and Figure 2C the definition proposed in Equation (20). The coincidence with the BS is evident (BS blue line, OCR green line). However, taking into account that Equation (19) does not introduce any additional material parameter, it represents a more elegant method for determining OCR. This relation is implemented in this work.  Figure 2 were overlapped and has now been corrected in this version.] Figure 2D represents for isotropic state q = 0, the evolution of the mean stress p and the Hvorslev mean pressure p ei . Obviously, a difference between these evolution rates is present for unloading, because p ei evolves per definition with slope and p with slope .

MECHANICAL FORMULATION IN THE MULTIAXIAL SPACE
This section provides the generalization of the triaxial formulation to the multiaxial stress space. Systematically, each equation from the triaxial space is considered and adopted to the multiaxial space, whereby herein also the description of fabric is introduced by transforming the (hypo)elastic stiffness tensor to a cross isotropic one.
The stress rate is related to the strain rate by the anisotropic stiffness tensor E trans :

Elastic stiffness tensor
The elastic stiffness tensor adopted herein is hypoelastic and follows similar assumptions and relations as in previous works 37,52-54 : which represents the multiaxial generalization of the stiffness introduced in Equation (1). The last term in Equation (22) corresponds to the off-diagonal terms in Equation (1) and considers the reduction of the stiffness close to critical states and the rotation of incremental response envelopes for initial states lying near the CSS. The Bulk modulus K is deduced in Section 2, Equation (5). The shear modulus G is expressed in terms of K and and holds the same as for triaxial space; see Equation (8).
It is well known that hypoelastic constitutive relations cannot be expressed in terms of an elastic strain energy function. Thus, the response is path dependent, and dissipation may occur even though the material is supposed to be elastic. However, the main point is that the elastic strains are assumed to be relatively small so that any error in the conservation of energy is very small, and so it can be neglected. The stress rate . must be objective; thus, the material derivative of the Cauchy stress cannot be used. Yet, other objective rates such as the Jaumann, Green-Naghdi, and Truesdell can be used instead. Usually in hypoplastic constitutive equations, as well as for the proposed model, the Jaumann rate is used.
To overcome these shortcomings and in order to obtain a conservative constitutive model for clays under purely elastic conditions, a hyperelastic stress-strain relation derived from an elastic potential is being developed by the authors.

Incorporation of fabric effects for clays
Most natural clays show anisotropic behaviour because of their mode of deposition and the elongated shape of the particles. This anisotropy, known as inherent anisotropy or fabric, resulting from the deposition process, which tends to induce a horizontal bedding plane in the soil layer, can be regarded as transverse isotropy. 42 A scanning electron microscope (SEM) study has been performed recently 55 on a Kaolin clay. Simulations of this material will be shown later on in Section 5. From the results of the SEM (see also Figure 3), a preferred direction of the particles' orientation can be observed. At a scale of 50 μm, the surface structure is oriented ( Figure 3A). Even at a scale of 200 μm, the particles tend to be aligned in one direction and to form somewhat stronger particle stacks (see Figure 3B).
This study witnesses the theory proposed in Graham and Houlsby 42 at least for this clay. Other SEM studies reported in, for example, Bohac et al. 56 investigated the fissuration of a Neogene clay (Brno clay) and found out that the fissuration and even its laminated fabric may be explained by the anisochoric transition from the marine material into clay. The inclination of the stress paths and the stiffness of Brno clay tested in Bohac et al. 56 are very similar to those reported for Kaolin clay in Wichtmann and Triantafyllidis, 11 which are simulated with the proposed model in Section 5. They both indicate transversal isotropic behaviour of the material.
This type of fabric influences the elastic stiffness of the material. Even though five elastic parameters are needed to describe the transversal isotropic elastic stiffness, Graham and Houlsby 42 deduced a simplified transversal isotropic elasticity for the special case of anisotropic clays, whereby only three parameters are needed: the Young modulus E = E v , the Poisson ratio = h and the scalar factor governing the following relations between the horizontal and vertical stiffness: introducing only one additional material parameter into the proposed model. With the y-axis as the vertical direction representing the direction of anisotropy and the x z-plane as the plan of isotropy (see Figure 3C), the stress-strain increment equation for a transversal isotropic material considering the simplifications provided by Graham and Houlsby 42 can be expressed through the relation given in Equation (24), whereby E trans is the transversal isotropic stiffness.
To incorporate this stiffness into a constitutive model, the present stiffness (Equation 22) needs to be adjusted, scaled and rotated. After some mathematical manipulations, one may notice that any hyperelastic stiffness E renders a transversal isotropic one by a special scaling transformation as denoted in Niemunis et al. 57 : whereby m s is the unit vector along the sedimentation axis; for example, if the sedimentation axis is vertical then m s = {0, 0, 1}. E ijkl is the isotropic elastic stiffness tensor from Equation (22) and the material parameter is herein denoted as the anisotropic coefficient. Some recent research works have shown the importance of the inherent anisotropy on the behaviour of clays. According to the investigations taken on Kaolin clay by Wichtmann and Triantafyllidis, 11 the samples cut out in horizontal direction ( Figure 3C) and show a more dilative response and a higher undrained shear strength than those taken vertically. Furthermore, depending on the cutting direction-if all other conditions remain the same-the horizontal samples can withstand a much larger number of stress cycles to failure than the vertical samples. 11 This postulates the importance of simulating the different elastic stiffness for different cutting directions.
Furthermore, an extensive review about the small strain stiffness anisotropy of natural sedimentary clays has been performed by Masin and Rott 3 and implemented by Masin. 43 In their, work the assumption and simplification of Graham and Houlsby given in Equation (23) is neglected and instead five material parameters are used. Nevertheless, we will show that even using Equation (23) in order to simplify the constitutive model and reduce the number of material parameters, the behaviour of some clays is reproduced very well.

Definition of the characteristic stress surfaces and void ratios
In this work, two characteristic surfaces are introduced, whereby the first one corresponds to the CSS. The CSS incorporated into the present model is defined through a general and very simple relation, which has first been proposed in Fuentes et al. 37 and is now modified to conform the triaxial condition q c = p c M c g( ) of the Mohr-Coulomb criterion: with the Lode angle dependency g( r ) of the stress ratio r. In general, the function g( ) introduced in the model is borrowed from Argyris et al. 58 : whereby c = M e /M c = 3/(3 + M c ) represents simply the ratio between the critical state slope for triaxial extension and compression according to the Mohr-Coulomb relation. Some researchers incorporate the Matsuoka-Nakai criterion as a critical state surface. 17,30 This is rendered as a special case of the introduced critical state surface for M( )=M c g( ) as the smallest positive solution of the equation 2∕27(9 + 8tan 2 ( )) cos(3 )M 3 + 1∕3(−6 − 8tan 2 ( ))M 2 + 8tan 2 ( ) = 0. A comparison between the Matsuoka-Nakai (M-N), MCC and the introduced F c surface is represented in the deviatoric space in Figure 1A,B. The CSS or the yield surface ( ) = 0 in the therms of the Niemunis hypoplastic theory 30 is imposed into a constitutive equation by the condition . = . Neglecting the viscous therm, this condition renders which is satisfied trivially by . = 0 and by The second characteristic surface, namely, the bounding surface (BS) F b ( ) = 0, introduced into this models formulation imposes a bounding of the void ratio e < e i and reads 37 with the triaxial slope of the bounding surface If one wants to consider the dilatancy in the formulation, a dilatancy surface F d ( ,e)=0 can be introduced. However, here, it is left out for improvement in future works.

Flow rule
The unit tensor m denoted as flow rule defines the direction of flow for . = 0. Obeying experimental observations documented in, for example, previous works, 11,17,30,31,49,59 the flow rule should have purely volumetric direction at isotropic stress states and deviatoric direction at the CSS F c ( ) = 0; see Equation (26). An interpolation relation between these states is more general and flexible than other hypoplastic flow rule formulations (e.g., Wolffersdorff 50 ). The flexibility consists in the possibility to simply adjust the flow rule for the description of K 0 paths or to implement any dilatancy rule not affecting the bounding surface or the elastic tangential stiffness tensor. A suitable candidate for this purpose considering also Equation (26) is the following function adopted from previous studies 37,53,54 : with the deviatoric stress ratio r and the critical stress ratio r c defined in Equation (26). This definition of flow rule imposes a consistency of K 0 with Jaky's formula. 60

Degree of nonlinearity and overconsolidation ratio
The generalization of the triaxial formulation to the multiaxial stress space of the degree of nonlinearity is straightforward.
In fact, the triaxial stress ratio will be replaced through the magnitude of the multiaxial stress ratio r, and in order to account for the bounding surface instead of g M b , the norm of the bounding stress ratio r b from Equation (30) is used: The definition of the overconsolidation ratio remains the same as depicted in Equation (19) or (20). The differences between these relations are described in detail in Section 2.2 and illustrated in Figure 2.

NUMERICAL IMPLEMENTATION
The constitutive model was implemented in a Fortran subroutine as "User Material" (UMAT), compatible with the software Abaqus Standard. 61 To avoid numerical issues, a substepping scheme was adopted. Very small substepping sizes in the order of ||Δ ||=10 −5 were selected to assure numerical convergence. In addition, some numerical recommendations for viscous models pointed out by Niemunis 30,31 were considered. 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 detail in Tafili et al. 27 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. During the calculations, no further numerical issues were detected.

INSPECTION INTO THE PERFORMANCE OF THE MODEL COMPARED WITH LABORATORY TESTS
Having defined the constitutive equations of the model, we proceed now with the evaluation of its performance. In order to validate the model, several numerical predictions for different triaxial and oedometric tests of various clays have been compared with results from previous works. 1,11,59,[62][63][64][65] Actually, 71 triaxial tests, one oedometric test and eight incremental response envelopes of five normally consolidated and overconsolidated soft and stiff cohesive soils covering a wide plasticity range are presented. A representation of these materials in Casagrande's plasticity diagram is shown in Figure 4. Unless otherwise specified, in the following figures, the dashed lines represent the experimental data and the solid lines illustrate the simulations with the proposed model. These simulations, descriptions and figures are also presented in the PhD thesis of the first author. 45 Lying below the A-line (I P = 12.2%) (see Figure 4) the Kaolin clay 11 has to be classified as a silt. Vallericca clay 1 is a stiff overconsolidated clay deposit of Plio-Pleistocene age, which lies slightly over the A-line (I P = 31.6%). Thus, it can be considered as a high plasticity clay. The Boom clay 63 is a tertiary clay of Oligocene age, similar to the London clay.   The overconsolidation ratio was estimated to be OCR = 2.6, and following the Casagrande plasticity diagram, it has to be classified as a high plasticity clay (I P = 47%). Beaucaire Marl 65 was deposited in a shallow marine environment during the Pliocene age. It is a medium plasticity clay (I P = 17%). Schwerin mud 59 is a clay with organic additives and lies outside the common Casagrande diagram (I P = 114%); see Figure 4B. The tests were performed on reconstituted samples of Kaolin, Beaucaire Marl and Schwerin mud and both reconstituted and natural samples of Boom clay and Vallericca clay. The calibrated material parameters are listed in Table 1. Its calibration process is explained in detail in Appendix B1.

Soft and/or normally consolidated clays
In this section, the predictions of the proposed model are compared with triaxial tests carried out on reconstituted normally consolidated Kaolin samples 11 with variation of the initial mean pressure p 0 and the cutting direction of the sample. The displacement rate was held constant to . s = 0.025 mm/min corresponding to an axial strain rate of . 1 = 0.05%/min. Figure 5 represents five undrained monotonic triaxial tests and their numerical simulations. The initial mean pressure varied between p 0 = 50, 100, 200, 300 and 400. In the case of p 0 = 50 kPa, the sample is slightly overconsolidated with the initial overconsolidation ratio OCR=1.33. All other samples are sheared from a normally consolidated state. Following these indications given in Wichtmann and Triantafyllidis, 11 the void ratio was initialized according to e 0 = e i0 − ln(OCR p 0 ).
The effective stress path inclination ( Figure 5A) indicates the existence of inherent anisotropy. This slope was possible to be captured by the proposed model by introducing Equation (25); thus, the anisotropic coefficient takes the value = 2. Experiments carried out on Brno clay 56 show similar results with comparable inherent anisotropy properties. Figure 5B,C show the stress-strain relationship and the pore water pressure against the axial strain, respectively. The model was calibrated to obtain the same peak magnitude q max , which appears at a smaller strain compared with the experimental result. This delayed peak can be attributed to the experimental technique, because at a strain of 15%, one cannot assure a homogeneous strain distribution. 1,66 On the other side, this could be an indication that the degree of nonlinearity Y converges too fast to the limit value Y → 1. All other results of the numerical calculations describe very well the experimental observations.
In Figure 6, the results and predictions of three undrained triaxial tests with a common initial pressure of p 0 = 100 kPa and different cutting directions are shown. The results of the tests demonstrate once more the effect of the inherent  Figure 6A compared with the results illustrated in Figure 5A, also cut out in vertical direction. The effective stress path of the horizontal sample is inclined to the right showing a more dilative response and a higher undrained shear strength. Also, Duncan and Seed 9 documented that the undrained shear strength s u may either increase or decrease with c , depending on the type of clay. The tests are well predicted by the proposed model and the inherent anisotropy is properly modeled. The prevented dilatancy of the effective stress path (increase of deviatoric stress) of the horizontal sample, which can be observed in Figure 6A, results in pore water relaxation; see Figure 6C.

Stiff and/or overconsolidated clays
This section represents simulations of experiments carried out on structured and unstructured stiff clays with different overconsolidation ratios. The samples were first subjected to drained isotropic preloading towards p max = 200,400 or 800 kPa. In order to achieve different overconsolidation ratios, a drained unloading followed to p 0 = 100 kPa. The displacement rate was kept constant during shearing . s = 0.025 mm/min corresponding to a strain rate of . 1 = 0.05%/min. The void ratio was initialized according to e 0 = e i0 − ln(OCR p 0 ).
With increasing OCR, the material response is encountered more dilative and higher undrained shear strengths are reached. Similar effects of OCR on dilatancy or on the undrained monotonic strength are also documented in previous works. 1,8,[68][69][70][71] The models simulations of the effective stress path ( Figure 7A) show an early ill-matched contractancy for OCR >1. In this sense, it is worth noting that, owing to the foregoing isotropic compression, these samples had experienced significant modifications of their structure prior to shearing. Thus, the bounding and fabric are likely to have been altered along such loading paths, and the anisotropic factor should evolve with the plastic strain. Other authors, for example, Masin,17 documented that the change in anisotropy during isotropic loading is negligible. In this case, a suitable relation describing the evolution equation of is proposed to read .
However, the development of the stiffness anisotropy remains to be clarified in future research. After reaching the failure surface, the models response is dilative and corresponds to the experimental behaviour ( Figure 7A). The undrained shear strength, the stress-strain relationships and the excess pore water pressure are well modeled; compare also Figure 7B and Figure 7C. The softening in Figure 7B may result as an effect of discontinuous slipping as reported in Amorosi and Rampello, 1 indicated with an asterisk (*) and dashed line in Figure 7B.
The failure surface shows a curvature instead of a straight line, as introduced in some constitutive models, in the p−q space, similar as reported in previous studies. 11,72,73 A closer look to the phase transformation line (PTL) is presented in Figure 7D, whereby the points on the PTL are plotted in the PTL versus OCR diagram. OCR denotes hereby the initial overconsolidation ratio. Whereas the model does not show any dependency of PTL on the initial OCR, the experimental data for Kaolin show a strong dependency up to OCR ≈ 5. It is indicated that after this, an asymptotic state is reached.
In Figure 8A, state paths from undrained triaxial compression and extension tests with effective pressures up to about 4 MPa on intact and reconstituted samples of Boom clay are shown. 63 Figure 8B shows the simulations of these tests. Hereby, we restrict our attention to the shape of stress paths normalized by the equivalent pressure at the isotropic normal compression line p ei = exp((e i0 − e)∕ ) (for detailed definition, see also Figures 1E,F and 2). The peak states for the intact samples exceed those for the reconstituted samples, which can be attributed to the inherent anisotropy. It can be observed that the proposed model predicts well dilatant/contractant behaviour for a wide range of overconsolidation ratios in compression as well as in extension and the increase of the peak shear strength for overconsolidated states is also well captured by the model. Note that all material constants for natural as well as for reconstituted samples resulted to same values; see Table 1, except the anisotropic coefficient , which results in = 1.5 for intact samples and = 1 for reconstituted ones. The difference on the inherent anisotropy between reconstituted and natural samples has been documented also for other clays of marine or lacustrine origin by previous studies. 74,75 In the following figures, some high-and medium-pressure (HP and MP, respectively) tests of isotropically and anisotropically consolidated (CI and CA, respectively) samples of intact Vallericca clay reported by Amorosi and Rampello 1 are compared with predictions of the proposed model. According to the interpretation criteria specified in Amorosi and Rampello, 1 stress-strain data in Figures 9A, 10A and 11A are shown by full lines for the portion of test results for which homogeneous straining was inferred and dotted lines for which discontinuous slipping was assumed. Figure 9A shows three undrained triaxial tests of natural samples prior isotropic compression to high effective stresses, and Figure 9B shows the simulations with the proposed model. Taking a look into the parameters listed in Table 1, = 1 indicates no inherent anisotropy for this intact material. In fact, both bonding and fabric are likely to be altered along extended compression paths these samples were exposed to. The normally consolidated samples show similar contractant behaviour, whereby the stress ratio approaches the residual state = q/p ′ = 1 at large strains. The model captures very well all the aforementioned effects of the experiments. The overconsolidated sample shows a dilatant behaviour with a reduction of the pore water pressure, and it reaches a peak stress ratio at 1 ≈7 %. The models response captures these observations well, yet reaches the peak stress ratio already at 1 ≈ 3.8%.  Figure 10A shows the experimental results for drained and undrained shearing of natural samples. Prior to shearing, the natural samples were anisotropically recompressed to achieve a normally consolidated state, and when necessary, anisotropically unloaded under K 0 conditions. 1 Shearing of overconsolidated samples shows that the undrained shear strength decreases, whereby the value of axial strain at which the peak is attained increases with the OCR. The models response, presented in Figure 10B, shows both an increasing value of axial strain at which the peak is attained as well as an increasing undrained shear strength for increasing OCR > 1, which was observed also by other experimental works. 8,11,[68][69][70][71] Shearing of normally consolidated samples shows a continuous contractant behaviour, whereas the increase of both the volumetric strain and the pore water pressure with the axial strain is evident. The residual or the critical state corresponding to = q ′ /p = 1 and stationary values of Δu and v are reached for both the overconsolidated and the normally consolidated samples.
The tests and simulations presented in Figure 11 were carried out on samples anisotropically compressed to achieve vertical effective stress ′ v,max ≈ 6750 kPa ≈ 2.5 ′ v , with ′ v being the vertical yield stress observed in oedometric tests. 1 Shearing of overconsolidated samples shows a similar behaviour to that of the CI-HP tests, but here, the undrained shear strength as well as the value of axial strain at which the peak is attained increases with the OCR. A peak is evident for OCR >3. All experimental observations where well captured by the proposed model; see Figure 11B. Furthermore, the models' response shows a maximum of the pore water pressure and a decrease thereafter to a constant value. The experimental data show a subsequent increase of Δu for lower OCRs, whereby at this stage, the discontinuities are well developed (dotted lines in Figure 11A). The normally consolidated samples show a contractant behaviour, which is well reproduced by the model.
Owing to the high level of imposed compression, the HP samples, both isotropic and anisotropic consolidated (CI-HP and CA-HP), are expected to have experienced significant modifications of their structure prior to shearing. Actually, the inherent anisotropy is expected to be altered and the behaviour of these samples to render the same as reconstituted samples. On the other side, less structure change prior to shearing is expected for the medium-pressure anisotropically consolidated samples (CA-MP). To investigate these effects, Figure 12A shows the effective stress paths of CA-MP and CA-HP test series and Figure 12B presents their simulations with the proposed model. It can be observed that the inclination of the effective stress path of normal and slightly overconsolidated samples on one side and overconsolidated samples on the other of both CA-MP and CA-HP test series does not differ, although the peak effective stresses of the two test series are substantially different. The peak of normal and slightly overconsolidated samples is reached after small changes in mean stress p ′ , and then true contractant behaviour is exhibited. The same effects render the simulations with the proposed model. Overconsolidated samples show dilatant behaviour, similar to that observed for other clays. 11,76 The stress paths reach the peak at the failure envelope, which is in good agreement with the simulations.

Time-dependent behaviour of clays
The higher the strain-rate dependency of clays, the higher the plasticity of the material. This effect is documented also in other works. 7,10,11,[77][78][79][80][81] In this section, we will restrict our attention to the strain-rate dependency of the Kaolin clay revised in Wichtmann and Triantafyllidis 11 and Schwerin mud reported in Krieg et al. 59 The calibrated parameters for these materials are listed in Table 1, whereby the Kaolin samples presented in Figure 13 are cut out in vertical direction (parallel to the sedimentation direction). The experiments 11 (dashed lines) and simulations (solid lines) of four undrained triaxial tests with variation of the axial strain-rate Figure 13A in the effective stress space, Figure 13B in form of stress-strain relations and Figure 13C in form of excess pore water pressure-axial strain relations. As expected, both experiments and simulations show an increasing undrained monotonic strength and decreasing prevented contractancy with increasing strain rate. These effects are stronger pronounced for materials with higher plasticity. Consequently, the higher the strain rate, the higher the excess pore water pressure; see Figure 13C. The simulations are in well agreement with the experimental behaviour.
Due to its higher plasticity, the strain-rate dependence is more pronounced for the Schwerin mud tested by Krieg et al. 59 A CU triaxial test carried out with strain rate jumps following the deviatoric strain rates  Figure 14A. The effective stress paths described by the experiment and by the response of the proposed model are practically identical. Note that the relation between the compression and swelling index for Schwerin mud renders / = 15; see Table 1. Other models are likely to have problems dealing with / 8 as pointed out in Masin. 17 For large ratios of / , that is, stiff response at unloading, the swept-out-memory (SOM) surface becomes nonconvex in the vicinity of isotropic stress states.  Figure 14B presents an oedometric test and its numerical prediction. The experiment consists of several constant-rate-of-strain steps with different rates and unloading-reloading cycles The experimental path includes also three creep and one relaxation steps. The evolution of the void ratio and stress and the strain-rate dependencies are very well reproduced by the model.

Directional response
In the following, the response of reconstituted Beaucaire Marl to the stress probing program detailed in Constanzo et al. 65 is simulated with the proposed model. The predictions of the model are described using the so-called incremental response envelopes (IRE). 82 A representation with response envelopes (RE) was first introduced by Gudehus 83 as a convenient tool for visualizing the tangential stiffness and other properties of nonlinear and rate-type constitutive equations. Roughly speaking, a response envelope is a polar diagram of stiffness plotted for different directions of stretching. In the general case, a response envelope is a hypersurface in a six-dimensional space. However, for the particular case of triaxial loading, the RE can be conveniently represented in the isometric Rendulic plane (Δ 1 − √ 2Δ 3 or Δ 1 − √ 2Δ 3 ). We start by choosing an initial stress 0 . The strain RE is obtained as a plot of the final strains calculated with normalized stress probes diag[ 1 , 3 , 3 ] (with || . ||Δt = const. or with || . || = 1) applied in different directions. A constitutive model may be understood as a  Test no.  35  0  Tx128  Tx115  46  21.91  -Tx130  90  71.57  Tx121  Tx132  126  90  Tx126  Tx119  154  104.49  -Tx116  180  123.69  Tx123  -215  180  Tx127  Tx134  226  201.91  -Tx129  270  251.57  Tx122  Tx117  305  270  Tx125  Tx113  0  303.69  Tx124  Tx118 a mapping that carries a circle plotted in the Δ 1 − √ 2Δ 3 space to the strain Δ 1 − √ 2Δ 3 space where it becomes an ellipse. The size of each strain increment vector defining the RE can be directly interpreted as a directional secant compliance of the material, for the respective loading direction and stress magnitude. By simply replacing rates with finite increments, the same definition stands for the IRE.
The testing program reported in Constanzo et al. 65 and simulated with the proposed model consists of 20 drained stress probes, starting from a common initial stress state (either A with q/p ′ = 0kPa/150kPa=0 or B q/p ′ = 60kPa/150kPa=0.4) and pointing in different directions in the triaxial plane; for details see Table 2 and Figure 15. The axial stress rate was maintained equal to . 1 = 2 kPa/h. Both initial states represent virgin states for the material, that is, OCR = 1. For a representation in form of IREs with finite-size increments, the interpretation of elastic and (visco)plastic properties of the model is necessary. The size of each strain increment vector defining the IRE can be directly interpreted as a directional secant compliance of the material for the respective loading direction and stress increment magnitude as depicted in Constanzo et al. 65 Figures 16 and 17 show the experimental and numerical IREs at R = 20,30,40 and 50 kPa for initial state A and B, respectively. The IREs show that the softest response of the material is reached for those paths characterized by a large deviatoric component, that is, Tx125, Tx122, Tx126 and Tx121 for initial state A and Tx119 and Tx113 for initial state B. The stiffest reponse of the material is rendered upon full stress path reversal with respect to the previous history, that is, Tx123 and Tx127 for initial state A and Tx134 and Tx129 for initial state B. This behaviour of soils is documented in many other works. 30,59,63,73 As expected, for isotropic loading, the response of the material is softer when the stress path points in the direction of continued loading, that is, Tx124 for initial state A and Tx130 for initial state B. It is interesting that due to the stiffness increase the IREs become nonconvex for isotropic loading corresponding to Tx128 (initial state A) and Tx115 (initial state B). An evidence of incremental nonlinearity can be observed due to the nonsymmetry of the IREs for both initial states A and B about the origin of the strain increment space. Figures 16 and 17 suggest that this effect is evident already at strains about 0.15% and stresses about R = 20 kPa. Besides the salient features discussed above, worthy to note is the progressively shift of the IREs downwards to the right for the initial state A and upwards to the left for initial state B (Figures 16 and 17, respectively). Arising with this behaviour, the IREs become increasingly nonsymmetric about the isotropic axis the higher R is. For initial state A, the secant compliance increases for those directions close to negative deviatoric direction, that is, axisymmetric extension (Tx125). As initial state B is closer to the critical state line for axisymmetric compression than extension, the secant compliance increases for those directions close to positive deviatoric direction, that is, axisymetric compression (Tx119). For the significant increase of the secant compliance (probes Tx125 for initial state A and Tx119 for initial state B), the question arises if the magnitude of the increase of the secant compliance is reliable considering that such an observation is based on the results of one single probe.
The predictions of the proposed model ( Figures 16B and 17B, black lines) are from both quantitative and qualitative standpoints in a perfect agreement with the features discussed above. The decrease of the secant stiffness for Tx125 (initial state A) and Tx119 (initial state B) is slightly underpredicted. All other characteristics, considering also the nonsymmetry of the IREs around the origin of the strain increment space, nonconvexity about the isotropic stress path and stiffness increase upon stress path reversal are well captured by the numerical simulations. Other models such as MCC, 3-SKH, CLoE, K-hypoplastic model with and without intergranular strain reviewed in Masin et al. 84 (only for the probe with initial state B) appear to significantly underpredict the stiffness of the material. This presents a fundamental shortcoming for the performance-based design of geotechnical engineering problems. From a qualitative standpoint, the REs of the two K-hypoplastic models and (to a much lesser extent) those of CLoE show some degree of nonconvexity in a region located around = 0.4 (i.e., upon continued loading with respect to the consolidation history). 17 The IREs predicted with the two elastoplastic models (MCC and 3-SKH) show a convex shape with some irregularities of the MCC envelope close to neutral loading in extension. For comparison purposes, in Figure 17B, also the IREs obtained with Masins model reported in Masin 84 (K-hypoplastic model with intergranular strain) are illustrated with dashed blue lines. The IREs predicted with the proposed model, however, appear from a qualitative and a quantitative point of view, all in a very good agreement with the features of the experimental observations discussed above.

Limitations of the model
The definition of the flow rule, Equation (31), imposes a consistency of K 0 with Jaky's formula. 60 It is, however, not tested in detail if the flow rule fits the experiments by Pradhan et al. 85 nor the dilatancy theory of Rowe et al., 86 inasmuch as they present works based on sands. The state-dependent dilatancy for clays proposed in Tafili and Triantafyllidis 4 is surely not predicted well with the model; see Figure 7. Improvements in this direction are left out for future research.
The behaviour of clays exposed to cyclic loading remains also to be developed in future work. Basic principles of this model are, however, given in previous works 28,37 based on the ISA plasticity.
It is well known that hypoelastic constitutive relations cannot be expressed in terms of an elastic strain energy function. The response is path dependent, and dissipation may occur even though the material is supposed to be elastic. However, the main point is that the elastic strains are assumed to be relatively small so that any error in the conservation of energy is very small and thus can be neglected. The stress rate . must be objective; thus, the material derivative of the Cauchy stress cannot be used. But other objective rates such as the Jaumann, Green-Naghdi and Truesdell can be used instead. Usually in hypoplastic constitutive equations, as well as for the proposed model, the Jaumann rate is used. To overcome these shortcomings and in order to obtain a conservative constitutive model for clays under purely elastic conditions, a hyperelastic stress-strain relation derived from an elastic potential is being developed by the authors. Clays are geomaterials ubiquitous in sedimentary soils and rocks and are therefore involved in a wide area of geotechnical applications. 34,87 Many of them include thermal loadings, for example, geothermal structures where heat exchanger piles are used or nuclear waste disposal. When thermal energy is stored in clays, the material behaviour will change as a result of changing temperature. Early studies [88][89][90][91] and recent ones 59,92,93 have documented that the compression index was found independent of temperature, yet the higher the temperature, the lower the void ratio at any given consolidation pressure: Δe = e 0 − e = C T ln(T∕T 0 ), whereby C T would represent a new material parameter, namely, the temperature coefficient. For a constant displacement rate . = const, the studies have shown that under different temperatures, unique lines with the same slope on the e−ln( ′ ), namely, isotherms, are reached. The void ratio difference between two isotherms corresponds to Equation (34) and a formerly heated soil behaves like an overconsolidated one until the corresponding isotherm is reached. Some other models, as for example, the one proposed by Masin and Khalili, 20 introduce explicitly a temperature-dependent compression index (T) as well as a temperature-dependent characteristic void ratio e 0 (T).
Interesting to note is that the rate dependency of the soil is found independent of the temperature I v = C / = const(T). 59,90,91,94 Consequently, the temperature and rate dependence of the soil can be superposed in the model formulation: with . temp the strain rate accounting for temperature dependence of the soil. According to the work of Lovisa et al., 92 the total increase in shear strength was about 40% at high temperature. Their results show also that the preconsolidation pressure follows a similar development compared with the shear strength. The strengthening effect appears for most part to depend on changes of the microstructure in the clay, which influences also the fabric. These results and Equations (34) and (35) present a basis for an enhancement of the model to account for temperature dependency. Thus, the model provides a platform where different effects can be integrated without changing the structure of the mathematical formulation.

FE-CALCULATION OF A SHALLOW FOUNDATION
In this section, the proposed model and its implementation within the commercial FE program Abaqus is validated by a numerical analysis of a shallow foundation as presented also in Tafili. 45 The width of the foundation is 1 m, and the size of the considered soil region is 20 × 10 m; see Figure 18A. Accounting for the symmetry with respect to the vertical axis, only half of the system has been modeled. A vertical displacement of 0.2 m is prescribed at the nodes under the foundation. The period of time for the punching was varied between t = {1, 10, 100, 1000, 10000, 100000} s, which renders u 2 = u 2 ∕t = 2 × {10 −1 , 10 −2 , 10 −3 , 10 −4 , 10 −5 , 10 −6 } m/s allowing to evaluate the rate dependence of the soil. A clay with low plasticity and a clay a high plasticity clay were used for the calculations, which correspond to Kaolin clay and Schwerin mud calibrated in Section 5, respectively. The parameters are listed in Table 1. The contact between soil and foundation is perfectly rough, that is, the horizontal movement of the nodes under the foundation is prohibited.
The FE mesh containing ≈12 000 elements is illustrated in Figure 18B. In order to prevent volumetric locking, elements with reduced integration (CPE4R) have been used. Although a displacement of 0.2 m was prescribed, the geometric nonlinearity was switched off in order to enable fast calculations for comparison purposes.
The simulations contain undrained conditions with variation of the displacement rate . u 2 and the anisotropic coefficient . In order to simulate the undrained condition, the bulk modulus of water K w = 2×10 6 MPa is multiplied by the volume strain increment to update the pore water pressure p w . The stiffness of water is considered for the calculation of the material Jacobian J = J mat + J w = J mat + K w 1 ⊗ 1. The initial stress conditions were adjusted to the vertical stress 1 = ′ h and the horizontal stress 2 = K 0 1 = (1 − sin( c )) 1 . The initial void ratio corresponds to OCR = 1 and has been calculated through the relation e 0 = e i0 − ln(OCR p 0 ). Figure 19 shows the load-settlement curves, whereby in (A) for Kaolin clay and (C) for Schwerin mud the displacement rate has been varied between . u 2 = u 2 ∕t = 2 × {10 −1 , 10 −2 , 10 −3 , 10 −4 , 10 −5 , 10 −6 } m/s and in (B) the anisotropic coefficient for Kaolin has been varied = {0.7, 1.0, 1.5, 2.0}. In each case, a bearing capacity was reached, whereby the displacement at the limit strength varied from 0.08 to 0.18 m. As expected, with increasing displacement rate, higher asymptotic reaction forces are reached. Whereas, due to the low plasticity of Kaolin, this effect is rather weak pronounced in the calculations in Figure 19A, calculations accounting for the high plasticity Schwerin mud show a much more pronounced difference between the bearing capacities at different displacement rates; see Figure 19C. Fur-thermore, the higher plasticity leads to an increasing displacement at which the bearing capacity is reached, which is of great importance for performance-based designs. It is worthy to note that the limit undrained strength corresponds for each calculation approximately to the limit analysis solution developed by Prandtl 1 = (2+ ) c u . Figure 19B indicates no dependence of asymptotic reaction force on the anisotropic coefficient. It is, however, interesting to observe that the displacement at which the asymptotic state is reached is substantially affected by the fabric of the soil. On that note, u 2 takes values from u 2 = {0.04, 0.05, 0.08, 0.11} m for = {0.7, 1.0, 1.5, 2.0}, respectively. Therefore, a 2.2 times larger displacement is necessary to reach the asymptotic reaction force for = 2 compared with an isotropic material response = 1 (on the example of Kaolin). This indicates that for a performance-based design, not only the viscosity represents an important feature but also the simulation of the fabric of the soil is important.

CONCLUDING REMARKS
A simple hypoplastic model was formulated to describe the fabric effects and rate dependency of clays, without imposing the limitation to only viscous clays. Thus, also time-independent behaviour of clays can be described with the same model. For this purpose, the strain rate is composed of the total strain rate, the hypoplastic strain rate (for instantaneous plastic strain) and the viscoplastic strain rate (accounting for viscous effects). The concept of the bounding surface in conjunction with the maximum void ratio e i , as well as the critical state surface in conjunction with e c , is introduced in the model. Following this, the loading surface and the overconsolidation ratio are reviewed and reformulated according to the maximum void ratio e i , which defines the isotropic normally consolidated state (OCR = 1). The anisotropic parameter is introduced to account for fabric by rotating the (hypo)elastic stiffness tensor to a cross transversal one. An evolution equation for with the hypoplastic strain rate is proposed as well.
Qualitative and quantitative evaluation as well as experimental validation have demonstrated the good performance of the proposed model in describing the mechanical behaviour, particularly the time and strain rate effects as well as the fabric of natural and reconstituted clays. No restriction with respect to the value of the initial overconsolidation ratio has been found. For each tested material, the same material parameter set has been used for different experimental tests witnessing the prediction accuracy of the proposed model. Finally, the model has been used to predict the bearing capacity of a shallow foundation with Abaqus. As expected, the undrained shear strength turned out to depend on the displacement rate. The simulations show furthermore that for a performance-based design, not only the viscosity represents an important feature but also the simulation of the fabric of the soil is important.
Some other experimental observations, such as the coupling with temperature effects and an accurate description of dilatancy, represent the focus of future research.

APPENDIX A: CALIBRATION OF MATERIAL PARAMETERS
The model requires the calibration of eight material parameters: five of them ( , , e i0 , h and M c ) are common parameters for clayey soils, for example, typical for the MCC theory as well, one parameter for the viscosity I v , one parameter for the fabric and the material constant f b0 controlling the maximum stress ratio for triaxial compression. They are listed in Table A1, whereby a subdivision into different groups according to their role within the model has been done and the approximative range for their values as well as the required tests for their calibration are noted. The following lines are devoted to explain the calibration process.
A compression path with constant strain rate is expected to approach asymptotically to a line with slope equal to in the e versus ln(p) space. The parameter is calibrated instead upon an unloading path.
The reference void ratio e i0 corresponds to the maximum void ratio at the reference pressure p = 1 kPa evaluated at the virgin consolidation line OCR = 1. According to the proposed model, the maximum void ratio curve is reached under infinite fast strain velocity || . || = ∞. Due to the experimental complexity when dealing with infinite strain velocity, an alternative calibration method for this parameter has been developed by Fuentes et al. 37 For instance, we consider an experimentally feasible velocity D r at which OCR = OCR r ≠ 1 holds. This curve is named reference isotach and D r reference creep rate. 30,37 After some simple mathematical manipulations of the constitutive relation (2), the following equation is obtained: (A1)  Conceiving the numerical solution for OCR r , the parameter e i0 can be computed e i0 = e r − ln (OCR r ) , whereby e r is extrapolated from the experimental curve with || . || = D r at p = 1 kPa.
The slope of the triaxial compression critical state line M c is calibrated connecting the endpoints of the available triaxial compression tests in the p−q plane. The bounding surface factor f b0 controls the compression triaxial stress ratio = q/p at higher overconsolidated states OCR > 2. When data are scarce, however, a value of f b0 = 1.3 is recommended according to the authors' experience.
The parameter describing the intensity of creep for normally consolidated samples, denoted as the viscosity index I v , has been examined also in former works. 30,95,96 One can show that considering two different isotachs with || . a || and || . b || corresponding to OCR a and OCR b , respectively, the following relation 30,97 holds: . (A3) The calibration of the parameters related to fabric, namely, and h , can be achieved very simply using the following self-explained code written in the commercial numerical analysis software Mathematica; see Figure A1. 30 Hereby, the effective stress inclination Δq/Δp can be measured from either monotonic or cyclic undrained triaxial tests and gives the ordinate axis value in Figure A1. After variation of the coefficient of anisotropy , one can read the appropriate Poisson ration h on the abszissa.