Double image stress point bounding surface model for monotonic and cyclic loading on anisotropic clays

Abstract

In this paper, an anisotropic non-associative constitutive model, based on the well-established critical state soil mechanics framework, is developed for the natural and reconstituted clays subjected to monotonic and cyclic loadings. A robust non-elliptical yield surface is implemented in the model which enables it to capture a wide variety of yield stress points. In addition to a realistic soil stiffness simulation, the employed flexible plastic potential surface equips the proposed model to predict a more accurate coefficient of earth pressure at rest. The model uses a combination of conventional volumetric and comprehensive rotational hardening rules to control the evolution of the yield surface caused by plastic strain increments. The adopted rotational hardening rule uses the concept of governing plastic strain increment which not only ensures the existence of the equilibrium state of anisotropy but also takes the effect of plastic strains at different constant stress ratios into account more sensibly. The proposed model is enhanced with a novel double-image stress point bounding surface plasticity type theory to capture nonlinear behaviour inside the yield surface, and also to simulate the cyclic responses. The detailed model formulation is discussed in both triaxial stress space and the $e-\mathrm{ln}\, p$ plane, and the important features of the model are elaborately explored. The capabilities of the model are also demonstrated by an extensive sensitivity analysis. In the end, the model is used to carry out simulations of monotonic and cyclic element tests on different clays and the results are compared with the available experimental data.

1 Introduction

As the earliest elastoplastic constitutive models based on the critical state soil mechanics (CSSM) theory, the original and modified Cam-clay (OCC and MCC) models [52,53,54] have paved the way to the advanced constitutive modelling of fine-grained cohesive soils. These models enhanced the prediction capability of various experimentally observed responses of isotropically consolidated soft clays. However, some of their intrinsic shortcomings, such as the inflexible and bulky shape yield surface (YS) and the overestimation of peak shear strength in the dry side of the critical state line (CSL), in addition to their inability to simulate the behaviour of anisotropically reconstituted soils with preferred fabric orientation, and also the incapability in capturing the nonlinear hysteretic responses of soils during loading–unloading-reloading cycles, have encouraged researchers to extend these models to more advanced ones that are able to simulate more realistic soil behaviour and also to resolve the abovementioned drawbacks [5, 7, 66].

The initial stress condition which is applied to a representative sample of soil during its deposition history results in an inherent anisotropic material with a preferential arrangement and orientation of soil fabric which, of course, is not constant and can be altered during stress variation and deformations [22, 23]. The anisotropy causes the directional dependent mechanical properties which may contribute to different stress–strain behaviour. Regarding experimental observations [16], the influence of fabric anisotropy has been implemented into models using a sheared/inclined YS which can rotate when the direction of the applied stress is not in line with the orientation of the YS [45,46,47,48, 51, 58]. Most of the existing constitutive models have only addressed the influence of anisotropy on the pre-failure behaviour of soils, while its effect on the failure (i.e., unrestrained shearing at constant stress level and volume which is known as the critical state) has been neglected. Nevertheless, some constitutive models have been developed based on the assumption that the failure is not a function of the anisotropy state, and a unique CSL can be reached independently of both anisotropy and loading direction [61]. Although the uniqueness of the CSL is still a source of debate from the constitutive modelling perspective, the experimental evidence is mainly in favour of the non-unique loading-direction-dependent CSL [68].

To enable elastoplastic models to simulate irreversible strains inside the YS, different theories have been proposed [1, 11, 21]. Among them, bounding surface plasticity (BSP) is perhaps the most recognized theory in the field of soil constitutive modelling that improves the performance of various basic models by enabling them to produce more realistic simulations. This concept relates the plastic modulus of the stress state within the YS to an imaginary stress point on the YS through a set of mapping rule (MR) and projection centre (PC) to allow the plastic strains to develop when the stress state varies inside the YS. The application of this theory in soil plasticity dates back to the works of Dafalias [11], Dafalias and Herrmann [8], and Anandarajah and Dafalias [2]. Since then, the relative simplicity and attractive features of this theory have persuaded researchers to carry out more related studies and developments. Subsequently, various bounding surface models have been proposed (e.g., [20, 26, 27, 35, 39, 41, 49, 50]) to tackle the difficulties associated with the simulation of overconsolidated soils and cyclic responses. This is while finding a suitable set of PC and MR, which are the two essential parts of the BSP theory, remains a challenging aspect. Additionally, simulation of cyclic behaviour is accompanied by mathematical complications related to transferring the PC with changing the loading direction [55] or the need to define more additional parameters [70]. A more detailed discussion on these issues is presented in Rezania and Dejaloud [49].

Using a novel double-image stress point (DISP) method, a robust elastoplastic constitutive model is proposed in this paper to capture the realistic behaviour of a wide range of normally consolidated to highly overconsolidated clays under the application of cyclic and monotonic loadings. The new model which is the second member of the adaptive anisotropic [15] family of models is named AA2-DISP. This method enhances the model to tackle the difficulties associated with finding a suitable set of projection centre and mapping rule to simulate the highly overconsolidated behaviour more accurately. Furthermore, with no need to define new positions for the PC, the proposed DISP method facilitates the model to simulate cyclic behaviour with the least amount of additional parameters. The adopted flexible YS provides the model with the ability to capture the yield stress points precisely both in the dry and wet sides of the CSL. A non-associated flow rule formulation is also integrated by utilising a flexible plastic potential surface (PPS) which satisfies the normality condition (i.e., zero plastic volumetric strain) at the critical state. It should be noted that the term flexible YS/PPS refers to the fact that the proposed constitutive model benefits from a YS equation that, based on the values of two model-specific shape parameters, produces a wide range of YS shapes from a simple ellipse to more complex tear or bullet shape YSs. Moreover, an efficient rotational hardening rule (RH) is adopted which guarantees a unique equilibrium state of anisotropy during plastic shearing which both YS and PPS tend to reach. Finally, the proposed DISP method is implemented to determine the plastic behaviour of clays, when the stress state lies inside the YS, as well as the hysteretic behaviour during loading–unloading–reloading cycles. By establishing a relationship between the plastic moduli of the current stress state and two image stresses on the YS, the proposed DISP method enables the model to recreate promising monotonic and cyclic simulations.

The layout of the paper is as follows; in Sect. 2, the formulation of the model is presented in the triaxial stress space and then it is extended to the $e-\mathrm{ln}\, p$ plane (where e is the void ratio and $p$ is the mean effective stress). An extensive sensitivity analysis is performed in Sect. 3 to determine the contribution of model parameters to the predictions. In the next section, the model performance is verified against the monotonic and cyclic test data related to a number of clays, including Lower Cromer till (LCT) [19], Boston blue clay (BBC) [33], and a compacted clay that was collected from a dam in Thailand [59]. Finally, a concise conclusion is provided in the last section.

2 Model formulation

According to experimental evidence [16, 38, 60], an inclined non-elliptical YS is employed to take the effect of soil anisotropy into account. Furthermore, backed up by the results of different studies [24, 69], in order to improve the stress–strain predictions and also to better capture the coefficient of lateral earth pressure at rest (${K}_{0}$), a non-associated flow rule is used in the model development. In addition to the conventional isotropic hardening law, a robust RH is employed to simulate the development and elimination of anisotropy during loading. Moreover, a new bounding-surface-type concept is introduced to capture both the cyclic response and the nonlinearity of the behaviour of overconsolidated clays inside the YS.

In the following sub-sections, the model formulation in the triaxial stress space is presented together with the rationale behind different model components. For an improved understanding, the model responses are also discussed in the $e-\mathrm{ln}\, p$ plane. It should be noted that all stresses-related parameters used in this paper represent effective stresses, and the prime sign is therefore neglected. Additionally, the mean effective and deviatoric ($q$) stresses are calculated using the effective stress (${\sigma }_{ij}$), deviatoric stress (${s}_{ij}$), and deviatoric fabric (${\alpha }_{ij}^{\text{d}}$) tensors as follows

$$p=\frac{1}{3}{\sigma }_{kk}$$

(1a)

$${q}^{2}=\frac{3}{2}\left({s}_{ij}-{\alpha }_{ij}^{\text{d}}p\right)\left({s}_{ij}-{\alpha }_{ij}^{\text{d}}p\right)$$

(1b)

2.1 Yield surface

Extensive experimental studies on finding the stress points associated with the initiation of plastic flow indicate that a flexible non-elliptical YS best represents the yield loci for different soil types [14, 16]. The flexibility of the YS should be considered an inseparable feature of an accurate soil constitutive model since frequent studies prove that not only YS shapes are different from one soil to another, but also, even for a specific soil, they are highly dependent on the consolidation condition of the soil [15]. To this end, a flexible anisotropic YS is developed in the proposed model which can produce different shapes using two model-specific shape parameters (${n}_{\text{y}}$ and ${r}_{\text{y}}$) in comparison with a typical anisotropic MCC model (e.g., [13, 28, 29, 61]). In the triaxial stress space, the mathematical formulation of the proposed YS reads as:

$$f={\left(q-\alpha p\right)}^{2}-{\left({N}^{2}-{\alpha }^{2}\right){p}^{2}\left(\frac{{C}_{\text{y}}\mathrm{ln}\, \frac{{p}_{0}}{p}}{\mathrm{ln}\,{r}_{\text{y}}}\right)}^{\frac{2}{{n}_{\text{y}}}}=0$$

(2)

where $N$ is a model parameter that represents the stress-ratio-type shape factor of the YS, ${p}_{0}$ is preconsolidation pressure, ${\alpha }^{2}=3/2 {\alpha }_{ij}^{\text{d}}{\alpha }_{ij}^{\text{d}}$, and ${C}_{\text{y}}$ is termed as the coefficient of anisotropy, which for a YS with specific values of $N$ and ${n}_{\text{y}}$, and it depends on its orientation.

