An anisotropic elastoplastic model for soft clays based on logarithmic contractancy

A new constitutive model for soft structured clays is developed based on an existing model called S‐CLAY1S, which is a Cam clay type model that accounts for anisotropy and destructuration. The new model (E‐SCLAY1S) uses the framework of logarithmic contractancy to introduce a new parameter that controls the shape of the yield surface as well as the plastic potential (as an assumed associated flow rule is applied). This new parameter can be used to fit the coefficient of earth pressure at rest, the undrained shear strength or the stiffness under shearing stress paths predicted by the model. The improvement to previous constitutive models that account for soil fabric and bonding is formulated within the contractancy framework such that the model predicts the uniqueness of the critical state line and its slope is independent of the contractancy parameter. Good agreement has been found between the model predictions and published laboratory results for triaxial compression tests. An important finding is that the contractancy parameter, and consequently the shape of the yield surface, seems to change with the degree of anisotropy; however, further study is required to investigate this response. From published data, the yield surface for isotropically consolidated clays seems ‘bullet’ or ‘almond’ shaped, similar to that of the Cam clay model; while for anisotropically consolidated clays, the yield surface is more elliptical, like a rotated and distorted modified Cam clay yield surface. © 2015 The Authors. International Journal for Numerical and Analytical Methods in Geomechanics published by John Wiley & Sons Ltd.


INTRODUCTION
Extensive experimental testing of soils under different stress paths and conditions as well as the increase in computing power has led to the development of advanced constitutive models that reproduce more accurately the mechanical behaviour of soils. Since the pioneering work of Roscoe and co-workers [1][2][3], many constitutive models have been proposed in the framework of critical state soil mechanics. The first critical state model was the (original) Cam clay (CC) model, whose plastic potential surface was obtained on the basis of assuming a simple frictional form for the plastic work. Associated flow conditions were assumed, and therefore, the yield and plastic potential surfaces coincide.
The original model was modified (modified Cam clay [MCC]) [4] using a different formulation of the dissipated energy during plastic straining to get an elliptical yield surface that overcomes some of the limitations of the original surface, for example, the singularity on the mean effective stress axis (q = 0). Further yield and plastic potential surfaces have been proposed in the literature since then, for example, [5,6]. Hence, the shapes of the yield and plastic potential surfaces vary from model to model. A constitutive model that is able to reproduce a variety of shapes (yield surfaces) could provide predictions that are more accurate. Following that idea, Lagioia et al. [7] developed a versatile expression for the yield and plastic potential surfaces based on a mathematical relation between the dilatancy and the stress ratio. However, this model is limited to isotropic soils and does not account for natural soil features such as fabric and some apparent bonding that will be progressively lost during loading.
Natural soft soil deposits exhibit inherent anisotropic behaviour because of its deposition history. Therefore, the extent of anisotropy can be modelled by a rotated and distorted elliptical yield surface (e.g. Dafalias [8], S-CLAY1 model [9], MIT-E3 model [10] or Sekiguchi-Ohta model [11]) by considering inter-dependence (coupling) of volumetric and deviatoric plastic strains in the plastic work equation. The major differences between these models are shape of the yield surface and rotational hardening rule. The S-CLAY1 model [9] was further extended to include soil structure (S-CLAY1S [12]) through an intrinsic yield surface. The S-CLAY1S model has proven its ability to reproduce the behaviour of normally or slightly overconsolidated structured soft clays [12][13][14][15]. Despite its good performance, especially for settlement prediction, horizontal displacements are generally not well matched, for example, [15]. Those differences may be attributed to the shape of the yield surface (i.e. associates with flow rule), or similarly, to the horizontal/vertical stress ratio predicted by the model for compression loading.
Ohta et al. [16] presented a unified framework for different shapes of the yield surface, assuming associated flow conditions. The framework is based on curve fitting of experimental results of the contractancy (compressive volumetric strain, ε v ) during drained shear at constant mean effective stress (p′) of normally consolidated clays. Those experimental results were first presented by Shibata [17]. Ohno et al. [18] proposed two categories of contractancy models, namely, exponential and logarithmic contractancy models, depending on the type of function used to fit the experimental results. (Original) CC and MCC models are particular cases of the general contractancy models. Ohta et al. [16] extended the contractancy models to anisotropic conditions using the stress parameter introduced by Sekiguchi and Ohta [11]. Therefore, they are called extended Sekiguchi-Ohta models [16]. These models do not account for evolution of anisotropy with loading and apparent bonding of natural soils. Further, the extended Sekiguchi-Ohta models fail to predict a unique critical state line (CSL) as for example in triaxial compression and extension in p'-q space.
According to Dafalias and Taiebat [19], introducing a non-associated flow rule can provide improved predictions regardless of the rotational hardening rule employed, while it is able to obtain a unique CSL. However, the introduction of a non-associated flow rule may be either not necessary or convenient. Therefore, development of an adequate but simple constitutive model for anisotropic and structured clays is still relevant.
This paper extends the S-CLAY1S model [12] using the framework of logarithmic contractancy [16] to include some flexibility in the shape of the yield surface. The new model introduces an additional parameter, called the contractancy parameter, which controls the shape of the yield surface. The contractancy parameter can be related to the coefficient of earth pressure at rest for normally consolidated conditions, K 0NC , the undrained shear strength, c u , or the stiffness under shearing stress paths. In this way, the proposed model, called E-SCLAY1S, extends the predictive capabilities of the S-CLAY1S model, while including just an additional parameter with clear physical meaning.
The paper presents the formulation of the new model (Section 2) and its numerical implementation (Section 3). Section 4 highlights the main features of the model, such as the slope and uniqueness of the CSL and the influence of the contractancy parameter on the coefficient of lateral earth pressure at rest, the undrained shear strength, the yield surface and the soil stiffness. The model is validated against laboratory tests on two clays, namely, Kaolin clay and Santa Clara clay (Section 5) and finally some discussion and conclusions are provided.

