# Comparison of Anisotropic Rate-Dependent Models for Modeling Consolidation of Soft Clays

Minna Karstunen, Mohammad Rezania, Nallathamby Sivasithamparam, Zhenyu Yin

► To cite this version:

Minna Karstunen, Mohammad Rezania, Nallathamby Sivasithamparam, Zhenyu Yin. Comparison of Anisotropic Rate-Dependent Models for Modeling Consolidation of Soft Clays. International Journal of Geomechanics, 2014, 15 (5), pp.115-119. 10.1061/(ASCE)GM.1943-5622.0000267. hal-01005979

HAL Id: hal-01005979

https://hal.science/hal-01005979

Submitted on 13 Jun 2014

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Comparison of Anisotropic Rate-Dependent Models for Modeling Consolidation of Soft Clays

M. Karstunen<fnref n="1" />; M. Rezania<fnref n="2" />; N. Sivasithamparam<fnref n="3" />; and Z.-Y. Yin<fnref n="4" />

**Abstract:** Two recently proposed anisotropic rate-dependent models are used to simulate the consolidation behavior of two soft natural clays: Murro clay and Haarajoki clay. The rate-dependent constitutive models include the EVP-SCLAY1 model and the anisotropic creep model (ACM). The two models are identical in the way the initial anisotropy and the evolution of anisotropy are simulated, but differ in the way the rate effects are taken into consideration. The models are compared first at the element level against laboratory data and then at the boundary value level against measured field data from instrumented embankments on Murro and Haarajoki clays. The numerical simulations suggest that at the element level, the EVP-SCLAY1 model is able to give a better representation of the clay response under oedometric loading than ACM, when the input parameters are defined objectively. However, at the boundary value level, the issue is not as straightforward, and the appropriateness of the constitutive model may depend heavily on the in situ overconsolidation ratio (OCR).

**keywords:** Rate-dependency; Soil anisotropy; Embankment on soft soil; Finite-element method

## Introduction

Creep and rate dependence have fascinated geotechnical engineers since the early 1920s, but it is only relatively recently that comprehensive rate-dependent constitutive models have tried to account for some fundamental features of natural soil behavior, such as anisotropy and the effect of apparent bonding (see Zhou et al. 2005; Leoni et al. 2008; Hinchberger and Qu 2009; Karstunen and Yin 2010; Yin et al. 2011). The aim of this paper is to compare objectively the predictions by two anisotropic rate-dependent models: the EVP-SCLAY1 model (Karstunen and Yin 2010), which is an overstress model (Perzyna 1963), and the anisotropic creep model (ACM) proposed by Leoni et al. (2008). First, the models are compared at the element level against laboratory data. The comparisons are then extended to the boundary value level, considering two instrumented

embankments on soft clay, namely Murro test embankment (Koskinen et al. 2002) and Haarajoki embankment (Finnish National Road Administration 1997) in Finland. For the sake of simplicity, the effects of apparent bonding and destructuration of the rate-dependent behavior are ignored; although, as demonstrated by several researchers, e.g., Hinchberger and Qu (2009) and Yin et al. (2010), they are necessary for modeling certain phenomena, such as strain softening, tertiary creep, and creep rupture.

## ACM and EVP-SCLAY1 Models

The two models, ACM and EVP-SCLAY1, are identical in the way the initial anisotropy and the evolution of anisotropy are modeled, but differ in the way the rate effects are accounted for. The complete formulations can be found in Appendix. The initial anisotropy is modeled by defining anisotropic surfaces, attributed to the history of the soil deposit, which define the boundary between large irrecoverable strains and relatively small strains. These surfaces in stress space are analogous to the anisotropic yield surface in the rate-independent S-CLAY1 model (Wheeler et al. 2003). As with the S-CLAY1 model, it is assumed that because of subsequent loading that produces irrecoverable strains, there is an evolution in the orientations of these surfaces as a function of irrecoverable volumetric and shear strain increments. As demonstrated by Karstunen and Koskinen (2008), this assumption offers a simple way to represent extremely well the evolution of anisotropy for reconstituted clays. When combined with a destructuration law, the models can be extended to account for the large strain behavior of natural clays (Karstunen and Koskinen 2004; Karstunen et al. 2005).

In the case of the EVP-SCLAY1 model (Karstunen and Yin 2010), this boundary surface is called a static yield surface [see Fig. 1(a)], and it is assumed that no rate-dependent behavior occurs when the stresses are within the surface. This basically means that there is a threshold value of effective stresses that needs to be exceeded before rate-dependent behavior starts, as discussed by Qu et al. (2010). In contrast, in the ACM model, the boundary between the large irrecoverable strains and the relatively small strains is called normal consolidation surface (NCS) [see Fig. 1(b)], and it is assumed

<footnote n="1">¹Professor, Dept. of Civil and Environmental Engineering, Chalmers Univ. of Technology, SE 412 96 Gothenburg, Sweden; Professor, Dept. of Civil and Environmental Engineering, Univ. of Strathclyde, John Anderson Building, 107 Rottenrow, Glasgow G4 0NG, U.K.; Docent, Aalto Univ., School of Science and Technology, FI-00076 Aalto, Finland. E-mail: minna.karstunen@chalmers.se</footnote>

<footnote n="2">²Lecturer, Nottingham Centre for Geomechanics, Dept. of Civil Engineering, Univ. of Nottingham, Coates Building, University Park, Nottingham NG7 2RD, U.K.; formerly, Research Associate, Dept. of Civil and Environmental Engineering, Univ. of Strathclyde, John Anderson Building, 107 Rottenrow, Glasgow G4 0NG, U.K. (corresponding author). E-mail: mohammad.rezania@nottingham.ac.uk</footnote>

<footnote n="3">³Marie Curie Postdoctoral Fellow, Dept. of Civil and Environmental Engineering, Univ. of Strathclyde, John Anderson Building, 107 Rottenrow, Glasgow G4 0NG, U.K.; Researcher, Plaxis bv, Delft, 2628 XJ, Netherlands. E-mail: ns@plaxis.com</footnote>

<footnote n="4">⁴Professor, Center for Marine Geotechnics Research, Dept. of Civil Engineering, Shanghai Jiao Tong Univ., Shanghai 200240, China. E-mail: zhenyu.yin@gmail.com</footnote>
that creep strains occur even in the overconsolidated range, i.e., when the current stress surface in Fig. 1(b) is smaller than the NCS.

The incremental viscoplastic strains in the EVP-SCLAY1 model (Karstunen and Yin 2010) can be calculated as