$${C}_{\text{y}}={\left(\frac{N-\alpha }{{\left({N}^{2}-{\alpha }^{2}\right)}^\frac{1}{2}}\right)}^{{n}_{\text{y}}}$$

(3)

The isotropic version of this YS, by assuming $\alpha =$ 0 and $N=M$, can be reduced back to the works of Yu and his co-workers [69,70,71] on a unified model for clays and sands, namely CASM. Also, Khalili and his colleagues [30, 31, 40, 56] developed a series of bounding surface constitutive models for clay and sand using an isotropic YS similar to CASM. Due to its flexibility, their proposed YS has turned into a suitable option for clays that are consolidated isotropically; however, it fails dramatically to capture the yield points of natural clays with inherent and/or induced anisotropy. The notable performance of the proposed anisotropic YS in capturing the yield stress points for four different clays is demonstrated in Fig. 1. The experimental data are taken from Díaz-Rodríguez et al. [16] which can be referred to for more details of different clays’ properties. As the figure shows, the capability of the proposed anisotropic YS equation in capturing different formations of yield stress points confirms the importance of implementing both flexibility and anisotropy into the YS formulation.

Fig. 1

Comparison of the proposed yield surface with experimental yield stress points for six different clays: a Winnipeg clay, b Atchafalaya clay, c Riihimäki, Saint-Louis, and Ottawa clays, and d Mexico clay. Data from Diaz Rodriguez et al. [16]

Neglecting $N$, which solely adjusts the bulkiness and volume of the YS, ${n}_{\text{y}}$ and ${r}_{\text{y}}$ are two model parameters that specify the shape of the YS. A well-balanced combination of these parameters enables the model to effectively capture the yield stress points. The flexibility postulated by the implementation of these parameters results in significant improvements in simulating clay behaviour in both normally consolidated and overconsolidated states. Figure 2 shows how the variations of YS shape parameters contribute to transforming the configuration of the surface. As the figure shows, while ${n}_{\text{y}}$ alters the overall shape of the YS, ${r}_{\text{y}}$ can be treated as a measure that changes the volume of the YS. As can be seen, ${r}_{\text{y}}$ affect the intersection point of the YS and $N$-line which for the case of $N=M$ it is known as the spacing ratio [64]. On the other hand, the intersection point remains fixed by the variation of ${n}_{\text{y}}$, which means that the spacing ratio is independent of ${n}_{\text{y}}$ only when $N=M$. Simple models with a fixed-shape YS consider a constant value for spacing ratio, while experimental observations promote that different spacing ratios should be adopted for a realistic representation of yield points [6]. Based on the existing evidence, not only the spacing ratio varies for different types of clays, but also it is not fixed for a specific soil with different overconsolidation ratio (OCR) values. Some experimental data are suggesting that the consolidation condition can also greatly influence the shape of the YS, which in turn can lead to different spacing ratios [15]. These points confirm the deficiencies associated with adopting a fixed-shape YS in soil constitutive models as they are unable to recreate adequate simulations over a range of different initial conditions.

Fig. 2

Variation of the proposed YS with its shape parameters: a ${n}_{\text{y}}$ and b ${r}_{\text{y}}$

As an additional point, it should be noted that in the isotropic state ($\alpha =0$), by setting $N=M$, ${n}_{\text{y}} = 1$ and ${r}_{\text{y}}\cong 2.718$, the YS of the proposed model is reduced to that of the OCC model.

2.2 Plastic potential and flow rule

Some previous studies suggest that using an associated flow rule by adopting identical YS and PPS is a reasonable assumption when the soil anisotropy is considered with an inclined surface (e.g., [44, 61]), while a number of works concluded later that in some cases associated flow rule may not lead to favourable predictions (e.g., [9, 25]). Adopting a non-associated flow rule can significantly improve the capability of the model in simulating the peak shear strength before reaching CSL, which is an important feature of the undrained normally consolidated natural clays that poses inherent anisotropy [69]. Additionally, implementing a flexible PPS ensures that the model can readily simulate ${K}_{0}$ conditions [18]. To this end, in the proposed model, a non-associated flow rule is adopted by introducing a flexible PPS, which is oriented in accordance with the YS. The adopted PPS can be defined in the triaxial stress space using the following expression:

$$g={\left(q-\alpha p\right)}^{2}-{C}_{\text{p}}\left({M}^{2}-{\alpha }^{2}\right){p}^{2}{\left(\frac{1}{{m}_{\text{p}}-1}\left[{\left(\frac{{p}_{g}}{p}\right)}^{\frac{{n}_{p}\left({m}_{\text{p}}-1\right)}{{m}_{\text{p}}}}-1\right]\right)}^{\frac{2}{{n}_{p}}}=0$$

(4)

where ${n}_{p}$ and ${m}_{\text{p}}$ are the model parameters that control the shape of the plastic potential surface (i.e., plastic potential shape parameters), $M$ is the slope of the CSL in stress space, and ${C}_{\text{p}}$ is an anisotropic coefficient that can be determined as:

$${C}_{\text{p}}=\frac{{\left({m}_{\text{p}}M{\left(M-\alpha \right)}^{{n}_{p}-1}-\left({m}_{\text{p}}-1\right){\left(M-\alpha \right)}^{{n}_{p}}\right)}^{\frac{2}{{n}_{p}}}}{{M}^{2}-{\alpha }^{2}}$$

(5)

It should be mentioned that the employed YS does not satisfy the zero plastic volumetric strain increment at the CSL (which is a prerequisite for plastic flow at the critical state). Hence, to overcome this limitation, a different surface is required as the PPS.

To calculate the gradient of the PPS at the stress state on the YS, it is necessary to find the size of the surface, ${p}_{g}$, by solving Eq. 4. Figure 3 shows how parameters ${n}_{p}$ and ${m}_{\text{p}}$ change the shape of the PPS. Similar to the YS shape parameters, the effect of these parameters on the model responses will be discussed in detail in the subsequent section on model parameters.

Fig. 3

Variation of the proposed PPS with its shape parameters: a ${n}_{p}$ and b ${m}_{\text{p}}$

As the PPS is defined, now the plastic volumetric ($d{\varepsilon }_{\text{v}}^{\text{p}}$) and deviatoric ($d{\varepsilon }_{\text{d}}^{\text{p}}$) strain increments can be readily determined using the following expressions:

$$d{\varepsilon }_{\text{v}}^{\text{p}}=<L>\frac{\partial g}{\partial p}$$

(6a)

$$d{\varepsilon }_{\text{d}}^{\text{p}}=<L>\frac{\partial g}{\partial q}$$

(6b)

where $< >$ represent the Macaulay brackets and $L$ is the so-called loading index.

Moreover, an investigation of the stress-dilatancy relation ($d=d{\varepsilon }_{\text{d}}^{\text{p}}/d{\varepsilon }_{\text{v}}^{\text{p}}$) helps to understand how the adopted PPS contributes to flexible simulations under the one-dimensional compression loading. By derivating the PPS and applying conditions associated with the material at ${K}_{0}$ state (${\eta }_{{K}_{0}}$ and ${\alpha }_{{K}_{0}}$), the dilatancy can be read as follows:

$$d=\frac{{m}_{\text{p}}{\left({\eta }_{{K}_{0}} -{\alpha }_{{K}_{0}}\right)}^{{n}_{p}-1}}{\left({m}_{\text{p}}-1\right){\left({\eta }_{{K}_{0}} -{\alpha }_{{K}_{0}}\right)}^{{n}_{p}}+{m}_{\text{p}}M{\left(M-{\alpha }_{{K}_{0}}\right)}^{{n}_{p}-1}-\left({m}_{\text{p}}-1\right){\left(M-{\alpha }_{{K}_{0}}\right)}^{{n}_{p}}-{m}_{\text{p}}{\eta }_{{K}_{0}}{\left({\eta }_{{K}_{0}} -{\alpha }_{{K}_{0}}\right)}^{{n}_{p}-1}}$$

(7)

where ${\alpha }_{{K}_{0}}$ refers to the inclination of the YS when the soil element is subjected to the constant ${\eta }_{{K}_{0}}=3(1-{K}_{0})/(1+2{K}_{0})$ loading condition. Neglecting the elastic strains with regard to plastic ones, the zero lateral strain condition during one-dimensional compression results in a constant ratio between plastic deviatoric and volumetric strain increments which is $d=2/3$. By equating this value with Eq. (7), the variation of the predicted ${\eta }_{{K}_{0}}$ with regard to the model parameters can be easily investigated (as is done in the model parameter section).

It should be mentioned that to consider the effect of Lode angle dependency on the variation of $M$ and $N$ between their corresponding values in compression (${M}_{\text{c}}$ and ${N}_{\text{c}}$) and extension (${M}_{\text{e}}$ and ${N}_{\text{e}}$), the proposed expression by Sheng et al. [57] is applied to YS and PPS in a similar manner

$$({M}^{2}-{\alpha }^{2})={({M}_{\text{c}}^{2}-{\alpha }^{2})\left(\frac{2{r}^{4}}{1+{r}^{4}-\left(1-{r}^{4}\right)\mathrm{sin}3\theta }\right)}^{1/4}$$

(8)

where $r=({M}_{\text{e}}^{2}-{\alpha }^{2})/({M}_{\text{c}}^{2}-{\alpha }^{2})$. Moreover, the Lode angle is defined by $\mathrm{sin}3\theta =3\sqrt{3}/2{\left(S/J\right)}^{3}$, where $J$ and $S$ are the second and third stress invariants, respectively. The Lode-angle-dependent shape of the CSL and YS, for both isotropic and anisotropic states, in the deviatoric plane of the general stress space is shown schematically in Fig. 4.