PROPOSED MODEL: E-SCLAY1S
The proposed model extends S-CLAY1S [12], which is a MCC-type model [4] that accounts for anisotropy and destructuration. Anisotropy of plastic behaviour is represented through an inclined and distorted yield surface and a rotational hardening law to model the development or erasure of fabric anisotropy during plastic straining; while interparticle bonding and degradation of bonds (structure) is reproduced using intrinsic and natural yield surfaces [20] and a hardening law describing destructuration as a function of plastic straining. For the sake of simplicity, the mathematical formulation is presented in the following in triaxial stress space, which can be used only to model the response of cross-anisotropic samples (cut vertically from the soil deposit) subject to oedometer or triaxial loading. The original inclined yield surface of the S-CLAY1S model is elliptical [8] (Figure 1): The preceding yield function often cannot describe experimental data of yield points with enough accuracy as well as undrained stress paths [19]. An improvement can be achieved by modifying the yield function. The proposed model (E-SCLAY1S) introduces a degree of freedom in the shape of the yield surface (Eq. (1)) using the framework of logarithmic contractancy [16] (Appendix A): where Ψ is an intermediate parameter to simplify Eq. (2) and n L is a new parameter (contractancy parameter) that controls the shape of the yield surface. The subscript L refers to logarithmic contractancy, following the notation by Ohta et al. [16]. So, the shape of the yield surface in the E-SCLAY1S model depends on the contractancy parameter, n L (Figures 2 and 3). It is worth mentioning that it corresponds to the shape of both the natural and intrinsic yield surfaces and also the plastic potential surface, because an associated flow rule is assumed. For the sake of brevity, hereafter, it is referred to as the yield surface. The E-SCLAY1S preserves the hierarchical development of S-CLAY1S, as the former reduces to the later for n L = 2, that is, Eq. (2) reduces to Eq. (1) as Ψ = 1.