$$ \dot{\epsilon}_{ij}^{vp} = \mu \left\langle \exp \left[ N \cdot \left( \frac{p'_{m}^{d}}{p'_{m}^{s}} - 1 \right) \right] - 1 \right\rangle \frac{\partial f_{d}}{\partial \sigma'_{ij}} \quad (1) $$

where $\mu$ = fluidity parameter; $\langle \rangle$ = McCauley brackets; and $N$ = strain-rate coefficient relating to the strain-rate effect on shear strength and preconsolidation pressure. $N$ and $\mu$ are therefore the key parameters controlling the viscoplastic strain rate. $p'_{m}^{d}$ and $p'_{m}^{s}$ = size of the dynamic loading surface and the static yield surface, respectively [see Fig. 1(a)]; and $f_{d}$ = viscoplastic potential function, represented by the dynamic loading surface.

In ACM, the incremental volumetric creep strains are calculated as

$$ \dot{\epsilon}_{vol}^{c} = \frac{\mu^{*}}{\tau} \left( \frac{1}{p'_{p}/p'_{eq}} \right)^{\beta} \quad (2) $$

where $\beta = (\lambda^{*} - \kappa^{*})/\mu^{*}$ = creep exponent; $\lambda^{*}$, $\kappa^{*}$, and $\mu^{*}$ = modified compression, swelling, and creep indices; $\tau$ = reference time (= 1 day for conventional oedometer test); and $p'_{p}$ and $p'_{eq}$ = sizes of the normal consolidation surface and the current stress surface [see Fig. 1(b)]. Therefore, by comparing Eqs. (1) and (2), it becomes clear

that with ACM, creep strains will be predicted even at overconsolidated range, in contrast to EVP-SCLAY1, whereas both models somehow relate the actual creep rate to the relative sizes of the respective surfaces shown in Fig. 1.

The advantage of ACM over EVP-SCLAY1 is that it is possible to derive the required soil constants for ACM directly from the experimental tests, whereas the viscosity coefficients of the EVP-SCLAY1 require calibration, as shown by Karstunen and Yin (2010). Both models have been implemented in the Plaxis 2D and Plaxis 3D 2010 finite-element code (Brinkgreve et al. 2010) as user-defined soil models, hence enabling the use of the models at both element and boundary value levels.

## Determination of Model Parameters

The EVP-SCLAY1 and ACM models involve a number of soil constants and state variables that can be classified in three groups, as follows:

* Parameters that are similar to the modified Cam clay parameters include soil constants $\nu'$ (Poisson's ratio), $M$ (stress ratio at critical state), $\lambda$ (slope of the normal compression line), and $\kappa$ (slope of the swelling/recompression line). In ACM, the modified compression and swelling indexes, where $\lambda^{*} = \lambda/(1+e)$ and $\kappa^{*} = \kappa/(1+e)$, are used. Furthermore, the initial values for two state variables, namely $e_{0}$ (initial void ratio) and $p_{m0}$ (initial size of the yield surface) are required. In the context of finite-element analyses, the initial value of $p_{m0}$ is calculated based on the OCR (vertical

<figure><img src="image_2.png" /><figcaption>Fig. 1. (a) Yield surfaces of the EVP-SCLAY1 model; (b) normal consolidation surface and current stress surface of ACM model</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 2. Yield surfaces of Murro clay at a depth of 4 and 7 m [data after Karstunen and Yin (2010)]</figcaption></figure>
overconsolidation ratio) or POP (preoverburden pressure), the normally consolidated $K_0^{NC}$ value (lateral earth pressure at rest, estimated by Jaky's formula), and the initial vertical effective stress.

* Parameters describing initial anisotropy and its evolution include soil constants $\omega$ (rate of rotation of the surfaces) and $\omega_d$ (relative rate of surface rotation). The latter can be theoretically derived based on M values [see Wheeler et al. (2003) for details].

Furthermore, the initial values of the fabric tensor describing the anisotropy of the fabric (i.e., the arrangement of particles and particle contacts) need to be defined. For deposits that have experienced one-dimensional (1D) consolidation under their self-weight without significant erosion or unloading, it can be assumed that the fabric is initially cross-anisotropic. In this special case, the components of the fabric tensor can be calculated using

<table><caption>**Table 1.** Parameter Values for Murro Clay at a Depth of 3.6 m</caption><thead><tr><th rowspan="2">λ</th><th rowspan="2">κ</th><th rowspan="2">e<sub>0</sub></th><th rowspan="2">POP (kPa)</th><th rowspan="2">ν'</th><th rowspan="2">M</th><th rowspan="2">α<sub>K0</sub></th><th rowspan="2">ω</th><th rowspan="2">ω<sub>d</sub></th><th colspan="2">EVP-SCLAY1</th><th colspan="2">ACM</th></tr><tr><th>N</th><th>μ (day<sup>-1</sup>)</th><th>μ*</th><th>τ (day)</th></tr></thead><tbody><tr><td>0.5</td><td>0.041</td><td>2.44</td><td>1.0</td><td>0.3</td><td>1.65</td><td>0.66</td><td>20</td><td>1.02</td><td>20</td><td>8.64 × 10<sup>-5</sup></td><td>1.92 × 10<sup>-3</sup></td><td>1</td></tr></tbody></table>

<table><caption>**Table 2.** Parameter Values for Haarajoki Clay at a Depth of 6.4 m</caption><thead><tr><th rowspan="2">λ</th><th rowspan="2">κ</th><th rowspan="2">e<sub>0</sub></th><th rowspan="2">POP (kPa)</th><th rowspan="2">ν'</th><th rowspan="2">M</th><th rowspan="2">α<sub>K0</sub></th><th rowspan="2">ω</th><th rowspan="2">ω<sub>d</sub></th><th colspan="2">EVP-SCLAY1</th><th colspan="2">ACM</th></tr><tr><th>N</th><th>μ (day<sup>-1</sup>)</th><th>μ*</th><th>τ (day)</th></tr></thead><tbody><tr><td>0.96</td><td>0.03</td><td>2.6</td><td>32.0</td><td>0.2</td><td>1.43</td><td>0.55</td><td>49</td><td>0.97</td><td>18</td><td>6.00 × 10<sup>-6</sup></td><td>3.47 × 10<sup>-3</sup></td><td>1</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Conventional oedometer test predictions of Murro clay</figcaption></figure>
a scalar value $α_0$, which, like $ω_d$, is theoretically simply a function of $M$ [see Wheeler et al. (2003) for details].

* Parameters describing viscosity require the input of two viscosity coefficients. These are soil constants $N$ (strain rate coefficient) and $μ$ (fluidity coefficient) for EVP-SCLAY1 [see Karstunen and

Yin (2010) for details], and $μ^*$ (modified creep index) and $τ$ (reference time, which is linked to the definition of vertical pre-consolidation stress) for ACM [see Leoni et al. (2008) for details]. Although these models have a relatively large number of parameters, most of them have a clear physical meaning, and hence

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Conventional oedometer test predictions of Haarajoki clay</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Fig. 5. Geometry of Murro test embankment</figcaption></figure>

they can be determined in a relatively straightforward manner (Karstunen and Yin 2010).

## Model Parameters for Murro and Haarajoki Clays

Murro clay (Koskinen et al. 2002) is a 23-m-thick sulfide-rich silty clay deposit located near the town of Seinäjoki in Finland that has been a subject of many studies (see Koskinen et al. 2002; Karstunen and Koskinen 2004; Karstunen et al. 2005; Karstunen and Yin 2010) because of a well-instrumented test embankment, built in 1993, and extensive programs of specialist laboratory testing. The deposit can be roughly divided into two main layers: a 1.6-m-thick overconsolidated dry crust and an underlying thick layer of almost normally consolidated soft clay. The most compressible layers are at depths between 1.6 and 7 m, which is the focus of the element-level tests presented in this paper. Murro clay is highly anisotropic with regard to yielding, as demonstrated by the yield points in the $p'-q$ plane in Fig. 2. Murro clay is moderately sensitive, with sensitivity around 6–14.

Haarajoki in southern Finland is the location of an instrumented embankment forming a noise barrier, which was the subject of an international competition organized by the Finnish National Road Administration (1997). Half of the 100-m-long embankment has been constructed on the natural clay, and the rest of the embankment is on vertical drains. In contrast to Murro clay, Haarajoki clay is lightly overconsolidated, with an estimated vertical preoverburden pressure >30 kPa. The deposit can be characterized by a high degree of anisotropy and some apparent interparticle bonding, with sensitivity varying from almost zero at the surface to more than 50 at a depth of 15 m (Yildiz et al. 2009).

Conventional oedometer test results, combined with the results from drained and/or undrained triaxial tests on Murro and Haarajoki clays, have been used for determination of average parameter values for the EVP-SCLAY1 and ACM models. As there is natural variation between clay samples, which is further magnified by potential sample disturbance, for each clay layer a number of tests of the same type were used to estimate the average values of parameters. A number of publications report the required parameter values (e.g.,

<table><caption>Table 3. Initial Values for the State Parameters for Murro Test Embankment</caption><thead><tr><th>Layer</th><th>Depth (m)</th><th>POP (kPa)</th><th>e<sub>0</sub></th><th>α<sub>K0</sub></th></tr></thead><tbody><tr><td>1–1</td><td>0.0–0.8</td><td>20</td><td>1.4</td><td>0.63</td></tr><tr><td>1–2</td><td>0.8–1.6</td><td>10</td><td>1.4</td><td>0.63</td></tr><tr><td>2</td><td>1.6–3.0</td><td>1</td><td>1.8</td><td>0.63</td></tr><tr><td>3</td><td>3.0–6.7</td><td>1</td><td>2.4</td><td>0.63</td></tr><tr><td>4</td><td>6.7–10.0</td><td>1</td><td>2.1</td><td>0.63</td></tr><tr><td>5</td><td>10.0–15.0</td><td>1</td><td>1.8</td><td>0.63</td></tr><tr><td>6</td><td>15.0–23.0</td><td>1</td><td>1.5</td><td>0.63</td></tr></tbody></table>

Karstunen et al. 2005; Yildiz et al. 2009; Karstunen and Yin 2010). However, the procedure for calibrating model parameters is presented briefly as follows:

* $\lambda$ and $\kappa$ are measured from conventional 1-day loading oedometer tests on intact samples (from which the corresponding $\lambda^*$ and $\kappa^*$ values for the ACM model can be calculated). For natural clays, $\lambda$ and $\lambda^*$ are not soil constants, and usually the highest values derived from oedometer tests should be used. $M$ is derived from either drained or consolidated undrained triaxial compression test results. The value for Poisson's ratio for each clay is selected according to previous experience (e.g., Karstunen et al. 2005; Karstunen and Yin 2010; Yin et al. 2011), as the triaxial shear tests did not involve unloading-reloading stages that would enable direct determination of purely elastic $v'$ values.

* The anisotropy parameters $\alpha_{K0}$ and $\omega_d$ can be calculated from $M$ values, and the values for $\omega$ can be obtained either by simulating specialist triaxial tests (see Wheeler et al. 2003) or by using the alternative semiempirical formulations suggested by Zentar et al. (2002), Leoni et al. (2008), and Yin et al. (2011). As discussed by Karstunen and Yin (2010), the simulations are rarely sensitive to the value of $\omega$ at the boundary value level.

* By using the value of $C_{\alpha e}$ (creep index) measured from conventional 1-day loading oedometer tests on intact clay samples, the values for $\mu^*$ for the ACM model can be calculated (Leoni et al. 2008). This may sound straightforward, but the interpretation can be difficult in the case of highly sensitive clays, when the creep index is not a soil constant (see Yin et al. 2011).
<table><caption>**Table 4.** Values for the Conventional Soil Parameters and Additional Soil Constants for Murro Test Embankment</caption><thead><tr><th rowspan="2">Layer</th><th rowspan="2">γ (kN/m³)</th><th rowspan="2">ν'</th><th rowspan="2">M</th><th rowspan="2">k<sub>x</sub> (m/day)</th><th rowspan="2">k<sub>y</sub> (m/day)</th><th rowspan="2">ω</th><th rowspan="2">ω<sub>d</sub></th><th colspan="2">EVP-SCLAY1</th><th colspan="2">ACM</th></tr><tr><th>λ</th><th>κ</th><th>λ*</th><th>κ*</th></tr></thead><tbody><tr><td>1–1</td><td>15.8</td><td>0.35</td><td>1.6</td><td>2.13 × 10<sup>-4</sup></td><td>1.64 × 10<sup>-4</sup></td><td>45</td><td>1.02</td><td>0.16</td><td>0.010</td><td>0.0667</td><td>4.20 × 10<sup>-3</sup></td></tr><tr><td>1–2</td><td>15.8</td><td>0.35</td><td>1.6</td><td>2.13 × 10<sup>-4</sup></td><td>1.64 × 10<sup>-4</sup></td><td>45</td><td>1.02</td><td>0.16</td><td>0.010</td><td>0.0667</td><td>4.20 × 10<sup>-3</sup></td></tr><tr><td>2</td><td>15.5</td><td>0.35</td><td>1.6</td><td>2.13 × 10<sup>-4</sup></td><td>1.64 × 10<sup>-4</sup></td><td>25</td><td>1.02</td><td>0.50</td><td>0.030</td><td>0.1786</td><td>1.07 × 10<sup>-2</sup></td></tr><tr><td>3</td><td>14.9</td><td>0.10</td><td>1.6</td><td>1.78 × 10<sup>-4</sup></td><td>1.34 × 10<sup>-4</sup></td><td>20</td><td>1.02</td><td>0.50</td><td>0.036</td><td>0.1471</td><td>1.06 × 10<sup>-2</sup></td></tr><tr><td>4</td><td>15.1</td><td>0.15</td><td>1.6</td><td>1.10 × 10<sup>-4</sup></td><td>9.07 × 10<sup>-5</sup></td><td>25</td><td>1.02</td><td>0.36</td><td>0.030</td><td>0.1161</td><td>9.70 × 10<sup>-3</sup></td></tr><tr><td>5</td><td>15.5</td><td>0.15</td><td>1.6</td><td>6.85 × 10<sup>-5</sup></td><td>5.48 × 10<sup>-5</sup></td><td>25</td><td>1.02</td><td>0.32</td><td>0.034</td><td>0.1143</td><td>1.21 × 10<sup>-2</sup></td></tr><tr><td>6</td><td>15.9</td><td>0.15</td><td>1.6</td><td>1.04 × 10<sup>-4</sup></td><td>8.22 × 10<sup>-5</sup></td><td>30</td><td>1.02</td><td>0.14</td><td>0.004</td><td>0.056</td><td>1.60 × 10<sup>-3</sup></td></tr></tbody></table>

* The values for N and μ for the EVP-SCLAY1 model cannot be determined directly from the experimental data. They can be calibrated by simulation of oedometer test results (see Karstunen and Yin 2010), and by their nature they are not necessarily unique. However, based on the authors' experience, the non-uniqueness is not an issue when boundary value problems such as embankments on soft clay are considered, as the predicted results with different, but yet appropriate, parameter combinations for a given clay are the same from a practical point of view.

* The permeability, k₀, is obtained from conventional 1-day loading oedometer tests on intact samples and is assumed to be the same in vertical and horizontal directions.

The selected values of parameters for Murro and Haarajoki clays used for oedometer test simulations are summarized in Tables 1 and 2. It is worth emphasizing that the values represent objective average values based on a whole suite of tests available at a given depth range, and hence have not been calibrated on a test-by-test basis.

## Simulations of Element Level Tests

The element level tests that are considered in this paper are conventional 1-day oedometer tests on Murro and Haarajoki clays. Modeling 1D compression tests has been performed by simulating conventional oedometer experiments with *Plaxis 2D 2010* and using EVP-SCLAY1 and/or ACM as the user-defined soil model to represent the material behavior.

Fig. 3 shows the model prediction against experimental data, by plotting the predicted and measured vertical strains versus time in semilogarithmic scale for Murro clay at a depth of 3.6 m, which represents the most compressible layer in the field. These are objective predictions, in that no attempt has been made to create a best possible match, and the values of the input parameters (see Table 1) are exactly the same as used further on in the corresponding layer at the boundary value level. Therefore, the soil constants have been determined independently of the test that is being modeled. For all loading stages simulated, EVP-SCLAY1 appears to give better predictions than ACM.

In Fig. 4, similar predictions have been made for Haarajoki clay, considering two different load increments at each depth. Here the results are not as conclusive: for most stress levels and depths, EVP-SCLAY1 appears to give better predictions than ACM, but at depths of 6.35 and 17.22 m, ACM gives better predictions than EVP-SCLAY1 for the vertical effective stress levels of 108 and 330 kPa, respectively. In the following, the models are compared at the boundary value level. Of course, it would be possible to get a better match with the ACM model by, for example, increasing POP, OCR, or reference time or indeed changing the values for the soil constants, to manipulate the value of the creep exponent. That, however, would be curve fitting rather than exploring objectively the model's predictive ability.

<table><caption>**Table 5.** Values of Viscosity Parameters for Murro Test Embankment</caption><thead><tr><th rowspan="2">Layer</th><th colspan="2">EVP-SCLAY1</th><th colspan="2">ACM</th></tr><tr><th>N</th><th>μ (day<sup>-1</sup>)</th><th>μ*</th><th>τ (day)</th></tr></thead><tbody><tr><td>1–1</td><td>25</td><td>8.64 × 10<sup>-8</sup></td><td>8.69 × 10<sup>-4</sup></td><td>1</td></tr><tr><td>1–2</td><td>25</td><td>8.64 × 10<sup>-8</sup></td><td>8.69 × 10<sup>-4</sup></td><td>1</td></tr><tr><td>2</td><td>11</td><td>2.59 × 10<sup>-7</sup></td><td>2.33 × 10<sup>-3</sup></td><td>1</td></tr><tr><td>3</td><td>20</td><td>8.64 × 10<sup>-5</sup></td><td>1.92 × 10<sup>-3</sup></td><td>1</td></tr><tr><td>4</td><td>9</td><td>6.91 × 10<sup>-4</sup></td><td>1.51 × 10<sup>-3</sup></td><td>1</td></tr><tr><td>5</td><td>25</td><td>4.32 × 10<sup>-5</sup></td><td>1.49 × 10<sup>-3</sup></td><td>1</td></tr><tr><td>6</td><td>12</td><td>8.64 × 10<sup>-4</sup></td><td>7.30 × 10<sup>-4</sup></td><td>1</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Time-settlement predictions for Murro test embankment [field data from Karstunen and Yin (2010)]</figcaption></figure>

## Simulations of Murro Test Embankment

To study model performance at the boundary value level, two-dimensional finite-element simulations have been performed with *Plaxis 2D 2010*. The Murro test embankment is 2 m high and 30 m long, and the slopes have a gradient of 1:2 (see Fig. 5 for the geometry and soil layers). The groundwater table is assumed to be at the depth of 0.8 m. The embankment material is crushed rock. The construction was completed within 2 days (Karstunen et al. 2005). The embankment was modeled with a simple linear elastic-perfectly plastic Mohr-Coulomb model using the following typical values for the embankment material: E = 40,000 kN/m², ν' = 0.35, φ' = 40°, ψ' = 0°, c' = 2 kN/m², and γ = 19.6 kN/m³. The values of the soil constant and state variables used are listed in Tables 3–5.

Fig. 6 plots the settlement versus time at the centerline [for the exact location of the instruments, see Karstunen and Yin (2010)].
EVP-SCLAY1 gives a rather good prediction, whereas in contrast ACM grossly overpredicts the vertical deformations. Fig. 7 gives more insight into the reasons for this gross overprediction by ACM. The predicted settlement by the ACM model [Fig. 7(b)] shows non-negligible vertical deformation also outside the embankment area, indicating that significant creep strains are predicted because of the in situ stresses only. This is not realistic, of course. This is also reflected in the predicted horizontal deformations under the crest [Fig. 8(b)] and toe [Fig. 9(b)] of the embankment. EVP-SCLAY marginally underpredicts the horizontal movements under the crest of the embankment [Fig. 8(a)], whereas under the toe of the embankment the predictions are rather good. However, in the field, the horizontal movements develop more slowly than predicted by the rate-dependent models. Under the crest of the embankment, the

location of the maximum horizontal movement is not well predicted, whereas under the toe the depth is predicted with good accuracy.

## Simulations of Haarajoki Embankment

The Haarajoki test embankment is 2.9 m high and 100 m long, with slope gradient of 1:2 (see Fig. 10) and lateral berms for stability. Half of the embankment (50-m-long section) was constructed on an area improved with prefabricated vertical drains, and the other half (Section 35840, considered for this study) was constructed on natural deposits without any ground improvement. The embankment itself was constructed in 0.5-m-thick layers, and each layer was applied and compacted within 2 days. The deposit consists of a 2-m-thick dry crust layer overlying about 22 m of soft clay. The groundwater

<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Surface settlements of Murro test embankment: (a) EVP-SCLAY1; (b) ACM [field data from Karstunen and Yin (2010)]</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Lateral displacement predictions under the crest of Murro test embankment: (a) EVP-SCLAY1; (b) ACM [field data from Karstunen and Yin (2010)]</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Lateral displacement predictions under the toe of Murro test embankment: (a) EVP-SCLAY1; (b) ACM [field data from Karstunen and Yin (2010)]</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Fig. 10. Geometry of Haarajoki embankment</figcaption></figure>

table is assumed to be at the ground surface. Based on the ground investigation data (Finnish National Road Administration 1997), the subsoil is divided into seven sublayers with different compressibility parameters and POP values, as summarized in Tables 6–8. The embankment, which was made of granular fill, was modeled with a simple Mohr–Coulomb model assuming the following material parameters: $E = 40,000 \text{ kN/m}^2$, $v' = 0.35$, $\phi' = 40^\circ$, $\psi' = 0^\circ$, $c' = 2 \text{ kN/m}^2$, and $\gamma = 21 \text{ kN/m}^3$.

The measured and predicted settlements at the centerline of Section 35840, as a function of time, are presented in Fig. 11. EVP-SCLAY1 clearly underpredicts the vertical deformations whereas, yet again, ACM, after a good prediction for a year or so, overpredicts the vertical deformations. Because of the presence of the vertical drains in the other section, it is likely that the field data might not be totally representative of the plane strain geometry assumed in the simulations. Namely, when looking at the vertical deformations along the longitudinal section (see Koskinen et al. 2002), the vertical deformations in Section 35840 appear to be influenced by the section with vertical drains. To confirm that, as discussed by Yildiz (2009), the embankment would need to be modeled in three dimensions, albeit using the symmetry along the centerline.

Unlike in the case of the Murro test embankment, ACM does not predict significant vertical deformations outside the embankment area [Fig. 12(b)]. Obviously, given that the POP is assumed to be >30 kPa throughout the deposit, there is now enough overconsolidation to scale down the creep rates in the nonloaded area. However, as in the case of the Murro test embankment, ACM significantly overpredicts the lateral deformations. This is true for both the crest and the toe [Figs. 13(b) and 14(b)]. In contrast, EVP-SCLAY1 gives good prediction under the toe of the embankment [Fig. 14(a)].

## Conclusions

The paper presents a comparison of two anisotropic rate-dependent models, ACM (Leoni et al. 2008) and EVP-SCLAY1 (Karstunen

<table><caption>Table 6. Initial Values for the State Parameters for Haarajoki Embankment</caption><thead><tr><td>Layer</td><td>Depth (m)</td><td>POP (kPa)</td><td>e<sub>0</sub></td><td>α<sub>K0</sub></td></tr></thead><tbody><tr><td>1</td><td>0–2</td><td>110</td><td>1.25</td><td>0.63</td></tr><tr><td>2</td><td>2–6</td><td>32</td><td>2.9</td><td>0.44</td></tr><tr><td>3</td><td>6–7</td><td>32</td><td>2.6</td><td>0.55</td></tr><tr><td>4</td><td>7–12</td><td>32</td><td>2.35</td><td>0.44</td></tr><tr><td>5</td><td>12–15</td><td>32</td><td>2.2</td><td>0.46</td></tr><tr><td>6</td><td>15–18</td><td>32</td><td>2.0</td><td>0.61</td></tr><tr><td>7</td><td>18–22</td><td>32</td><td>1.25</td><td>0.61</td></tr></tbody></table>

and Yin 2010). Both rate-dependent models are able to account for initial anisotropy and the evolution of anisotropy in a similar manner, analogous to the S-CLAY1 model (Wheeler et al. 2003). However, they have differences in the mathematical formulations for calculating the creep strain rates. Simulations were done at both the element level (1-day oedometer tests) and the boundary value level for Murro and Haarajoki clays. Based on the results, the EVP-SCLAY1 model, which is an overstress model, is able to yield good predictions at both element and boundary value levels, in particular for the case of the Murro test embankment on almost normally consolidated clay. In contrast, ACM severely overpredicted both vertical and horizontal deformations at Murro, as unrealistic vertical deformations were triggered outside the embankment area because of the in situ stresses only. In the case of the Haarajoki embankment, the results were inconclusive: EVP-SCLAY1 underpredicted deformations, whereas ACM overpredicted them. However, the field measurements at Haarajoki are influenced by the section of the embankment that was founded on vertical drains (see Koskinen et al. 2002). Therefore, the embankment would need to be remodeled in three dimensions to take that effect into account.

The analyses presented in this paper ignored the effects of apparent bonding and destructuration for the sake of simplicity. For relatively sensitive clays, such as Murro clay and Haarajoki clay, this may not be appropriate. Indeed, certain phenomena, such as the
<table><caption>**Table 7.** Values for the Conventional Soil Parameters and Additional Soil Constants for Haarajoki Embankment</caption><thead><tr><th rowspan="2">Layer</th><th rowspan="2">γ (kN/m³)</th><th rowspan="2">ν'</th><th rowspan="2">M</th><th rowspan="2">k<sub>x</sub> (m/day)</th><th rowspan="2">k<sub>y</sub> (m/day)</th><th rowspan="2">ω</th><th rowspan="2">ω<sub>d</sub></th><th colspan="2">EVP-SCLAY1</th><th colspan="2">ACM</th></tr><tr><th>λ</th><th>κ</th><th>λ*</th><th>κ*</th></tr></thead><tbody><tr><td>1</td><td>17.5</td><td>0.2</td><td>1.6</td><td>3.46 × 10<sup>-4</sup></td><td>1.73 × 10<sup>-4</sup></td><td>37</td><td>1.02</td><td>0.2</td><td>0.010</td><td>0.089</td><td>8.89 × 10<sup>-3</sup></td></tr><tr><td>2</td><td>14.3</td><td>0.2</td><td>1.15</td><td>1.04 × 10<sup>-4</sup></td><td>5.18 × 10<sup>-5</sup></td><td>33</td><td>0.70</td><td>1.33</td><td>0.010</td><td>0.341</td><td>1.54 × 10<sup>-2</sup></td></tr><tr><td>3</td><td>14.3</td><td>0.2</td><td>1.43</td><td>8.64 × 10<sup>-5</sup></td><td>4.32 × 10<sup>-5</sup></td><td>49</td><td>0.97</td><td>0.96</td><td>0.030</td><td>0.267</td><td>1.17 × 10<sup>-2</sup></td></tr><tr><td>4</td><td>15.1</td><td>0.2</td><td>1.15</td><td>8.64 × 10<sup>-5</sup></td><td>4.32 × 10<sup>-5</sup></td><td>44</td><td>0.70</td><td>0.96</td><td>0.036</td><td>0.287</td><td>1.25 × 10<sup>-2</sup></td></tr><tr><td>5</td><td>15.1</td><td>0.2</td><td>1.20</td><td>8.64 × 10<sup>-5</sup></td><td>4.32 × 10<sup>-5</sup></td><td>35</td><td>0.76</td><td>1.06</td><td>0.030</td><td>0.331</td><td>8.75 × 10<sup>-3</sup></td></tr><tr><td>6</td><td>15.7</td><td>0.2</td><td>1.55</td><td>8.64 × 10<sup>-5</sup></td><td>4.32 × 10<sup>-5</sup></td><td>36</td><td>1.01</td><td>0.45</td><td>0.034</td><td>0.15</td><td>1.23 × 10<sup>-2</sup></td></tr><tr><td>7</td><td>17.5</td><td>0.2</td><td>1.55</td><td>3.46 × 10<sup>-4</sup></td><td>1.73 × 10<sup>-4</sup></td><td>37</td><td>1.01</td><td>0.10</td><td>0.004</td><td>0.044</td><td>8.89 × 10<sup>-3</sup></td></tr></tbody></table>

<table><caption>**Table 8.** Values for Viscosity Parameters for Haarajoki Embankment</caption><thead><tr><th rowspan="2">Layer</th><th colspan="2">EVP-SCLAY1</th><th colspan="2">ACM</th></tr><tr><th>N</th><th>μ (day<sup>-1</sup>)</th><th>μ<sup>*</sup></th><th>τ (day)</th></tr></thead><tbody><tr><td>1</td><td>20</td><td>8.00 × 10<sup>-4</sup></td><td>1.16 × 10<sup>-3</sup></td><td>1</td></tr><tr><td>2</td><td>18</td><td>8.00 × 10<sup>-4</sup></td><td>4.44 × 10<sup>-3</sup></td><td>1</td></tr><tr><td>3</td><td>18</td><td>6.00 × 10<sup>-6</sup></td><td>3.47 × 10<sup>-3</sup></td><td>1</td></tr><tr><td>4</td><td>18</td><td>6.00 × 10<sup>-6</sup></td><td>3.73 × 10<sup>-3</sup></td><td>1</td></tr><tr><td>5</td><td>30</td><td>8.00 × 10<sup>-5</sup></td><td>4.32 × 10<sup>-3</sup></td><td>1</td></tr><tr><td>6</td><td>60</td><td>5.00 × 10<sup>-4</sup></td><td>1.95 × 10<sup>-3</sup></td><td>1</td></tr><tr><td>7</td><td>60</td><td>5.00 × 10<sup>-4</sup></td><td>5.79 × 10<sup>-4</sup></td><td>1</td></tr></tbody></table>

<figure><img src="image_12.png" /><figcaption>**Fig. 11.** Time-settlement predictions of Haarajoki embankment (Section 35840) [field data from Yildiz et al. (2009)]</figcaption></figure>

reduction of undrained shear strength measured under the Murro test embankment after 8 years of consolidation, reported by Koskinen and Karstunen (2006), could be explained only by rate-dependent models that account for destructuration. The same applies to phenomena such as tertiary creep and creep rupture.

Our estimations of compressibility, apparent preconsolidation stress, and the tendency of a natural clay to creep are highly dependent on sample quality. The Swedish and Norwegian piston samplers that were used to extract undisturbed samples of Haarajoki and Murro clays are unlikely to result in high-quality samples. This is likely to be true in particular for Murro clay because of its high silt content. Further investigations on the effect of sample disturbance on rate effects and deformations are needed.

# Appendix. Model Formulations

In the following, compression has been assumed as positive, following the geotechnical tradition.

## Formulation of EVP-SCLAY1 Model

The constitutive model EVP-SCLAY1 (Karstunen and Yin 2010) is based on the overstress theory of Perzyna (Perzyna 1963) and the elastoplastic model SCLAY1 (Wheeler et al. 2003). According to Perzyna's overstress theory, the total strain rate is additively composed of the elastic strain rates and viscoplastic strain rates

$$ \dot{\epsilon}_{ij} = \dot{\epsilon}_{ij}^{e} + \dot{\epsilon}_{ij}^{p} \quad (3) $$

where $\dot{\epsilon} = (i, j)$ component of the total strain rate tensor; and the superscripts $e$ and $p$ = elastic and viscoplastic components. The elastic behavior in the proposed model is assumed to be isotropic in a way similar to the modified Cam clay model (Roscoe and Burland 1968).

An elliptical surface, relating to the current state of preconsolidation, is adopted as the static yield surface

$$ f = \frac{3}{2} [\sigma'_{d} - p'\alpha_{d}] : [\sigma'_{d} - p'\alpha_{d}] - \left(M^2 - \frac{3}{2}[\alpha_{d}]:[\alpha_{d}]\right) (p'_{m} - p')p' = 0 \quad (4) $$

where $\sigma'_{d}$ = deviatoric stress tensor ($\sigma'_{d} = \sigma'_{ij} - p'\delta_{ij}$), with Kronecker's delta ($\delta_{ij} = 1$ for $i = j$ and $\delta_{ij} = 0$ for $i \neq j$); $p' = \text{mean effective stress}$; $\alpha_{d}$ = deviatoric fabric tensor, a dimensionless vector with the same form as the deviatoric stress vector (Wheeler et al. 2003); and $M = \text{slope of the critical state line}$. For the special case of a cross-anisotropic sample in triaxial stress space, the scalar parameter $\alpha = \sqrt{[3/2(\alpha_{d}:\alpha_{d})]}$ defines the inclination of the ellipse of the yield curve in the $p' - q$ plane as illustrated in Fig. 1(a). The dynamic loading surface has an identical shape, but a different size $p'_{m}d$ [see Fig. 1(a)] compared with the static yield surface. To represent the dynamic loading surface, Eq. (4) can be rewritten as