Fig. 4

CSL and YS in the deviatoric plane of the general stress space

2.3 Hardening rules

The proposed model is enhanced with two different hardening rules, namely isotropic and RH rules, that change the size and the orientation of the YS with plastic volumetric and plastic shear strains.

2.3.1 Isotropic hardening

A conventional volumetric hardening rule is employed to control the expansion and contraction of the YS size. The isotropic hardening rule is defined in the same manner as in the MCC model, where the changes in the size of the YS are related to the plastic volumetric strain increments as:

$$d{p}_{0}=\frac{v{ p}_{0}}{\lambda -\kappa } d{\varepsilon }_{\text{v}}^{\text{p}}$$

(9)

where $v=1+e$ is the specific volume, and $\lambda$ and $\kappa$ are the slopes of the normal compression and swelling/unloading lines in the $e-\mathrm{ln}\,p$ plane.

2.3.2 Rotational hardening

The proposed constitutive model is categorised as an anisotropic model where the plastic anisotropy is incorporated through the inclination of the YS, and the plastic strains-induced development and erasure of anisotropy are controlled through a RH rule. The adopted RH rule should pose some key features (i.e., in terms of governing plastic strain and equilibrium anisotropy state) to produce reasonable simulations. In the following, the details of how these aspects form the backbone of the proposed RH rule are explained.

2.3.2.1 Governing plastic strain increment

Governing plastic strain increments refer to how the components of the plastic strains contribute to evolve the anisotropy during a desirable loading condition. A diverse combination of the plastic strain increment constituents is believed to contribute to soil anisotropy. For example, some studies relate the rotation of the YS to only the plastic volumetric strain increment (e.g., [3, 12, 13, 25, 36, 62]) or only the plastic deviatoric strain increment (e.g., [7]), while other studies believe in the contribution of both plastic volumetric and deviatoric strain increments in the evolution of anisotropy (e.g., [42, 44, 61]). Despite these discrepancies among researchers, there is almost no convincing answer to the questions about “how the contribution from different components of the plastic strain increments varies during loading?” and “how to implement this variation into the rotational hardening rule?”. In this study, the proposed hardening rule by Dejaloud and Rezania [15] is reformulated to address these questions by adopting a governing strain increment that varies when the stress state tends to reach to or depart from the CSL.

Both plastic volumetric and deviatoric strain increments are contributing to the proposed governing strain increment ($d{\varepsilon }_{\text{g}}^{\text{p}}$) to control the rotation of the YS during the shearings. Using a hyperbolic function, the dominance of plastic volumetric strain at low values of stress ratio ($\eta =q/p$) is satisfied, while this function signifies the role of plastic deviatoric strain simultaneously with the increment of stress ratio. In other words, for isotropic compression condition ($\eta =0$), the governing strain increment turns to the plastic volumetric strain increment ($d{\varepsilon }_{\text{g}}^{\text{p}}=d{\varepsilon }_{\text{v}}^{\text{p}}$), and when the stress state touches the CSL ($\eta =M$), the governing strain increment changes to the plastic deviatoric strain increment ($d{\varepsilon }_{\text{g}}^{\text{p}}=d{\varepsilon }_{\text{d}}^{\text{p}}$). To this end, the proposed governing strain increment reads as

$$d{\varepsilon }_{\text{g}}^{\text{p}}=\left(\mathrm{A }d{\varepsilon }_{\text{v}}^{\text{p}}+\left[1-\mathrm{A}\right] |d{\varepsilon }_{\text{d}}^{\text{p}}|\right)$$

(10)

where $A$ is a hyperbolic function that is expressed as

$$A=\mathrm{tanh}\left(a<1-\frac{\left|\eta \right|}{M}{>}^{b}\right)$$

(11)

where $a$ and $b$ are the parameters that control the contribution of the strain increment components. Although calibrated values can be assigned to these parameters, for practicality, constant values of $a=$ 5 and $b=$ 2 are proposed based on the authors’ experience, which provides a smooth transition between two components of plastic strain increments.

2.3.2.2 Equilibrium state of anisotropy

Another important aspect of clay anisotropy that should be taken into account by a RH rule is the concept of the equilibrium state of anisotropy as an ultimate target that attracts the YS in the stress space during non-hydrostatic loading. This bounding value also acts as a limiting constraint controlling the excessive rotation of the YS. In the physical perspective, the equilibrium state represents a specific configuration of soil fabric with a preferable inter-particle structure under the application of a constant-stress-ratio loading that affects the shear strength and stress–strain behaviour. Dafalias [12] considered a fraction of $\eta$ as the bounding value by employing a model constant $x$ (Eq. 12 in the reference). Although this simplified expression produces satisfactory results, it might be unable to restrict the rotation of the YS when $\eta >M$. To mitigate the excessive rotation of the YS, Dafalias and Taiebat [9] redefined parameter $x$ as a function of $\eta$ through an exponential function and illustrated that the newly modified version provides more conformity with the experimental observations. Dafalias and Taiebat [10] later improved their last hypothesis and introduced an expression to deal with the isotropic fabric orientation during critical state failure. It was for the equilibrium state for the range of $|\eta |\ge M$ with a declining trend for bounding value to zero as the quantitative representation of the isotropic state. Wheeler et al. [61] defined the equilibrium state of fabric anisotropy implicitly through their RH rule which can be attained by solving $\dot{\alpha }=0$. In consequence, in their model, the equilibrium value of anisotropy varies between $1/3 \eta$ and $3/4\eta$, according to the value of a model constant $\beta$ (Fig. 5 in the reference). However, their proposed RH rule does not guarantee the limit of the YS evolution, since there are loading conditions (i.e., $\eta >3M$) that might lead to unreasonable rotations.

Fig. 5

Equilibrium state of anisotropy for Otaniemi clay for constant stress ratio loading. Data from Wheeler et al. [61]

Based on the contribution of plastic strain components, an equilibrium state of anisotropy is defined in this model using hyperbolic function $A$. The proposed equilibrium state varies between two limiting values (${\chi }_{v}$ and ${\chi }_{d}$) which are model parameters and represent the orientation that the YS tries to reach due to pure plastic volumetric and deviatoric strains under constant stress ratio loading condition when $\eta =M$. It can be concluded when the stress state reaches the CSL (where all plastic strains are deviatoric, i.e., $d{\varepsilon }_{\text{g}}^{\text{p}}=d{\varepsilon }_{\text{d}}^{\text{p}}$), the YS rotates along with ${\chi }_{d}$. The proposed equilibrium state of anisotropy can be expressed as:

$${\alpha }_{\text{e}}=\eta \left[A \left({\chi }_{v}-{\chi }_{d}\right)+{\chi }_{d}\mathrm{exp}\left(-c<\frac{\left|\eta \right|}{M}-1>\right)\right]$$

(12)

where $c$ is a model parameter. Since the first term in the square brackets determines the equilibrium state of anisotropy before touching the CSL (i.e., $\eta \le M$), the second term restrains the rotation of the YS when the stress ratio passes the CSL ($\eta >M$) and prevents the possibility of excessive rotation. This is while the RH of the S-CLAY1 model [61] does not restrict the rotation of the YS which might lead to numerical errors (i.e., analysis termination) when $\alpha \to M$, because the elliptical YS turns into a line. Based on the authors’ experience, adopting $c>$ 1 ensures that not only the rotation of the YS is restricted, but also it leads to satisfactory simulations. In addition, this aspect of the proposed RH rule addresses one of the challenging topics about the isotropy/anisotropy state of the soil when it reaches the CSL [10].

Based on experimental observations, and backed by numerical inspections, ${\chi }_{v}$ is considered to be equal to unity (${\chi }_{v}=1$). By choosing ${\chi }_{d}$ as an ultimate target that the YS tries to incline along with when the stress state reaches the CSL, the variation of the equilibrium state of anisotropy for different stress ratios is determined simultaneously with governing plastic strain increment by function $A$. In other words, function $A$ creates a correlation between the governing plastic strain and the equilibrium state based on the stress ratio level. Figure 5 illustrates the comparison of actual and simulated values of the equilibrium state of anisotropy for Otaniemi clay [61]. By considering the calibrated values of ${\chi }_{d} = 0.4$, $a=5.5$ and $b=3.5$ to fit with experimental data at low and high stress ratios, the suitability of Eq. (12) in reproducing experimental data is obvious.

The set of governing plastic strain increment and equilibrium state of anisotropy, plus a simple expression that controls the absolute rate of the YS rotation, form the rotational hardening rule of the proposed model, as:

$$d\alpha =\mu \frac{p}{{p}_{0}} \left({\alpha }_{\text{e}}-\alpha \right)d{\varepsilon }_{\text{g}}^{\text{p}}$$

(13)

where $\mu$ is a model parameter that controls the absolute rate of the YS rotation tending to reach the equilibrium state, and $p/{p}_{0}$ is considered to decrease the rate of the YS rotation for over-consolidated conditions.

To achieve accurate model simulations, it is crucial to properly define the initial state of the soil that in addition to the values of preconsolidation pressure and initial void ratio, it should include the initial inclination of the YS, representing the inherent anisotropy of soil fabric (i.e., ${\alpha }_{0}$). During the sedimentation process of natural clays and also through the preparation of reconstituted samples, soil fabric orientation alters to lie alongside the anisotropy equilibrium state corresponding to a specific stress ratio ${\eta }_{0}$. However, due to the limited applied strain during the above-mentioned procedures, the inclination might not fully reach the equilibrium state. Therefore, the domain between ${\chi }_{d}$-line and ${\alpha }_{\text{e}}$ is assumed as the permissible range of ${\alpha }_{0}$ and defines it as a simple function of ${\eta }_{0}$ by

$${\alpha }_{0}={\omega \eta }_{0}$$