598
For the E-SCLAY1S model, and in general for logarithmic contractancy models, the yield surface is elliptical for n L = 2, 'bullet' shaped for n L < 2 and 'tear' shaped for n L > 2. Most of the anisotropic models (e.g. S-CLAY1S and MIT-E3) use elliptical surfaces, but empirical evidences show the  limitations associated with elliptical yield surfaces (e.g. [21]). Similarly to (original) CC model, the yield surfaces for n L < 2 have a singularity at η = α ( Figure 2). As in the S-CLAY1S model, the effect of bonding in the E-SCLAY1S model is described by an intrinsic yield surface [20], which has the same shape and inclination of the natural yield surface but with a smaller size ( Figure 1). The size of the intrinsic yield surface is specified by the state parameter p ′ mi , which is related to the size p ′ m of the natural yield surface by the state parameter χ as the current amount of bonding The last letter of the model ('S') refers to the soil structure. So, when the hierarchical version of the model without destructuration is used, the model is simply called E-SCLAY1.
The three hardening rules of the original S-CLAY1S model, namely, isotropic hardening, rotational hardening and degradation of bonds rule, are kept as the original [9,12]. The first rule relates the increase in the size of the intrinsic yield surface to the increments of plastic volumetric strain (dε p v ) where ν is the specific volume, λ i is the gradient of the intrinsic normal compression line (NCL) in the compression plane (ln p′ À ν space) and κ is the slope of the swelling line in the compression plane. The second hardening law is the rotational hardening law, which describes the rotation of the yield surface with plastic straining where η is the tensorial equivalent of the stress ratio defined as η = σ d /p′, dε p d is the increment of plastic deviatoric strain, α d is the deviatoric fabric tensor, which has the same form as the deviatoric stress vector [9], and ω and ω d are additional soil constants that control rotational hardening.
The third law for destructuration is formulated in such a way that both plastic volumetric strains and plastic shear strains tend to decrease the value of the bonding parameter χ towards a target value of zero; it is defined as where ξ and ξ d are additional soil constants. As full details of the hardening laws and determination of the soil constants is beyond the scope of this paper, they can be found in [9,12].
The extension of the model from triaxial stress space to general (multiaxial) stress space is also equivalent to that of S-CLAY1S [22]. The model has been implemented using the Euler backward implicit integration scheme [22], in such a way that it can be incorporated directly into finite element codes (e.g. PLAXIS [23]) for engineering applications. The implementation is presented in the following section.

DISCRETIZATION AND NUMERICAL IMPLEMENTATION
The decomposition of total strains in classical elasto-plastic theory using an additive rule can be expressed in terms of elastic and plastic components of strains as where d remarks an incremental operator, the boldface characters are used to denote tensor quantities and superscript e denotes the elastic component and p denotes the plastic component. The plastic strain increment can be calculated using the plastic multiplier (dΛ) The plastic potential is the yield function because an associated flow rule is assumed.
To derive the plastic multiplier of the E-SCLAY1S model, the consistency condition of the yield function (df y = 0) is used By substituting stress increment and isotropic and rotational hardening rules as well as destructuration law, the plastic multiplier is derived as For a small increment in implicit integration scheme, it can be further simplified in terms of the value of the yield function ( f 0 y ) as where superscript T corresponds to a matrix transpose and hardening moduli Η 0 , Η α and Η χ are derived as Using the Euler backward implicit integration scheme, the trial stress is modified under consideration of the occurring plastic strains as long as convergence is reached. The convergence criterion is fulfilled when the iterative stress state returns to the yield surface ( f y < 10 À7 ). If plasticity is associated with a given strain increment, it is essential to calculate the following system of equations By using the plastic multiplier in Eq. (13), given the strain increment is applied to arrive at the elastic predictor, the stress increment dσ′ can be calculated as In this implementation, size of load (strain) increment is controlled by sub-stepping within the subroutine in order to obtain solutions. The maximum strain increment used to simulate the results presented in this paper is |dε| < 0.1 %. Figure 4 presents a summary of the Euler backward implicit algorithm to implement the proposed model.