$$ f_d = \frac{(3/2) [\sigma'_{d} - p'\alpha_{d}] : [\sigma'_{d} - p'\alpha_{d}]}{(M^2 - (3/2)[\alpha_{d}]:[\alpha_{d}]) p'} + p' - p'_{m}d = 0 \quad (5) $$

The viscoplastic strain rate in Eq. (3) is assumed to obey an associated flow rule with respect to the dynamic loading surface. The incremental viscoplastic strains in the EVP-SCLAY1 model (Karstunen and Yin 2010) can be calculated as
<figure><img src="image_13.png" /><figcaption>Fig. 12. Surface settlements of Haarajoki embankment: (a) EVP-SCLAY1; (b) ACM [field data from Yildiz et al. (2009)]</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 13. Comparison of lateral displacement predictions at the crest of Haarajoki embankment [field data from Yildiz et al. (2009)]</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 14. Comparison of lateral displacement predictions at the toe of Haarajoki embankment [field data from Yildiz et al. (2009)]</figcaption></figure>

$$ \dot{\varepsilon}_{ij}^{vp} = \mu \left\langle \exp \left[ N \cdot \left( \frac{p_m^{td}}{p_m^{ts}} - 1 \right) \right] - 1 \right\rangle \frac{\partial f_d}{\partial \sigma'_{ij}} \quad (6) $$

where $\mu$ = fluidity parameter; $\langle \rangle$ = McCauley brackets; and $N$ = strain-rate coefficient relating to the strain-rate effect on shear strength and preconsolidation pressure. $N$ and $\mu$ are therefore the key parameters controlling the viscoplastic strain rate. $p_m^{td}$ and $p_m^{ts}$ = size of dynamic loading surface and the static yield surface, respectively [see Fig. 1(a)]; and $f_d$ = viscoplastic potential function, represented by the dynamic loading surface.

The expansion of the static yield surface, which represents the kinematic hardening of the material, is assumed to be due to the viscoplastic volumetric strain rate $\dot{\varepsilon}_v^{vp}$, similarly to the critical state models

$$ \dot{P}_m^s = P_m^s \left( \frac{1+e}{\lambda - \kappa} \right) \dot{\varepsilon}_v^{vp} \quad (7) $$

where $\lambda$ = slope of the normal compression curve in the $e-\ln \sigma'_v$ plane; $\kappa$ = slope of the swelling line; and $e$ = void ratio.

The rotational hardening law that controls the development or erasure of anisotropy caused by viscoplastic strains is defined based on the formulation proposed by Wheeler et al. (2003). Both volumetric and deviatoric viscoplastic strains influence the rotation of the yield surface

$$ \dot{\alpha}_d = \omega \left[ \left( \frac{3\eta_d}{4} - \alpha_d \right) \langle \dot{\varepsilon}_v^{vp} \rangle + \omega_d \left( \frac{\eta_d}{3} - \alpha_d \right) \dot{\varepsilon}_d^{vp} \right] \quad (8) $$

where dots again refer to rates; $\eta_d$ = tensorial equivalent of the stress ratio defined as $\eta_d = \sigma_d/p'$; $\omega$ and $\omega_d$ = additional soil constants ($\omega$ controls the absolute rate of the rotation of the yield and loading surfaces toward their current target values; $\omega_d$ controls the relative effectiveness of viscoplastic volumetric and deviatoric strains in rotating the yield and loading surfaces); and $\dot{\varepsilon}_d^{vp}$ = viscoplastic deviatoric strain rate.
## Formulation of ACM Model

In the ACM model, the total strain rate is additively composed of the elastic strain rates and (irrecoverable) creep strain rates, using superscript c to denote creep

$$ \dot{\epsilon}_{ij} = \dot{\epsilon}_{ij}^{c} + \dot{\epsilon}_{ij}^{c} \qquad (9) $$

In the ACM model (Leoni et al. 2008), a so-called equivalent mean stress $p'_{eq}$ is defined as

$$ p'_{eq} = p' + \frac{(3/2) [\sigma'_{d} - p'\alpha_{d}] : [\sigma'_{d} - p'\alpha_{d}]}{(M^2 - (3/2)[\alpha_{d}] : [\alpha_{d}])p'} \qquad (10) $$

which represents the size of the current stress surface (CSS). This is analogous to the static yield surface in Eq. (4). CSS cannot be called a yield surface in a classic sense, as with the ACM model it is possible to have irrecoverable creep strains everywhere in the stress space. However, within the normal consolidation surface (NCS) of the ACM model, defined with Eq. (10) by substituting $p'_{eq}$ with $p'_{p}$ [see Fig. 1(b)], the creep strains are albeit possible but they are small indeed. In the ACM model, the preconsolidation pressure $p'_{p}$ evolves with the volumetric creep rate according to the hardening law, which is analogous to Eq. (7), but is expressed in terms of modified compression and swelling indexes $\lambda^{*}$ and $\kappa^{*}$, where $\lambda^{*} = \lambda/(1+e_{0})$ and $\kappa^{*} = \kappa/(1+e_{0})$ are the modified compression and swelling indexes, respectively. This can be expressed in integrated format as

$$ p'_{p} = p'_{p0} \exp\left(\frac{\dot{\epsilon}_{v}^{c}}{\lambda^{*} - \kappa^{*}}\right) \qquad (11) $$

When $p'_{eq} = p'_{p}$, the current stress state lies on the normal consolidation surface [Fig. (1b)] and the soil state is normally consolidated, and as a consequence large creep strains occur. The ratio $p'_{p}/p'_{eq}$ (referred to as OCR$^{*}$), gives a measurement of the distance between the current stress surface and the normal consolidation surface, being a generalization of OCR. In the ACM model, the normally consolidated surface is treated as the contour of constant volumetric creep, so the incremental volumetric creep strains are calculated as

$$ \dot{\epsilon}_{v}^{c} = \frac{\mu^{*}}{\tau} \left( \frac{p'_{eq}}{p'_{p}} \right)^{(\lambda^{*} - \kappa^{*})/\mu^{*}} \qquad (12) $$

where is $\mu^{*}$ = modified creep index (slope of $\epsilon_{v}$-ln t plot); and $\tau$ = reference time. As discussed by Leoni et al. (2008), the concept of a reference time makes this type of model very attractive, as it enables scaling of the input for the overconsolidation ratio dependent on the strain rate in the test. The deviatoric component of the creep strain rate vector can be obtained from the flow rule, which is assumed to be associated. The rotational hardening law in the ACM model is similar to that of EVP-SCLAY1 model

$$ \dot{\alpha}_{d} = \omega \left( \left[ \frac{3\eta_{d}}{4} - \alpha_{d} \right] \langle \dot{\epsilon}_{v}^{c} \rangle + \omega_{d} \left[ \frac{\eta_{d}}{3} - \alpha_{d} \right] \dot{\epsilon}_{d}^{c} \right) \qquad (13) $$

## Acknowledgments

The work presented was carried out as part of a Marie Curie Industry-Academia Partnerships and Pathways projects GEO-INSTALL (PIAP-GA-2009-230638) and CREEP (PIAP-GA-2011-286397). The experimental work was done at Aalto University sponsored by the Academy of Finland (Grant 128459).

## References

* Brinkgreve, R. B. J., Swolfs, W. M., and Engin, E. (2010). *Plaxis 2010 reference manual*, Plaxis, Delft, Netherlands.

* Finnish National Road Administration. (1997). "Competition to calculate settlements at the Haarajoki test embankment." Competition programme, competition materials, Finnra, Helsinki, Finland.

* Hinchberger, S. D., and Qu, G. (2009). "Viscoplastic constitutive approach for rate-sensitive structured clays." *Can. Geotech. J.*, 46(6), 609-626.

* Karstunen, M., and Koskinen, M. (2004). "Anisotropy and destructuration of Murro clay." *Proc., Advances in Geotechn. Eng. Skempton Conf.*, Vol. 1, Thomas Telford, London, 476-487.

* Karstunen, M., and Koskinen, M. (2008). "Plastic anisotropy of soft reconstituted clays." *Can. Geotech. J.*, 45(3), 314-328.

* Karstunen, M., Krenn, H., Wheeler, S. J., Koskinen, M., and Zentar, R. (2005). "Effect of anisotropy and destructuration on the behavior of Murro test embankment." *Int. J. Geomech.*, 10.1061/(ASCE)1532-3641 (2005)5:2(87), 87-97.

* Karstunen, M., and Yin, Z.-Y. (2010). "Modelling time-dependent behaviour of Murro test embankment." *Geotechnique*, 60(10), 735-749.

* Koskinen, M., and Karstunen, M. (2006). "Numerical modelling of Murro test embankment with S-CLAY1S." *Proc., 6th European Conf. Numerical Methods Geotechnical Engineering*, Taylor & Francis, London, 11-19.

* Koskinen, M., Karstunen, M., and Wheeler, S. J. (2002). "Modelling destructuration and anisotropy of a natural soft clay." *Proc., 5th European Conf. on Numerical Methods in Geotechnical Engineering (NUMGE)*, Presses de l'ENPC, Paris, 11-20.

* Leoni, M., Karstunen, M., and Vermeer, P. A. (2008). "Anisotropic creep model for soft soils." *Geotechnique*, 58(3), 215-226.

* Perzyna, P. (1963). "The constitutive equations for work-hardening and rate sensitive plastic materials." *Proc. Vibr. Prob.*, 4(3), 281-290.

* Qu, G., Hinchberger, S. D., and Lo, K. Y. (2010). "Evaluation of the viscous behaviour of clay using generalised overstress viscoplastic theory." *Geotechnique*, 60(10), 777-789.

* Roscoe, K. H., and Burland, J. B. (1968). *On the generalized stress-strain behaviour of 'wet' clay*, Vol. I, Cambridge University Press, Cambridge, U.K., 553-609.

* Wheeler, S. J., Näätänen, A., Karstunen, M., and Lojander, M. (2003). "An anisotropic elasto-plastic model for soft clays." *Can. Geotech. J.*, 40(2), 403-418.

* Yildiz, A. (2009). "Numerical modeling of vertical drains with advanced constitutive models." *Comput. Geotech.*, 36(6), 1072-1083.

* Yildiz, A., Karstunen, M., and Krenn, H. (2009). "Effect of anisotropy and destructuration on behavior of Haarajoki test embankment." *Int. J. Geomech.*, 10.1061/(ASCE)1532-3641(2009)9:4(153), 153-168.

* Yin, Z., Karstunen, M., Chang, C., Koskinen, M., and Lojander, M. (2011). "Modeling time-dependent behavior of soft sensitive clay." *J. Geotech. Geoenviron. Eng.*, 10.1061/(ASCE)GT.1943-5606.0000527, 1103-1113.

* Yin, Z.-Y., Chang, C. S., Karstunen, M., and Hicher, P. Y. (2010). "An anisotropic elastic-viscoplastic model for soft clays." *Int. J. Solids Struct.*, 47(5), 665-677.

* Zentar, R., Karstunen, M., and Wheeler, S. J. (2002). "Influence of anisotropy and destructuration on undrained shearing of natural clays." *Proc., 5th European Conf. on Numerical Methods in Geotechnical Engineering (NUMGE)*, Presses de l'ENPC, Paris, 21-26.

* Zhou, C., Yin, J., Zhu, J., and Cheng, C. (2005). "Elastic anisotropic viscoplastic modeling of the strain-rate-dependent stress-strain behavior of K0-consolidated natural marine clays in triaxial shear tests." *Int. J. Geomech.*, 10.1061/(ASCE)1532-3641(2005)5:3(218), 218-232.