(14)

where for simplicity $\omega$ is defined as the mean of ${\chi }_{d}$ and ${\alpha }_{\text{e}}$

$$\omega =\frac{1}{2}\left({\chi }_{d}+\mathrm{tanh}\left(a<1-\frac{|{\eta }_{0}|}{M}{>}^{b}\right) \left(1-{\chi }_{d}\right)\right)$$

(15)

2.4 Plastic modulus on the YS

In the case of establishing a complete stress–strain relationship for normally consolidated states, the consistency condition should be satisfied, which guarantees that the stress state does not cross the YS and remains on it during the plastic flow ($\dot{f}=0$). By determining the loading index,$L$, from the routine plasticity procedure as

$$L=\frac{1}{{K}_{\text{p}}}\frac{\partial f}{\partial {\sigma }_{ij}}d{\sigma }_{ij}$$

(16)

and following a standard differentiation of the YS with regard to its state variables and substituting the hardening rules and loading index in the resultant differential equation, the plastic modulus, ${K}_{\text{p}}$, is obtained as:

$${K}_{\text{p}}=-\frac{\partial f}{\partial \alpha }\left(\frac{\partial \alpha }{\partial {\varepsilon }_{\text{v}}^{\text{p}}}d{\varepsilon }_{\text{v}}^{\text{p}}+\frac{\partial \alpha }{\partial {\varepsilon }_{\text{d}}^{\text{p}}}|d{\varepsilon }_{\text{d}}^{\text{p}}|\right)-\frac{\partial f}{\partial {p}_{0}}\frac{\partial {p}_{0}}{\partial {\varepsilon }_{\text{v}}^{\text{p}}}d{\varepsilon }_{\text{v}}^{\text{p}}$$

(17)

2.5 Plastic modulus inside the YS–DISP method

In this paper, a novel and efficient method, named the double-image stress point (DISP), is proposed to simulate the nonlinear behaviour associated with the shearing of overconsolidated clays, as well as reproducing the cyclic responses. This method can be interpreted as an enhanced extension of the BSP where it tries to relate the plastic modulus of the current/actual stress state to the plastic moduli of two imaginary points on the YS through a set of PC and MR using a weighted average technique. The fundamental assumption in this approach is that the stress state inside the YS inherits its specifications from two points by considering the relative distances between them and the actual stress state, which is different from the conventional BSP method that obtains the plastic modulus from only one imaginary stress state on the YS. In addition to better monotonic simulations, this assumption enables the model to simulate the loading–unloading–reloading cyclic capacity of soils with good accuracy, using a minimum number of additional parameters. By neglecting the additional surfaces in its formulation (i.e., loading surface in the BSP theory), the double-image stress point method cannot be categorised into the multi-YS methods. However, in a special condition that will be explained in the sequel, this method can turn into a conventional BSP model.

Figure 6 demonstrates the concept behind the proposed DISP method. The radial mapping rule passing from the projection centre on the $\alpha$-line and the actual stress state crosses the YS on two points. Using a weighted average and a shape hardening function (${S}_{l}$), the plastic modulus of the stress state inside the YS (i.e., ${\overline{K} }_{p}$) is related to the plastic moduli of the two image stress points

$${\overline{K} }_{p}={\phi }_{1}^{2}{K}_{\text{p}}^{I{M}_{1}}+{\phi }_{2}^{2}{K}_{\text{p}}^{I{M}_{2}}+{S}_{l}$$

(18a)

$${\phi }_{1}=\frac{\left|I{M}_{2}-SS\right|}{\left|I{M}_{2}-I{M}_{1}\right|}$$

(18b)

$${\phi }_{2}=\frac{\left|I{M}_{1}-SS\right|}{\left|I{M}_{2}-I{M}_{1}\right|}$$

(18c)

where the absolute sign in the above equations indicates the distance between the points. The shape hardening function is defined using the following equation:

$${S}_{l}=h {p}^{3}\mathrm{cot}\left(\frac{\pi }{2} {\phi }_{1}^{2\sqrt{h}}\right)$$

(19)

where $h$ is a model parameter that controls the relationship between the plastic moduli of the actual stress state and the image stress points; its influences on the overall response of the model for both overconsolidated and cyclic conditions will be illustrated in the model parameters section. As is clear, when the stress state reaches the $I{M}_{1}$, then ${\phi }_{1}=1$, ${\phi }_{2}=0$ and the value of the shape hardening function tends to zero (${S}_{l}=0$), which means ${\overline{K} }_{p}={{K}_{\text{p}}}_{1}^{IM}$. Moreover, it should be noted that the definitions of $I{M}_{1}$ and $I{M}_{2}$ depend on both loading direction and the position of stress state with respect to the PC. Figure 7 illustrates the steps of defining image points during a cyclic loading procedure for a simple isotropically consolidated sample. In the loading procedure ($L>0$), the image stress point near the actual stress state is termed as $I{M}_{1}$ (Fig. 7a, c), while during the unloading ($L<0$) $I{M}_{1}$ jumps to the opposite side (Fig. 7b). It should be mentioned that gradients of YS and PPS associated with point $I{M}_{1}$ are utilised in the model calculations, and $I{M}_{2}$ just participates in calculations of plastic modulus.

Fig. 6

Schematic of adopted projection centre and mapping rule in the proposed DISP method

Fig. 7

Steps of defining image stress points during the application of cyclic loading stages, a loading, b unloading, and c reloading

To find the image stress points on the YS, based on the experimental evidence [32], a moving PC is adopted in this model [49]. The rationale behind the PC is that during the loading procedure, the $I{M}_{1}$ related to the stress states on the subcritical area should be projected to the subcritical side of the YS, and the one related to the supercritical stress states should be projected on the supercritical side of the YS. To this end, the projection centre is defined in such a way that it moves freely on the $\alpha$-line with the variation of the stress state. Therefore, the location of the PC (${\sigma }_{ij}^{\text{c}}$) can be defined as a fraction of the size of the YS, as

$${\sigma }_{ij}^{\text{c}}=\gamma \left({p}_{0}-p\right)\left({\alpha }_{ij}^{\text{d}}+{\delta }_{ij}\right)$$

(20)

In this equation, ${\delta }_{ij}$ is the Kronecker delta, and $\gamma$ controls the location and the pace of evolution of the projection centre and can be read as:

$$\gamma =\frac{p}{{p}_{0}}\frac{{{r}_{\text{y}}^{\left(\frac{1-{\chi }_{d}}{N/M-{\chi }_{d}}\right)}}^{{n}_{\text{y}}}}{{{r}_{\text{y}}^{\left(\frac{{\eta }_{0}-{\alpha }_{0}}{N-{\alpha }_{0}}\right)}}^{{n}_{\text{y}}}}$$

(21)

It should be considered that the PC must always lay inside the YS. Therefore, the parameter $\gamma$ should be limited to the maximum value

$${\gamma }_{max}=\frac{{p}_{0}}{{p}_{0}-p}$$

(22)

Furthermore, if $\gamma =0$, the projection centre sticks to the origin of the stress space and the double-image stress point method degrades to the conventional bounding surface method. It is to be noted that in the case of normally consolidated clay (OCR = 1), where the stress state lies on the yield surface, the model response is only governed by the conventional elastoplasticity which means the DISP method is not involved. Similarly, during dilation as the stress state lies on the YS (due to the plastic nature of dilation deformations), the DISP method is also not involved. In other words, the proposed DISP method only plays its role when the stress state moves inside the YS.

2.6 Normal consolidation and critical state lines (NCL and CSL)

To provide a better understanding of the model, it is necessary to elaborate on its formulation in both stress space and $e-\mathrm{ln}\, p$ plane. For example, reaching the CSL in the stress space (i.e., $\eta =M$) is not sufficient for experiencing the critical state condition (i.e., $\dot{p}=$ 0, $\dot{q}=$ 0, $d{\varepsilon }_{\text{v}}=$ 0, but $d{\varepsilon }_{d}\ne$ 0), since the classical critical state theory believes that it is necessary for the void ratio of the sample to also reach to the critical state void ratio simultaneously (i.e., $e={e}_{\text{CSL}}$) in the $e-\mathrm{ln}\, p$ space. However, for granular soils it is accepted that an additional condition related to the fabric orientation and anisotropy should be met in order to reach the critical state condition [34].

In the proposed model, NCL and CSL are represented by two parallel lines in the $e-\mathrm{ln}\, p$ plane, which in the case of clays they can be defined using logarithmic expressions as:

$${e}_{\text{NCL}}={e}_{\text{N}}-\lambda\, \mathrm{ln}\, p$$

(23a)

$${e}_{\text{CSL}}={e}_{\Gamma }-\lambda\, \mathrm{ln}\, p$$

(23b)

where ${e}_{\text{N}}$ and ${e}_{\Gamma }$ are the void ratio of the normally consolidated clay and the critical state void ratio at $p=1$ kPa, and as mentioned before, $\lambda$ is the slope of the NCL and CSL lines in the $e-\mathrm{ln}\, p$ space (Fig. 8). The NCL line represents the loci of material states of normally consolidated soils under the application of isotropic stress. However, in the general case of anisotropically normally consolidated samples, the material states deviate from the NCL and form a parallel anisotropic compression line (ACL) in the $e-\mathrm{ln}\, p$ plane, that can be expressed, similar to the NCL and CSL, as:

$${e}_{\text{ACL}}={e}_{\alpha }-\lambda\, \mathrm{ln}\, p$$

(24)

where ${e}_{\alpha }$ represents the void ratio associated with the ACL at $p=1$ kPa, and is related to both consolidation stress ratio (${\eta }_{0}$) and the YS configuration (i.e., shape and inclination), and it can be expressed with the following expression:

$${e}_{\alpha }={e}_{\text{N}}-(\lambda -\kappa)\, \mathit{ln}{{r}_{\text{y}}^{\left(\frac{{\eta }_{0}-{\alpha }_{0}}{N-{\alpha }_{0}}\right)}}^{{n}_{\text{y}}}$$

(25)

where ${e}_{\alpha }={e}_{\text{N}}$ when the sample is consolidated under isotropic condition (${\eta }_{0}={\alpha }_{0}=0$).

Fig. 8

NCL, ACL and CSL lines in $e-\mathrm{ln}\, p$ plane

The distance between the current void ratio and the void ratio at the CSL in the same mean effective stress was first defined as the state parameter ($\psi$) by Been and Jefferies [4]. Since the state parameter distinguishes the overconsolidated clays or dense sands ($\psi <0$) from lightly overconsolidated clays or loose sands ($\psi >0$) straightforwardly, several constitutive models took advantage of this parameter to adapt their models for different material states (e.g., [37, 43, 66, 67, 69]). The state parameter for a material state is a function of basic soil parameters and the current stress state as:

$$\psi =e-{e}_{\text{CSL}}=\left(\lambda -\kappa \right)\mathrm{ln}\, p/{p}_{c}$$

(26)

where ${p}_{c}$ is the mean effective stress at the intersection point of the unloading line and the CSL (Fig. 8). By extending this definition to the vertical distance between the ACL (as the general form of NCL) and CSL in the $e-\mathrm{ln}\, p$ space, the reference state parameter (${\psi }_{\text{R}}$) can be formulated as:

$${\psi }_{\text{R}}={e}_{\alpha }-{e}_{\Gamma }=\left(\lambda -\kappa \right)\mathrm{ln}\, {p}_{\alpha }/{p}_{c}$$

(27)

where ${p}_{\alpha }$ represents the consolidation pressure on the YS where $\eta ={\eta }_{0}$, which is identical to the preconsolidation pressure (${p}_{0}$) for ${\eta }_{0}={\alpha }_{0}=0$. The ratio between the consolidation pressure and the effective mean stress at the CSL is called spacing ratio ($R={p}_{\alpha }/{p}_{c}$), and it plays an important role in calibrating the position of the CSL in $e-\mathrm{ln}\, p$ plane. In constitutive models with a flexible YS, the spacing ratio plays a vital role and can be considered as their cornerstone. In other words, flexible YSs are developed to address different spacing ratios, in addition to providing a realistic representation of the yield stress points. In this work, based on the developed YS equation of the proposed model, the spacing ratio can be determined as:

$$R=\frac{{{r}_{\text{y}}^{\left(\frac{M-{\alpha }_{\text{e}}}{N-{\alpha }_{\text{e}}}\right)}}^{{n}_{\text{y}}}}{{{r}_{\text{y}}^{\left(\frac{{\eta }_{0}-{\alpha }_{0}}{N-{\alpha }_{0}}\right)}}^{{n}_{\text{y}}}}$$

(28)

where ${\alpha }_{\text{e}}={\chi }_{d}M$ is the equilibrium state of anisotropy that the YS reaches when the stress state meets the critical state. Regarding this equation, one can easily conclude that the CSL is moving during loading procedure and its position can vary with the evolution of the fabric anisotropy (i.e., by variation of the YS inclination), until the YS reaches the equilibrium state of anisotropy at ${\alpha }_{\text{e}}$. This might seem false at first glance because it can lead to a non-unique CSL during strainings under a certain loading conditions. However, by ensuring that the YS always reaches an equilibrium state inclination at the CSL, this variation does not violate the CSL uniqueness. This is while the convergence of the YS inclination to an equilibrium state of anisotropy might not always be satisfied with some RH rules (i.e., [9, 14]) that can lead to inconsistencies between model responses in the stress space and the $e-\mathrm{ln}\, p$ plane. In other words, the uniqueness of the CSL is tied closely to the uniqueness of the equilibrium state of anisotropy. Regarding this, in the case of considering anisotropy, it seems crucial to add the term of ${\alpha }_{\text{CSL}}={\alpha }_{\text{e}}$ to the conditions associated with the critical state, which means at the critical state, all the following conditions should be satisfied:

$$\eta = \eta_{\text{CSL}} = M,e = e_{\text{CSL}} , \alpha = \alpha_{\text{CSL}} = \alpha_{\text{e}}$$

(29)

Therefore, if the RH rule is unable to guarantee a unique equilibrium state of anisotropy, it certainly cannot simulate a unique CSL. The uniqueness of the equilibrium state of anisotropy for the proposed model can be proved straightforwardly. Considering Eq. (13), the YS stops rotating when the stress state reaches CSL as soon as either ${\alpha }_{\text{e}}=\alpha$ or $d{\varepsilon }_{\text{g}}^{\text{p}}=0$. As the governing plastic strain increment includes the deviatoric strain increment ($d{\varepsilon }_{\text{d}}^{\text{p}}$), it does not tend to zero at CSL ($\partial g/\partial q \ne$ 0). Therefore, the YS rotation continues until it aligns in the direction of the equilibrium state of anisotropy.

Figure 9 shows the variation of the spacing ratio with regard to all its influential parameters. Since ${n}_{\text{y}}$ and ${r}_{\text{y}}$ change the spacing ratio almost linearly, increment of $N/M$ ratio decreases the $R$ value dramatically. It means that high values of $N/M$ decreases the distance between ACL and CSL, and a hypothetical extreme state of $N/M\to \infty$ leads to $R=1$ which makes the ACL and CSL coincident. Besides the YS shape parameters and $N/M$ ratio, ${\chi }_{d}$ can alter the position of both ACL and CSL in the $e-\mathrm{ln}\, p$ plane by changing the initial (${\alpha }_{0}$) and final inclination (${\alpha }_{\text{e}}$) of the YS. As can be seen, for high values of ${\chi }_{d}$, which means more anisotropic material at both initial condition and CSL, the distance between ACL and CSL soars significantly. In simple words, since the material tends to reach a more anisotropic state, it takes more time for fibres to align in the direction that is represented by ${\alpha }_{\text{e}}$. Therefore, this process is associated with a more reduction in the mean effective stress. It might be the answer to the question that “how fabric anisotropy could affect the critical state?” which is not explicitly addressed in the existing constitutive models, while in most cases only the contribution of anisotropy to the pre-failure behaviour of clays has been demonstrated.

Fig. 9

Variation of spacing ratio with regard to its parameters

By combining Eqs. (25b), (27) and (28), the critical state line can be determined using the following equation:

$${e}_{\text{CSL}}={e}_{\text{N}}-\left(\lambda -\kappa \right)\mathrm{ln}\, R-\lambda \mathrm{ln}\, p$$

(30)

For a unique CSL in compression and extension loading directions, the equilibrium state of anisotropy should be considered as a function of the Lode angle (i.e., ${\alpha }_{\text{e}}={\alpha }_{\text{e}}(\theta )$) in the same way that $N$ and $M$ have been defined using Eq. (8). In other words, they should share the same Lode-angle-dependent factor, which in the triaxial space it leads to

$$\frac{{N}_{\text{e}}}{{N}_{\text{c}}}=\frac{{M}_{\text{e}}}{{M}_{\text{c}}}=\frac{{\alpha }_{\text{e}}^{\text{e}}}{{\alpha }_{\text{e}}^{\text{c}}}$$

(31)

where ${\alpha }_{\text{e}}^{\text{e}}$ and ${\alpha }_{\text{e}}^{\text{c}}$ are the equilibrium state of anisotropy corresponding to extension and compression loadings, respectively. To investigate this condition, an isotropic sample can be assumed subjected to both compression and extension loadings. Figure 10 compares the model simulations for different ratios of ${\alpha }_{\text{e}}^{\text{e}}/{\alpha }_{\text{e}}^{\text{c}}$. As can be seen, by adopting Eq. (31), the CSLs related to compression and extension become identical, and in both cases, the stress paths reach the critical state at an equal mean effective stress value, as shown by the red arrows in Fig. 10a. This is while considering more (i.e., ${\alpha }_{\text{e}}^{\text{e}}/{\alpha }_{\text{e}}^{\text{c}}>{M}_{\text{e}}/{M}_{\text{c}}$) or less (i.e., ${\alpha }_{\text{e}}^{\text{e}}/{\alpha }_{\text{e}}^{\text{c}}<{M}_{\text{e}}/{M}_{\text{c}}$) anisotropic states for extension shift the corresponding CSL further or closer to NCL, respectively, and leads to different mean effective stress values at CSL with respect to the compression.

Fig. 10

Comparison of the undrained stress paths in compression and extension by considering the different values of the equilibrium state of anisotropy ratios

3 Model parameters

To provide a more clear understanding basis about the contribution of different components of the proposed model, a detailed evaluation and calibration of model parameters is conducted through a consecutive process. The investigation of each category of model parameters (i.e., YS and PPS shape parameters, RH parameters and DISP parameters) is associated with a sensitivity analysis to clarify the influence of a broad range of parameter values on the model responses.