Slope and uniqueness of the CSL
The strength at the ultimate state after large strains have developed, that is, at critical state, is controlled by the plastic potential surface. For isotropic contractancy models, the slope of the CSL is equal to M (e.g. [18]). However, for anisotropic contractancy models, the slope of the CSL in stress space ( p′q) is not usually M (e.g. extended Sekiguchi-Ohta models [16]), and it depends on both M and n L . The E-SCLAY1S model has been developed to preserve M as the slope of the CSL in stress space, which is a main advantage with respect to the extended Sekigucha-Ohta models. The slope of the CSL is demonstrated in the following.
At critical state condition, there are no plastic volumetric strains, so and rearranging terms in Eq. (20) the slope of the CSL is given by η Eq. (22) may be developed substituting the value of Ψ (Eq. (3)).
From Eq. (23), it is clear that the slope of the CSL for triaxial compression is η = M. For triaxial extension, the slope of the CSL seems to depend not only on M but also on n L , see for example, the shape of the yield surface in p′-q plots ( Figure 2). However, for triaxial extension, the yield surface rotates towards the extension side, and at critical state, the inclination of the yield surface (α) is on the extension side, that is, α is negative, and Eq. (23) gives the same result for extension and for compression (η = M). Only for the unrealistic case of initial fabric anisotropy but without evolution of that fabric during plastic straining, that is, deactivating the rotational hardening rule, the slope of the CSL would be different for triaxial extension depending on n L in stress space: The constant slope of the CSL in stress space for any value of n L , if rotational hardening is allowed, is also valid for any other direction in the π-plane, that is, for any Lode's angle. In the current implementation of the model, no attempt has been made to distinguish between compression and extension, that is M c =M e , or include any dependency on Lode's angle, so it corresponds to the Drucker-Prager criterion ( Figure 3).
The E-SCLAY1S model also preserves the uniqueness of the CSL as illustrated in Figure 5, where the solid lines with arrows on p′-axis show the uniqueness of the CSL in stress space. To highlight and confirm the uniqueness of the CSL in stress space, undrained triaxial compression and extension tests were simulated from a K 0 consolidated state ( Figure 6). The rotation of the initial yield surface to the extension side during triaxial extension tests leads to the same slope both in compression and in extension, which is equal to M (Figure 6(a)). If the yield surface is not allowed to rotate, that is, no rotational hardening (ω = 0), the slope of the CSL at extension depends on the n L value as previously mentioned ( Figure 6(b)). However, this is a very unrealistic case, which should not be modelled and has been presented only for illustrative purposes.
The rotational hardening rule of the S-CLAY1 model was developed to predict a unique CSL in the e-ln p′ space. In the case of the proposed E-SCLAY1S model, the uniqueness of the CSL in the e-ln p′ space is preserved for a given n L value too. The vertical separation from the isotropic NCL to the CSL in the e-ln p′ space predicted by the E-SCLAY1S model is given by where e N is the void ratio on the NCL that corresponds to a unit stress and e Г is the void ratio of the CSL at a unit stress. χ d (M) is the predicted unique inclination of the yield curve at critical states, which is equal to M/3 [9]. For n L = 2 (S-CLAY1), the previous equation reduces to equation (6) of Wheeler et al. [9], and for n L = 2 and α = 0 (MCC), the previous equation reduces to e N À e Γ = (λ À κ)ln 2. The normalized vertical separation e N Àe Γ λ i Àκ to the NCL in the e-ln p′ space at critical state with the n L value is presented in Figure 7. For comparison, (original) CC, MCC and S-CLAY1 are also presented in the figure. It can be seen that the proposed model has flexibility in predicting the CSL in the e-ln p′ space compared with the previously developed models. According to Wheeler et al. [9], the experimental data from tests on Otaniemi clay do not provide evidence for a unique CSL in the e-ln p′ space.  The (original) CC model clearly overpredicts the coefficient of lateral earth pressure at rest for normally consolidation conditions, K 0NC [24]. Although the MCC model predicts more realistic values than the (original) CC model, it is well known to still overpredict K 0NC (e.g. [24]). Therefore, some authors (e.g. Federico et al. [25]) have used plastic potential surfaces with higher degrees of freedom to fit the desired K 0NC value. Federico et al. [25] also present a wide analysis of the analytical expressions that give K 0NC for isotropic critical state models. The analytical expression that gives the value of K 0NC for E-SCLAY1S is derived in Appendix B. The variation of the K 0NC stress path with n L is illustrated in Figure 8 for an initially isotropically consolidated soil sample. The isotropic version of E-SCLAY1 (α = ω = 0) is a hierarchical extension of MCC that introduces the additional parameter n L , which may be correlated with K 0NC using Eq. (B.7). The variation of K 0NC with n L is shown in Figure 9(a). Although a perfect fit of n L could be applied for each case, a value of n L around 3.5 gives similar values to those of Jaky's empirical formula (K 0NC = 1-sin ϕ). However, the yield function for that value (n L = 3.5) could lead to unrealistic high undrained shear strengths (Figure 9(b)). To avoid that, a non-associated flow rule could be proposed, using a n L value lower than 2 for the yield surface (e.g. 1.3) and a high n L value (e.g. 3.5) for the plastic potential surface. Nevertheless, the isotropic version of E-SCLAY1 has been presented to show its similar capabilities to other previous studies (e.g. [25]), but even initially remoulded soils show some fabric under one dimensional compression (e.g. [26,27]), the authors believe that trying to fit K 0NC values with isotropic plastic potential surfaces is not realistic. Once soil anisotropy is introduced, the proper K 0NC can be fitted adjusting the inclination of the yield surface, α 0 (e.g. [9]), because it is difficult to have enough accurate data to determine the initial inclination of the yield surface (α 0 ), and this inclination seems to be related to one dimensional compression of the soil during its deposition through K 0NC . The E-SCLAY1 model, through the n L parameter, introduces more flexibility in the possibilities to fit K 0NC and the initial inclination of the yield surface (α 0 ). However, in practical situations, there are not enough data about the initial inclination of the yield surface, and the practical approach here proposed is to fit all the parameters similarly to S-CLAY1 [9] and the additional parameter of the E-SCLAY1 model (n L ), using the undrained shear strength.

Undrained shear strength
One of the most important risks associated with the numerical simulation of geotechnical engineering problems under undrained conditions using soil constitutive models formulated in effective stresses is the possible unrealistic prediction of the undrained shear strength. That occurred, for example, in the numerical simulation of the deep excavation near the Nicoll Highway, Singapore, which collapsed in 2004 (Whittle and Davies [28]). The additional parameter (n L ) of the E-SCLAY1 model allows for a perfect matching of the undrained shear strength (c u ) as it will be shown in the comparison with laboratory measurements in Section 5. The variation of c u with n L may be seen in Figure 6(a).
For an initially isotropically and normally consolidated soil sample, c u may be normalized by the initial mean effective pressure, p′ 0 . Figure 10

Yield loci
The E-SCLAY1 model, through the n L parameter, provides an additional flexibility in comparison with the conventional S-CLAY1 model to fit the initial yield surface. However, in practical situations, there is little information about the initial yield surface, and even for the well-documented cases, there are some problems associated with its determination, for example, regarding the homogeneity of the natural soil samples and the identification of soil yielding (e.g. [27]). Most published yield loci (e.g. McGinty [29]) have been determined using bilinear interpretation of e-ln p′. This methodology is similar to the Casagrande method used to calculate the preconsolidation pressure and is summarized, for example, in Graham et al. [30]. The methodology involves some ambiguity because the behaviour of clay is non-linear except at very small strains and some engineering judgement is necessary. Besides, for soils with evolving anisotropy and different loading stress ratios, the methodology is not appropriate (e.g. McGinty [29]). For loading stress ratios that deviate from the initial loading stress ratio, yield curve rotation starts to develop under small volumetric strains before the strains get larger because of isotropic hardening. Consequently, the bilinear interpretation of e-ln p′ curves tend to overestimate yield stresses for stress paths that notably deviate from the initial one (e.g. Figure 11). To overcome these limitations, arithmetic stress scale (e.g. e-p′) is generally used (e.g. [27,31]).
Being aware of the limitations of most published yield loci and as an example of the improved capabilities of the E-SCLAY1 model, the yield surfaces of some well-documented soils are fitted in   Figure 12. The initial inclination of the yield surface was determined as proposed by Wheeler et al. [9] based on M, that is, K 0NC determined using Jaky's expression. The experimental data are taken from Graham et al. [32], Wheeler et al. [9] and Díaz-Rodríguez et al. [33]. No attempt to get the best fit was made, and only for illustration, the yield surfaces using n L = 2.5, which match better the experimental data than those using n L = 2.0, are presented in Figure 12. The size of the yield surface, p′ m , was kept constant.