The proposed model has 8 additional parameters compared to the MCC model. These parameters can be determined using the results of conventional undrained triaxial and one-dimensional loading tests. Table 1 shows a brief definition of model parameters and the most convenient tests to determine their values. The calibration procedure is commenced with the tuning of parameters that control the shape of the YS and PPS using the results of undrained triaxial test on the normally consolidated samples and also one-dimensional loading tests. Calibration of the RH rule parameters, using either one-dimensional loading or undrained triaxial tests, comes after the optimization of the YS and PPS. The DISP parameter is the last one that should be calibrated using undrained triaxial test data on the overconsolidated samples. In the following subsections, more details about each model parameter are provided. It should be noted that the sensitivity analysis is carried out on the normally ($\psi /{\psi }_{\text{R}}=1$), lightly ($\psi /{\psi }_{\text{R}}=0.5$), and highly overconsolidated ($\psi /{\psi }_{\text{R}}=-0.5$) samples that are consolidated anisotropically (${K}_{0}=0.5$) in $p=200$ kPa. Also, model responses are compared in both stress path and stress–strain spaces as well as in the $e-\mathrm{ln}\, p$ plane. It is noteworthy that the conventional elastic and critical state parameters are considered to be equal to the values related to Lower Cromer till [19] in Table 2. An important point that should be considered is that the preconsolidation pressure (${p}_{0}$) is changed by variation of each parameter during sensitivity analysis to maintain the value $\psi /{\psi }_{\text{R}}$ constant while the value of OCR is variable for all cases with $\psi /{\psi }_{\text{R}}\ne 1$. However, although the values of $p/{p}_{0}$ are different, all the cases related to $\psi /{\psi }_{\text{R}}=1$ are normally consolidated and their initial stress state is located on the YS at different stress ratios.

Table 1 Determination of the parameters that are original to the proposed model

Table 2 Parameters used in the model validations

A step-by-step procedure for calibrating new model-specific parameters is presented in the Appendix section. The new parameters’ calibration can be carried out with data from at least one conventional oedometer test combined with one undrained triaxial test. Given the natural variations between clay samples, a number of tests of the same type can be used to increase the reliability of estimated parameter values.

3.1 YS shape parameters

The flexibility of the adopted YS is provided via two shape parameters (${r}_{\text{y}}$ and ${n}_{\text{y}}$) and $N$ that specifies the bulkiness of the surface. These parameters provide the flexibility to adjust both the shape of the YS and the position of CSL relative to the NCL in the $e-\mathrm{ln}\, p$ plane. Consequently, their variations result in a wide range of model responses. Each parameter has a distinctive impact on the simulation capability of the model and will be discussed separately.

3.1.1 Parameter ${r}_{\text{y}}$

As shown in Fig. 2, this parameter mostly contributes in adjusting the intersection point of the YS and CSL in the stress space, and consequently the position of the CSL in the $e-\mathrm{ln}\, p$ plane, while the overall shape of the YS is almost unchanged. Therefore, it can be concluded that the variation of model responses with different values of ${r}_{\text{y}}$ is mainly due to the changes in the position of the CSL. For instance, regarding Eq. (28) and Fig. 9, higher values of ${r}_{\text{y}}$ move the CSL away from the ACL line, and simultaneously displace the intersection point to the descending section of the YS. As a result, the proposed model simulates a hardening behaviour that is followed by softening prior to the failure for the cases of normally and lightly overconsolidated samples (Fig. 11a and d). As can be seen, the softening part becomes more profound for higher values of ${r}_{\text{y}}$. For instance, in the case of the normally consolidated sample, by considering ${r}_{\text{y}}$= 6, the hardening part can be neglected and a fully softening behaviour is produced by the proposed model. However, there are no such consecutive hardening–softening trends for highly overconsolidated simulations (Fig. 11g), since during the loading procedure, unlike some cases of the normally consolidated and lightly overconsolidated stress states, the stress state or its related image point on the YS ($I{M}_{1}$) remains constantly on one side of the YS intersection point with the CSL (i.e., the highest point of the YS) and pure hardening responses are achieved.

Fig. 11

Sensitivity analysis for the yield surface shape parameter ${r}_{\text{y}}$ in different spaces for a–c normally consolidated, d–f lightly overconsolidated, and g–i highly overconsolidated samples

Another aspect that should be mentioned here is the fact that by increasing the ${r}_{\text{y}}$ value and considering constant $\psi /{\psi }_{\text{R}}$ and consolidation pressure (i.e., 200 kPa), the void ratio is decreased significantly (that might be unrealistic in some cases) that lead to the portrayal of more stiff samples, and coincidentally, the distance between the ACL and CSL is expanded which means, for example in the case of normally consolidated samples (Fig. 11b), more reduction in mean effective stress is required for the sample to reach the failure.

As a final point, parameter ${r}_{\text{y}}$ also affects the shear strength predicted by the proposed model for both normally and overconsolidated samples since any change in the parameter ${r}_{\text{y}}$ alters the width of both dry and wet sides of the YS. Hence, this parameter, besides other parameters that change the shape of the YS, improves the predicted shear strength which is one of the drawbacks of the elliptical YSs.

By considering that the location of the ACL is known (i.e., from one-dimensional consolidation test), ${r}_{\text{y}}$ can be calibrated using the undrained triaxial test results. This can be done through two different methods: (1) by calibrating the stress path, and/or (2) by calibrating the CSL in the $e-\mathrm{ln}\, p$ plane using data related to the void ratio and effective stress at the end of undrained triaxial tests when the critical state is achieved. Although both methods lead to reasonable initial assumptions, the authors suggest the first method while it implicitly specifies the location of the CSL.

According to the fundamentals of the CSSM, $\partial f/\partial p =0$ should be satisfied at the intersection point of the YS and CSL when an associated flow rule ($f=g$) is adopted, which of course is not the case in the proposed model. However, if an associated flow rule is assumed, the parameter ${r}_{\text{y}}$ can be defined as a function of other YS shape parameters

$${r}_{\text{y}}=\mathrm{exp}\left[\frac{{\left(M-\alpha \right)}^{{n}_{\text{y}}} }{{n}_{\text{y}}{\left(N-\alpha \right)}^{{n}_{\text{y}}-1}N}\right]$$

(32)

Although this assumption reduces the number of model parameters, it affects the model capabilities and restricts its flexibility in capturing the CSL and/or yield stress points.

3.1.2 Parameter ${n}_{\text{y}}$

Parameter ${n}_{\text{y}}$ mainly contributes in changing the overall shape of the YS in a large variety of shapes from teardrop- to bullet-shaped surfaces. However, regarding Eq. (28) and Fig. 9, it can also specify the location of the CSL in the $e-\mathrm{ln}\, p$ plane. In contrast to the case of parameter ${r}_{\text{y}}$ that mainly calibrates the position of the CSL, a mixing effect of both different YS shapes and spacing ratio ($R$) that comes with parameter ${n}_{\text{y}}$ is responsible for the diversity of model responses. Similar to parameter ${r}_{\text{y}}$, as ${n}_{\text{y}}$ increases the CSL moves away from the ACL and accordingly a softening behaviour emerges prior to failure. This softening behaviour is obvious for normally consolidated and lightly overconsolidated samples because of the position of the actual or image stress state relative to the YS apex (Fig. 12a and d).

Fig. 12

Sensitivity analysis for the yield surface shape parameter ${n}_{\text{y}}$ in different spaces for a–c normally consolidated, d–f lightly overconsolidated, and g–i highly overconsolidated conditions

An important point is that the variation of void ratio does not follow a specific trend as the parameter ${n}_{\text{y}}$ varies for different consolidation states. Since the void ratio increases by increasing ${n}_{\text{y}}$ for the normally consolidated state (Fig. 12b), it follows an inverse trend for highly overconsolidated samples (Fig. 12h). Moreover, in the case of lightly overconsolidated samples (Fig. 12e), the value of the void ratio increases as the parameter ${n}_{\text{y}}$ increases, and then it declines by further increasing the value of ${n}_{\text{y}}$. These diverse trends can be explained as follows. Since the void ratio remains constant during the undrained loadings, the void ratio related to the consolidation state can be determined as

$$e=\psi +{e}_{\text{CSL}}=\psi +{e}_{\Gamma }-\lambda \mathrm{ln}\, p$$

(33)

Considering Eqs. (27), (28) and (30), ${n}_{\text{y}}$ affects both ${\psi }_{\text{R}}$ (and similarly $\psi$) and ${e}_{\Gamma }$ but oppositely (i.e., increments of ${n}_{\text{y}}$ increases $\psi$ while decreases ${e}_{\Gamma }$). The participation of these parameters are dissimilar for different consolidation states, and the dominance of each of them controls the void ratio trend. For example, for highly overconsolidated samples, in general the void ratio is decreasing as the reduction of ${e}_{\Gamma }$ is dominating the trend, rather than the increase in $\psi$.

The value of ${n}_{\text{y}}$ can be determined using the undrained triaxial test data by trying to fit the stress path. Moreover, using experimental data related to the yield stress points in the stress space is another option to calibrate the value of ${n}_{\text{y}}$. However, due to uncertainties associated with determining the definitive yield stress points, the authors recommend fitting the triaxial undrained isotropic and anisotropic stress paths using the variation of ${n}_{\text{y}}$.

3.1.3 Parameter $N$

By changing the volume of the YS, parameter $N$ provides the model with a more desirable flexibility to capture the yield stress points and simulate the stress path and stress–strain responses. While this parameter also changes the position of the CSL in the $e-\mathrm{ln}\, p$ plane, different forms of simulations can be obtained. Figure 13 shows a range of model responses with regard to different values of $N/M$. As can be seen, this parameter has a comparable effect to the parameter ${r}_{\text{y}}$ on the model responses. Therefore, various modes of hardening–softening behaviour can be simulated by adopting suitable values of $N$. In addition, as is obvious in Fig. 13c, f and i, the value of $N$ has a great influence on the predicted shear strength. Moreover, the increase of $N$ leads to more stiff samples with lower void ratios, and more reduction in mean effective stress is required to reach the CSL (Fig. 13b, e and h).

Fig. 13

Sensitivity analysis for the yield surface shape parameter $N$ in different spaces for a–c normally consolidated, d–f lightly overconsolidated, and g–i highly overconsolidated conditions

After calibrating the parameters ${r}_{\text{y}}$ and ${n}_{\text{y}}$ to capture the overall trend of experimental test data, the model response can be fine-tuned by adjusting the value of $N$.