Soil deformation
For normally consolidated conditions under constant stress ratios (η), the parameter λ controls compressive volumetric strains because E-SCLAY1 is as a CC model. For shearing stress paths, compressive volumetric strains depend on the additional contractancy parameter (n L ). As an illustrative example, Figure 13(a) shows the results of simulated drained triaxial compression tests for different n L values. The soil response is stiffer for higher n L values. The influence of n L on the soil stiffness during shearing is summarized in Figure 13( A higher n L value reduces the space between the NCL and the CSL (Eq. (24) and Figure 7), and, therefore, the soil response is stiffer for shearing stress paths. The proposed model (E-SCLAY1S) is a logarithmic contractancy model [18], which means that it uses a logarithmic description of the compressive volumetric strains, ε v , during drained shear at constant mean effective stress and Figure 12. Yield surfaces for several clays: (a) Winnipeg clay (data after [32]), (b) Otaniemi clay (data after [9]), (c) Drammen clay (data after [33]) and (d) Pornic clay (data afer [33]). CSL, critical state line.
normally consolidated conditions (Eqs. (A.1) and (A.4)). So, n L controls the volumetric strains during shearing and its variation with the stress ratio, η. Ideally, n L could be calibrated by fitting experimental laboratory results of drained triaxial shear tests at constant mean effective stress. The influence of n L on the contractancy results is shown in Figure 14. Volumetric strains (ε v ) are normalized by the volumetric strain at critical state (ε vM ), and the stress ratio (η) is normalized by the stress ratio at critical state (M) to isolate the influence of n L from other model parameters. As an example, laboratory data on isotropically consolidated remoulded kaolin clay by Hattab and Hicher [34] are also presented in Figure 14.

Destructuration
For the sake of completeness, the model has been formulated including soil structure and loss of bonding. This model feature performs similarly to that of the S-CLAY1S model [12], so it has not been considered necessary to include here any specific simulation or comment about it. the well-studied S-CLAY1S model, the focus here is on the improvement provided by the additional contractancy parameter (n L ). All the parameters of the model but n L coincide with the S-CLAY1S model, so the approach proposed by Wheeler et al. [9], which gives satisfactory results for most cases (e.g. [12][13][14][15]), is here used to determine those parameters. The additional parameter n L has been vary to get a better fit of the experimental results, usually, of the undrained shear strength.

Kaolin clay
Stipho [26] conducted a series of undrained triaxial tests on isotropically and anisotropically consolidated specimens of Kaolin clay. Several degrees of overconsolidation (from 1 to 4 or 12) and initial anisotropy (K 0 ranging from 1 to K 0NC = 0.57) were used. The tests were stress controlled, and consequently, failure may not have been properly captured. Several researchers have used the results for verifying their constitutive models (e.g. [35,36]). The parameters used in the numerical simulations are summarized in Table I. They have been directly taken either from previous studies [35,36] for standard critical state parameters or following Wheeler et al. [9] for anisotropy. The parameter that controls the rotation of the yield surface with the plastic strains, ω, was set equal to a very low value (0.5), typical for remoulded Kaolin clay. Only the contractancy parameter was fitted to get a better agreement with the experimental results, particularly with c u . The best fit of n L is compared with S-CLAY1 (n L = 2) in Figure 15. For isotropically consolidated samples, the agreement is very good for n L = 1.3. In addition to stress paths, the stress-strain curves and the generated pore pressures are also well predicted by the model (Figure 16).  For anisotropically consolidated samples (K 0 < 1), it was necessary to gradually increase the value of n L to get a good matching of the experimental results ( Figure 15). For normally consolidated samples at K 0 = 0.67 and 0.57, the matching was not possible, and by observing the stress paths, it can be deduced that the waiting times after consolidation could have slightly overconsolidated the soil samples because of ageing or creep effects. A good agreement was found for S-CLAY1 (n L = 2), using the best fit value of overconsolidation ratio (OCR) because of ageing (1.1 for K 0 = 0.67 and 1.2 for K 0 = 0.57). For the sake of comparison, the stress paths for OCR = 1 are also included in Figure 15, and the starting points are the same because the applied initial stresses were p′/p′ 0 = 1.
The initial rotation of the yield surface (α 0 ) for anisotropically consolidated samples was obtained by simulation of the K 0 stress paths with n L = 2. The proposed model (E-SCLAY1) is an extension of S-CLAY1, and therefore, it shares some of the limitations, such as a good behaviour only for normally or slightly overconsolidated soft soils. For high degrees of overconsolidation, the flexibility introduced by n L improves the numerical predictions only for isotropically consolidated samples.

Santa Clara clay
Venda Oliveira and Lemos [27] present laboratory results of a sandy lean clay from Santa Clara dam area, Portugal. The clay was reconstituted prior to testing. They performed triaxial tests on isotropically and K 0 consolidated clay samples to evaluate several elastoplastic models. Stress path controlled drained triaxial tests were performed to determine the position of the initial yield surface and the direction of plastic strain increments (dε p ). The specimens used to study the isotropic behaviour were initially normally consolidated to an isotropic effective stress (p′ m ) of 200 kPa and subsequently unloaded and consolidated to an effective isotropic pressure of 100 kPa, which resulted in an OCR (p′ m /p′) of 2.0. Four drained triaxial tests with stress paths, dq/dp′, equal to 1.0, 2.5, 3.0 and 5.0 were then performed ( Figure 17).   The specimens used to analyse the anisotropic behaviour were initially subjected to K 0 consolidation (K 0 = 0.47) to reach a normally consolidated state of a vertical effective stress of σ′ vc = 200 kPa (σ′ hc = 94 kPa). Then, the specimens were unloaded and consolidated to a vertical effective stress of σ′ v0 = 160 kPa (σ′ h0 = 80 kPa), corresponding to an OCR (σ′ vc /σ′ v0 ) of 1.25. Four drained triaxial tests with stress paths, dq/dp′, equal to À1.5, À0.5, 1.0 and 3.0 were then performed ( Figure 18).
Careful evaluation of the yield loci based on both ε 1 -p′ and ε v -p′ plots was performed, and the probable limits of the yield zone are provided (inverted triangles). They also performed undrained compression triaxial tests to evaluate undrained stress paths (blue crosses) under normally consolidated conditions. Using the proposed model (E-SCLAY1), the results of those undrained compression triaxial tests were numerically simulated. The soil parameters for the numerical model are shown in Table I and were directly taken from [27]. Figure 17 compares the laboratory results with the numerical predictions for the isotropically consolidated Santa Clara clay. For the numerical simulations, two n L values, namely, n L = 2 (S-CLAY1) and n L = 1.3 (best-fit value), were used. Although the undrained stress path measured in the laboratory (blue crosses) is slightly irregular at the beginning, the best fit value (n L = 1.3) provides a good match of the laboratory results and notably improves the results for n L = 2. The initial yield surfaces of the model for n L = 1.3 ('bullet' shape) and 2 (elliptical) are also shown for evaluation against the yield zone (inverted triangles) and plastic strain increments (arrows). The soil exhibits some rotational hardening (ω = 200 and ω d = 0.91), which causes some deviation of the numerically simulated undrained stress paths from the initial yield surfaces. The agreement between the limits of the yield zone and the initial yield surface for n L = 1.3 is not as good as for the undrained stress path, which may be explained by the difficulties associated with the determination of yield loci. Figure 18 compares the results for the K 0 consolidated Santa Clara clay. In this case, it is difficult to get a good fit of the experimental undrained stress path (blue crosses), and the best fit value (n L = 1.8) was determined to match the undrained shear strength. As for the Kaolin clay, it was necessary to use a different n L value (1.8) than that of the isotropic case (1.3). For anisotropic conditions, the differences with S-CLAY1 (n L = 2) are not very important. As for the initially isotropic case, the initial yield surfaces of the model are also shown for evaluation against the yield zone and plastic strain increments. Plastic strain increments are plotted in (Figures 17, 18) for completion, but it is worth noting that plastic strain vectors are difficult to determine in practice because it is necessary to assume an elastic law and the strain increments need to be small but at the same time large enough to eliminate noises in the measurements. Figure 18. K 0 -consolidated Santa Clara clay (data after [27]). CSL, critical state line.

DISCUSSION ON THE SHAPE OF THE YIELD SURFACE
The contractancy parameter (n L ) controls the shape of the yield surface, yet, in practical situations, n L may be conveniently calibrated to fit c u or the stiffness along shearing paths. From the comparison with experimental results (Figures 15-18), it seems that the n L value may not be constant and may depend on the degree of anisotropy because for isotropically consolidated soils, its value is around 1.3, while for anisotropic conditions it is close to n L = 2, which corresponds to S-CLAY1. Although there may be some uncertainties related to the quality of the experimental data, the shape of the yield function for isotropic conditions seems to be closer to the original CC model than to the MCC model, while for anisotropic conditions a rotated and distorted yield surface seems appropriate (e.g. S-CLAY1). As an additional example, Figure 19 shows experimental values of yield stresses for Bothknennar clay (data after McGinty [29]). For intact soil samples (Figure 19 (a)), the yield surface of the S-CLAY1 model fitted through the experimental data points using M = 1.4 and α K0 = 0.31 (after McGinty [29]). The size of the yield surface ( p ′ m = 85 kPa) was obtained by McGinty [29] optimizing the best fit to the experimental data. As explained in Section 4.4, the yield points (black dots) have been determined using bilinear interpretation of e-ln p′ curves, and this may lead to an overestimation of the yield stress for those stress paths that cause significant rotation of the yield surface, for example, triaxial extension for this case.  Consequently, a good fitting of the yield points for triaxial extension could only be possible, introducing a different slope of the CSL for extension (M e ), that is, introducing a Lode's angle dependent failure criterion.
McGinty [29] also checked the yielding points of Bothkennar clay after isotropic consolidation (Figure 19(b)). The single square point indicates the maximum stress in the common first loading stage (210 kPa), while circular points represent yield points identified from individual second loading stages. The yield points in Figure 19(b) are reasonably symmetric about the p′-axis, suggesting that, as expected, the isotropic loading in the first stage had rotated the yield curve clockwise to an isotropic orientation, that is, symmetrical about the p′-axis. To improve the accuracy of the yield points, the authors have reinterpreted McGinty [29] data using arithmetic stress scale and volumetric and axial engineering strains ( Figure 20). Only compression tests and isotropically consolidated samples have been reinterpreted because those are particularly relevant for the comparison here presented about the shape of the yield surface. In the reinterpretation, instead of a yield point, a yield zone has been identified. This yield zone has been included in Figure 19(b) as a line between crosses. In most tests, an initial non-linear stress-strain behaviour is observed for low stresses. After that, a linear part may be identified, Figure 20. Reinterpretation of yield zone using arithmetic stress scale. Bothkennar clay, isotropically consolidated samples, compression tests (data after [29]). and later, the data show the initiation of an exponential curve, which marks the yield stress ( Figure 20). The linear part has been fitted by a straight line (dashed) to precisely identify the initiation of the exponential behaviour. A yield zone was chosen, as the results for volumetric and axial strains seem to be slightly different. Although McGinty [29] yield points tend to generally underestimate the yield stress, for example, for isotropic consolidation (Figure 19(b)), the reinterpretation confirms that the yield stress is overestimated for stress paths that notably differ from the initial one, if a logarithmic stress scale is used, for example, for η 2 = 1.3. The yield curve expression for S-CLAY1 (n L = 2) with α K0 = 0.0 and p ′ m = 210 kPa (corresponding to the MCC yield curve expression) is a poor match to the experimental data ( Figure 19(b)). The E-SCLAY1S model using n L = 1.4 gives a good matching of the experimental data and confirms the differences in the shape of the yield surface between anisotropic and isotropic consolidated samples. The experimental yield points in extension give lower values that those predicted by the model as expected, because no dependency on the Lode's angle, has been introduced in the model.

CONCLUSIONS
An anisotropic model for structured clays (E-SCLAY1S) has been formulated to extend a previous model (S-CLAY1S [12]) for normally or slightly overconsolidated soft clays by introducing the framework of logarithmic contractancy.
In addition to a complete description of the proposed model in triaxial stress space, an implicit Euler backward algorithm for the stress integration has been presented. A major advantage of E-SCLAY1S is that by suitable adjustment of the parameter (n L ), a wide range of yield surface shapes can be achieved. It is important to acknowledge the fact that the proposed model requires only an additional parameter (n L ) and it can be determined from conventional laboratory tests (drained or undrained triaxial tests). As compared with the non-associated flow rule presented by Dafalias and Taiebat [19] for improved predictions, the present model has the advantage of being simpler to calibrate and to implement into a finite element code.
The main features of the model can be summarized as follows: 1. Uniqueness of the CSL at stress space and constant slope of CSL as M, independent of n L valuethis is a major advantage of the proposed model compared with previous anisotropic logarithmic contractancy models (e.g. extended Sekiguchi-Ohta models [16]). 2. K 0 prediction and yield locithe proposed model through the additional parameter (n L ) introduces more flexibility in predicting K 0 and the yield points in p'-q space. 3. Undrained shear strength (c u )an improved prediction of c u can be obtained by adjusting the logarithmic contractancy parameter n L . 4. Stiffnessthe additional parameter (n L ) may be also used to fit the stiffness under shearing stress paths.
The comparisons with laboratory test data of two remoulded clays for different stress ratios and OCRs under undrained shearing revealed the predictive capabilities of the proposed model. The experimental data on Bothkennar clay, Santa Clara clay and Kaolin clay suggest that the model parameter n L , which controls the shape of the yield surface, may not be a soil constant and it can be a hardening parameter that varies with the amount of fabric (degree of anisotropy) of clays. However, further experimental studies on yield points are required to conclude that the parameter n L varies with the degree of anisotropy.
The model verification is limited to soils that do not exhibit bonding and destructuration behaviour with plastic straining. Because triaxial tests were used for validation, additional work should be conducted to verify the model for other stress paths, and also in boundary value problems. Further extension of the model to account for rate-dependent (creep) natural soft soil response can be made along the lines presented by Sivasithamparam et al. [37].