3.2 Plastic potential shape parameters

As discussed before, the proposed model is enhanced with a PPS that uses two shape parameters (${n}_{p}$ and ${m}_{\text{p}}$) to bring the desired flexibility into the model to capture both the stress–strain behaviour and the one-dimensional compression loading response more realistically.

The parameters associated with the PPS mostly contribute to the adjustment of the stress–strain response while they have an almost insignificant effect on the stress path when compared to the YS shape parameters (i.e., compare Figs. 14 and 15 with 11 and 12), especially in the case of normally consolidated and highly overconsolidated samples. As can be seen, by increasing ${n}_{p}$, the proposed model tends to predict less brittle response by decreasing the peak shear strength. However, ${m}_{\text{p}}$ affects the model differently and amplifies the peak shear strength (i.e., compare Figs. 14 and 15). Furthermore, these two parameters adjust the rate of deviatoric stress reduction after reaching the peak shear strength. For instance, high values of ${n}_{p}$ cause a more rapid reduction in the deviatoric stress (Fig. 14b and d), while an opposite trend can be observed by increasing ${m}_{\text{p}}$ (Fig. 15b and d).

Fig. 14

Sensitivity analysis for the plastic potential shape parameter ${n}_{p}$ in different spaces for a, b normally consolidated, c, d lightly overconsolidated, and e, f highly overconsolidated conditions

Fig. 15

Sensitivity analysis for the plastic potential shape parameter ${m}_{\text{p}}$ in different spaces for a, b normally consolidated, c, d lightly overconsolidated, and e, f highly overconsolidated conditions

In addition to adjusting the stress–strain responses, the PPS shape parameters enable the proposed model to simulate the ${K}_{0}$ values more accurately. Figure 16 demonstrates how different values of ${n}_{p}$ and ${m}_{\text{p}}$ regulate the model response under the application of one-dimensional loading. As can be seen, variation of the PPS shape parameters changes the stress ratio that the stress path tends to reach (Fig. 16a and c) and its corresponding equilibrium state of anisotropy (Fig. 17b and d). Along with this, Fig. 17 shows the variation of ${K}_{0}$ simulated by the proposed model based on Eq. (7) for different friction angles compared with Jacky’s empirical formula (${K}_{0}=1-\mathrm{sin}\varphi$). As can be seen, the capability of the adopted PPS and also the flexibility associated with the RH rule that pushes the YS to a prescribed equilibrium state of anisotropy ensure that the ${K}_{0}$ value can be adequately predicted [15].

Fig. 16

Variation of model responses under one-dimensional compression condition for different values of a, b ${n}_{p}$ and c, d ${m}_{\text{p}}$

Fig. 17

Variations of ${\eta }_{{K}_{0}}$ and ${K}_{0}$ predicted by the proposed model for different values of friction angle $\phi$ by considering different values of a, b ${n}_{p}$ c, d ${m}_{\text{p}}$

From the calibration point of view, these parameters can be adjusted using the experimental data from one-dimensional compression tests. Alternatively, one can use the stress–strain response from the undrained triaxial tests to calibrate the PPS parameters.

3.3 Rotational hardening rule parameters

The proposed RH rule is versatile enough to consider various aspects of clay anisotropy, despite the constant values of 5, 2, 1, and 1 that can be assigned to its parameters $a$, $b$, $c$, and ${\chi }_{v}$, respectively [15]. The two remaining parameters, ${\chi }_{d}$ and $\mu$, conserve the overall capabilities that should be provided by a RH rule. The role of each of these two parameters is discussed in the following.

3.3.1 Parameter $\mu$

Parameter $\mu$ controls the absolute pace of YS rotation toward its equilibrium state. Although this parameter controls how the YS rotates during shearings, as illustrated in Fig. 18, it has an insignificant effect on the overall prediction of the triaxial behaviour. As a recommendation, values between 3 ${p}_{\text{atm}}$ to $7{p}_{\text{atm}}$ can be considered for this parameter.

Fig. 18

Sensitivity analysis for the RH parameter $\mu$ in different spaces for a, b normally consolidated, c, d lightly overconsolidated, and e, f highly overconsolidated samples

3.3.2 Parameter ${\chi }_{d}$

This parameter defines the equilibrium state of anisotropy which the YS attempts to incline along with when the stress state reaches the CSL. Since ${\chi }_{d}$ changes the position of both ACL and CSL in the $e-\mathrm{ln}\, p$ plane, it has a similar effect as the YS shape parameters on the simulations. Figure 19 shows the modelling of undrained triaxial responses using various ${\chi }_{d}$ values. As shown, the variations of the CSL in the $e-\mathrm{ln}\, p$ plane lead to different pre-failure stress states, while at the low-stress levels all stress paths follow a similar trend. As discussed earlier, higher values of ${\chi }_{d}$ lead to wider distances between the ACL and CSL, and a softening behaviour appears in the model response which is associated with a reduction in the mean effective stress. Since ${\chi }_{d}$ can change the initial inclination of the YS, part of the differences in simulation results can be associated with changes of ${\alpha }_{0}$.

Fig. 19

Sensitivity analysis for the rotational hardening parameter ${\chi }_{d}$ in different spaces for a–c normally consolidated, d–f lightly overconsolidated, and g–i highly overconsolidated conditions

Figure 20 demonstrates how the variation of ${\chi }_{d}$ can reproduce different $\eta$ under the application of one-dimensional loading. Since ${\chi }_{d}$ changes the equilibrium state of anisotropy, the corresponding stress ratio is readily achievable by calibration of this parameter. Additionally, as demonstrated in Fig. 21, the variation of this parameter also changes the predicted ${K}_{0}$ value through the orientation of the YS and PPS under the constant ${\eta }_{{K}_{0}}$ stress ratio (i.e., ${\alpha }_{{K}_{0}}$ in Eq. 7).

Fig. 20

Variation of model responses under one-dimensional compression condition for different values of ${\chi }_{d}$

Fig. 21

Variations of ${\eta }_{{K}_{0}}$ and ${K}_{0}$ predicted by the proposed model for different values of friction angle $\phi$ by considering different ${\chi }_{d}$

With the above explanations, data from either one-dimensional loading or undrained triaxial tests can be used to calibrate ${\chi }_{d}$. It is noteworthy that the equilibrium inclination at the CSL (${\chi }_{d}M$) must be less than $N$ to have a valid YS ($\left({N}^{2}-{\chi }_{d}^{2}{M}^{2}\right)>0$).

3.4 DISP parameters

The novel DISP method enables the model to simulate the behaviour of overconsolidated samples, as well as cyclic responses. Using only one additional parameter ($h$), the proposed method enhances the model to capture the nonlinear behaviour inside the YS, and also to simulate the cyclic behaviour. This parameter clearly affects the model predictions in both lightly and highly overconsolidated conditions (Fig. 22). Regarding Eqs. (18) and (19), when parameter $h$ moves toward large numbers, the plastic modulus of the stress state inside the YS also increases, and ultimately, when $h$ reaches infinity, the proposed model turns into a classical elastoplastic model.

Fig. 22

Sensitivity analysis for the DISP parameter $h$ in different spaces for a, b lightly overconsolidated and c, d highly overconsolidated conditions

4 Model validation

The simulative performance of the proposed model is evaluated using available experimental data of two different anisotropically consolidated clays subjected to conventional drained and undrained triaxial loading. In addition, two cyclic hollow cylinder tests on normally consolidated and overconsolidated clay samples are modelled by the proposed constitutive model to demonstrated the efficiency of the DISP method. Prior to simulations, and following the procedures described in the previous section, the model parameters are first calibrated. Table 2 summarises the model parameter values for each soil type.

4.1 Simulation of monotonic behaviour

4.1.1 Lower Cromer till

LCT can be categorised as a low plasticity sandy clay with a liquid limit of 25%, a plastic limit of 13% and a plasticity index of 12%, which has been considered as clay in modelling literature (e.g., [5, 7]). The mineralogical analysis shows that the main constituent of LCT is quartz (more than 50%) and about 17% is composed of clay minerals (i.e., calcite, illite, smectite, kaolinite, chlorite, etc.) with the clay activity of 0.71. Moreover, the specific gravity of the coarser particles (${G}_{s}$) was determined to be equal to 2.65 [65]. Gens [19] conducted drained and undrained triaxial tests on isotropically and anisotropically reconstituted consolidated samples of LCT, with different OCRs, which were taken from uniform specimens.

Figure 23 compares the simulation results against the experimental data related to the anisotropically consolidated samples. Both predicted stress path and stress–strain response are in outstanding agreement with the experimental data. The remarkable simulation of the normally consolidated sample behaviour confirms the efficiency of the model features. Also, in case of the overconsolidated samples, besides the effectiveness of different model components (i.e., YS, PPS and RH), the proposed DISP method clearly leads the model to reproduce noticeably good predictions.

Fig. 23

Comparison of data and simulations for undrained triaxial tests on anisotropically consolidated samples of Lower Cromer till (data from Gens [19]); a stress path, and b stress–strain responses

An interesting observation in Fig. 23 simulations is the sickle-shaped stress path related to the sample with OCR = 2, where given the marked agreement between the model response and the test data, one can consider it as a distinctive improvement of the proposed model over the existing anisotropic clay models that generally fail to reproduce such a stress path (e.g., refer to Fig. 18 of [5] and Figs. 5, 6 and 7 of [9]). These types of simulations can be justified by considering the variation of the CSL with regard to anisotropy during straining, which was earlier mentioned as the moving CSL. At the early stages of the loading, the stress path moves toward the position of the CSL related to the initial state of anisotropy (i.e., ${\alpha }_{0}$). By the progression of loading, the anisotropy of the sample grows and the CSL tends to move away from its initial position. The turning point of the sickle-shaped path is where the CSL crosses the stress state (consider it in the $e-\mathrm{ln}\, p$ plan), and afterwards, the mean effective stress starts to decrease according to its tendency to approach and touch the CSL. Finally, the stress state and the CSL intersect with each other when the anisotropy reaches the equilibrium state. This aspect of soil behaviour highly depends on the initial material and stress states, where the initial position of the CSL (regarding ${\alpha }_{0}$) is on the right side of the stress state, but the final position (regarding ${\alpha }_{\text{e}}$) lies on the left side. For the case of highly overconsolidated samples, the initial and final positions of the CSL are located on the right side of the stress state, and this is why the stress path only tends to go to the right side.

Furthermore, Figs. 24 and 25 show simulation of drained triaxial tests on isotropically and anisotropically consolidated samples of LCT. As can be seen, the model shows good agreement with the test data, specifically in terms of the volume change and peak shear strength. Figures 23b and 24b show that by increasing the OCR, the tendency to contractive behaviour is reduced, and simultaneously more dilative responses are observed (e.g., fully dilative behaviour of isotropically consolidated sample with OCR = 10).

Fig. 24

Comparison of data and simulations for drained triaxial tests on isotropically consolidated samples of Lower Cromer till (data from Gens [19]); a stress–strain, and b volumetric strain–axial strain response

Fig. 25

Comparison of data and simulations for drained triaxial tests on anisotropically consolidated samples of Lower Cromer till (data from Gens [19]); a stress–strain, and b volumetric strain–axial strain response

4.1.2 Boston blue clay

Boston blue clay (BBC) is a low plasticity type marine clay composed of illite and quartz. This soil has a saturated unit weight of 18.26–19.22 kN/m³ and its dry unit weight is between 12.81 and 13.93 kN/m³. The liquid limit, plasticity index and liquidity index of BBC are 41%, 21% and 0.8%, respectively. Ladd and Varallyay [33] carried out a series of undrained triaxial tests on samples of BBC with different consolidation states. The initial void ratio of the specimens was varied between ${e}_{0}$ = 0.84–0.89 and with the preconsolidation pressure ranging from 273 to 785 kPa. The reported experimental data over this clay have been widely used by researchers in the soil constitutive modelling field (e.g., [5, 7, 25, 63]). Here, the values of conventional critical state parameters for BBC are taken from Ling et al. [36].

Figure 26 shows the model predictions of stress paths and stress–strain responses for the anisotropically consolidated BBC samples. Similar to the case of LCT, the figure confirms the capability of the proposed model in reproducing very good simulations of the observed behaviour in both normally consolidated and overconsolidated conditions. For the case of OCR = 2, the concept of moving CSL governs the stress path, through which the mean effective stress initially increases followed by a reduction prior to touching the CSL.

Fig. 26

Comparison of data and simulations for undrained triaxial tests on anisotropically consolidated samples of Boston blue clay (data from Ladd and Varallyay [33]); a stress path, and b stress–strain responses

In summary, regarding the simulation of monotonic responses of soil, the following remarks can be highlighted:

• The remarkable simulations of the normally consolidated samples confirm that the set of YS and PPS combined with the versatile RH rule provides the proposed model with the capability to capture the behaviour of anisotropically consolidated samples, including realistic prediction of the peak shear resistance and pre-failure softening behaviour.

• The proposed DISP method, using one parameter ($h$), enables the model to predict the behaviour of highly overconsolidated samples with very good accuracy. The position of the projection centre, which is determined based on the stress and anisotropy states, as well as the shape of the YS (Eq. 21), plays a vital role in the enhanced modelling of this condition.

• The concept of the moving CSL has a considerable effect on the simulation of lightly overconsolidated samples where the CSL may cross the stress state (in the $e-\mathrm{ln}\, p$ plane). In both conducted simulations of lightly overconsolidated samples, the moving CSL enables the model to capture the soil behaviour more accurately compared to the existing models in the literature that try to simulate this condition (e.g., [5, 7, 9]).

4.2 Simulation of cyclic behaviour

Using hollow cylinder apparatus, Soralump and Prasomsri [59] carried out undrained cyclic multistage strain-controlled tests on compacted clays obtained from the core zone of an earth dam in Thailand. The specimens were prepared at the optimum moisture content and maximum dry density and consolidated isotropically at different OCRs ranging from 1 to 4. Each test included several cyclic strain-controlled stages with shear strain levels.

In order to demonstrate the capability of the proposed DISP method in predicting the cyclic behaviour, the responses of two samples of normally consolidated and overconsolidated (OCR = 4) conditions are simulated in this section. Both samples were prepared at 5 layers with the initial void ratio of ${e}_{0}=0.546$ and then consolidated at the maximum confining pressure of 400 kPa. The normally consolidated sample was subjected to multi-stage cyclic strains ranging between 0.004 and 0.5%. The confining pressure of the other sample was reduced to 100 kPa to reproduce OCR = 4, and then, it was subjected to multi-stage cyclic strains ranging between 0.004 and 1.5%. Here, the values of conventional soil parameters are taken from Elia and Rouainia [17].

Figures 27 and 28 compare the experimental data with model responses (in terms of stress amplitudes and stress–strain loops) for samples subjected to the multi-stage cyclic strains in normally consolidated and overconsolidated conditions, respectively. As can be seen, using a single additional parameter, the proposed DISP method provides the model with the capability to predict the cyclic behaviour with reasonable accuracy. Figures 27d and 28d show that the model simulations, in terms of normalised cyclic shear stress amplitude and its degradation with the number of cycles, are in good agreement with the test data (in Figs. 27c and 28c).

Fig. 27

Comparison of cyclic shear strength and stress–strain hysteretic loops for normally consolidated compacted clay (data from Soralump and Prasomsri [59]); a, c, and e experimental data, and b, d, and f model simulations

Fig. 28

Comparison of cyclic shear strength and stress–strain hysteretic loops for overconsolidated (OCR = 4) compacted clay (data from Soralump and Prasomsri [59]); a, c, and e experimental data, and b, d, and f model simulations

In addition, Figs. 27f and 28f compare the simulated stress–strain loops with the observed responses (Figs. 27e and 28e). As can be seen, the proposed model is able to capture the overall trend of soil response well; however, there is an apparent discrepancy between the predictions and experimental data. Elia and Rouainia [17] reported the same observation and related this difference to the nature of the specimen that was used in the experiment, where it must have contained a high percentage of non-fine materials. The presence of these coarse grains leads to cyclic mobility response, similar to what might be expected when a granular material is subjected to cyclic loading with large strain levels. Moreover, based on the shape and size of the predicted stress–strain loops in Figs. 27f and 28f, it is shown that the proposed model is able to replicate the reduction of shear stiffness, and also consider the increment of damping as a function of shear strain levels.

As an important point, by considering a constant rate for monotonic and cyclic loading, the proposed model is supposed to simulate both responses with the same set of model parameter values. However, since the test data for a specific clay under both cyclic and monotonic loading conditions are not available in the literature, the accuracy of model response for the simulations with a unique set of parameters is yet to be proven.

5 Conclusion

A non-associative adaptive anisotropic constitutive model, named AA2-DISP, has been developed within the critical state soil mechanics framework to simulate the behaviour of natural and reconstituted clays under the application of cyclic and monotonic loadings. A flexible rotated YS is incorporated into the model to reproduce a wide range of shapes that are consistent with the experimental yield stress points related to the clays with preferred fabric orientations. The formulation of the proposed model also takes advantage of an inclined flexible PPS which improves the accuracy of both stress–strain response and predicted ${K}_{0}$ values.

Moreover, the proposed model uses an enhanced RH rule to determine the rotations of both YS and PPS. According to experimental evidence, a hyperbolic function is considered in the proposed RH rule to define the governing plastic strain increment that varies between the volumetric and deviatoric plastic strain increments for low and high values of stress ratio, respectively. In addition, the adopted RH rule guarantees a unique equilibrium state of anisotropy during plastic shearing that both YS and PPS tend to reach, which is crucial for the stress state to reach the CSL in the $e-lnp$ plane. Regarding this, the concept of moving CSL has been introduced in the model framework, which refers to the dependency of the CSL on the soil fabric anisotropy that has been ignored by the available constitutive models for clays. By formulating and regulating the movement of the CSL in the $e-lnp$ plane, the simulation capability of the model is markedly enhanced, especially for lightly overconsolidated samples.

Finally, an elegant DISP method, which is analogous to the well-known bounding surface theory, is developed and incorporated into the model formulation to capture the small strain nonlinearity and also to reproduce cyclic responses. This method facilitates the definition of a moving PC to simulate the hysteretic behaviour during loading–unloading–reloading cycles, as well as the monotonic behaviour of highly overconsolidated samples more accurately with the least amount of additional parameters.

In total, the AA2-DISP model has eight additional parameters compared to the idealised modified Cam-clay model. An extensive sensitivity analysis has been carried out, that can be considered as a clear guide for model parameter determination. It also allows the readers to understand the model formulation, its predictive capabilities and its advantages over a number of well-established anisotropic models for clays. As an original approach, the parametric study of the proposed model was also elaborated in the $e-lnp$ plane, in addition to the conventional stress path and stress–strain spaces which are more customary in the soil constitutive modelling field. This helps the readers to comprehend the governing mechanisms of the constitutive responses at the pre-failure and failure stages. The model was also employed for the prediction of the responses of three different clayey soils subjected to monotonic and cyclic loadings. The comparison of model predictions against experimental data illustrated that the model is remarkably successful in capturing clay responses under different conditions.

Appendix

The calibration of new model-specific parameters can be performed by using the step-by-step procedure presented in the following flowchart (Fig. 

Fig. 29

Model parameters calibration procedure

29).