Adaptive anisotropic constitutive modeling of natural clays

In this paper, an adaptive anisotropic constitutive model, namely AA1‐CLAY, is developed for clays based on the critical state framework. The model has a non‐associated flow rule, and it can also reduce to basic critical state constitutive models. A versatile yield surface (YS) is implemented in the model which can generate a broad range of shapes to improve the possibility of capturing experimental yield points of different clay types with high accuracy. In addition to the isotropic hardening rule, an innovative rotational hardening (RH) rule is also incorporated into the model to control the YS rotation rate. The proposed RH rule uses a transitional function to govern the effect of plastic strains at different constant stress ratios more realistically. Furthermore, the equilibrium state of anisotropy, which the YS tends to reach during plastic straining, is defined uniquely based on the experimental findings. The detailed model formulation is presented, and important features of the model are elaborately discussed. The capabilities of the model are also demonstrated by comparing the simulation results of several element tests on a number of different clays against the available experimental data.


INTRODUCTION
Many developments in advanced constitutive modeling of clays are associated with the critical state soil mechanics (CSSM) theory. As the earliest elastoplastic constitutive model based on the CSSM concept, the original Cam-clay (OCC) model has improved the prediction of isotropically consolidated soft soils. 1,2 In OCC, a bullet-shaped surface was adopted for both yield surface (YS) and plastic potential surface (PPS), and a simple isotropic hardening rule was employed to represent the evolutions associated with plastic strains. This model made a significant amount of progress in predicting various experimentally observed responses of isotropically consolidated soft clays. However, as OCC cannot distinguish between isotropic and one-dimensional normal loading, and to solve the computational difficulties associated with the presence of an acute tip in its bullet-shaped surface, an ellipsoid was later adopted for the YS and PPS which led to the development of modified Cam-clay (MCC) model. 3 Ignoring the issues related to the application of single YS in the basic CSSM-based models (such as inability in capturing the nonlinear responses in the small strain domain or inability in simulating the cyclic behavior), Cam-clay family of models still suffer from two other major deficiencies. Firstly, these models fail to simulate the behavior of natural and anisotropically reconstituted soils. And secondly, due to their bulky YS shape in the dry side of the critical state line (CSL), their performance in capturing the yield points of the overconsolidated clays is unsatisfactory. In other words, they can overestimate the peak strength of overconsolidated clays.
Many experimental data have confirmed that yield stress envelope of highly overconsolidated clays in triaxial stress space forms an approximately straight line on the dry side of the CSL rather than an MCC type ellipse. 4 This yield locus, in the general stress space, is termed as the Hvorslev surface in honor of the pioneering work of Hvorslev. 5 Several constitutive models have incorporated the Hvorslev surface as a straight or slightly curved line into their YS to predict the behavior of highly overconsolidated clays more precisely. [6][7][8][9][10][11] Since piecewise YSs cause computational difficulties due to the discontinuous joining points which violate the uniqueness of YS gradient vector condition, adopting continuous YS shapes which can approximate the Hvorslev surface is considered as the most coherent approach. This being the case, different CSSM model families with improved YSs have been developed in recent years with a wide range of theoretical backgrounds, for example, by experimental data fitting or employing thermodynamics theories. [12][13][14][15][16][17] These models are able to produce a wider range of YS shapes; hence, not only they can capture yield points more appropriately, but also they open up opportunities to develop unified constitutive models for both clay and sand. 18,19 The arrangement and orientation of soil particles or fabric, and their evolution during straining affect the behavior of soils through both elastic and plastic deformation domains. Many studies have confirmed that stress-strain relationship and shear strength of soils depend primarily on their fabric. By definition, the anisotropic configuration of solid particles leads to directional dependence of mechanical properties which may result in planes with lower shear strengths and even different stress-strain responses. Due to deposition history or the initial stress conditions, natural clays are considered as intrinsically anisotropic materials. This type of anisotropy, which is known as inherent anisotropy, might be adjusted during stress variation and deformations. The evolved fabric composition, which shows the influence of current stress conditions, is called stress-induced anisotropy.
It is a well-established fact that the yield points obtained from tests on undisturbed specimens of natural clay constitute YSs that are inclined in the representative stress spaces. This inclination of the YSs is largely agreed to be due to the inherent fabric anisotropy that exists in the structure of the clay. [20][21][22][23][24] A popular approach in considering the effects of inherent anisotropy within an elastoplastic modeling framework is development of elastoplastic soil constitutive models that involve inclined yield curves. [25][26][27][28][29] This methodology for representing the inherent anisotropy has its roots in the works of Sekiguchi and Ohta, 30 and Hashiguchi,31 and were subsequently followed by other researchers. 17,21,[32][33][34][35][36][37][38][39][40] In these works, the inclined YS is either fixed 30,41 or variable by adopting a rotational hardening (RH) rule to simulate the development or erasure of anisotropy during plastic straining. 17,[36][37][38][39][42][43][44][45][46] How plastic strain components participate in rotating the YS is one of the key aspects in proposing different RH rules. While various anisotropic constitutive models have been developed by correlating the evolution of the YS inclination with the volumetric plastic strain increment, 21,[32][33][34][35]38 Wheeler et al. 36 and Yang et al. 39 introduced constitutive models which adopt both volumetric and deviatoric plastic strain increments in their RH rules. Defining the RH rule as a function of only deviatoric plastic strain increment has also been examined in the work of Chen and Yang 17 .
In the present study, an adaptive anisotropic constitutive model based on the CSSM framework is presented for simulating the behavior of clays. By introducing two new shape parameters, a flexible YS is defined and incorporated into the model to produce various YS shapes. The proposed YS equation shows very promising capabilities in predicting the yield points both in dry and wet sides of the CSL. A non-associated flow rule formulation is also adopted by using an inclined PPS ellipse which satisfies the condition of zero volumetric plastic strain increment at the critical state. While the proposed model does not use any RH rule for PPS and considers identical inclinations for both YS and PPS during the evolution of plastic strains, they are generally not similar in shape and size. Moreover, the model uses an innovative RH rule, with a new expression, for both YS and PPS rotations, which also ensures the equilibrium state of anisotropy (which can also be termed as steady state) during plastic shearing which YS and PPS tend to reach. Based on observations from experimental results, a transitional function is implemented into the proposed RH rule to employ the governing plastic strains (volumetric and deviatoric) and equilibrium state more realistically at different constant stress ratio loading. Additionally, by choosing suitable model constants the development and erasure of fabric anisotropy (until reaching isotropic fabric state) at critical state can be modeled with the proposed RH rule. It is noteworthy to mention that the proposed model can be reduced to simpler well-documented constitutive models, as its YS shape and RH rule can be transformed into those of the previous models by choosing appropriate parameter values.
The layout of the paper is as follows. In Section 2, the general formulation of the model is presented in the triaxial stress space. The distinct features of the proposed model are discussed in Section 3. To determine the contribution of model parameters on the predictions, the calibration procedure is performed in addition to a detailed sensitivity analysis in Section 4 using the experimental results of Taipei silty clay (TSC). In the following section, the model performance is then verified against the test data of several clays including Kaolin clay (KC), Lower Cromer till (LCT), Boston blue clay (BBC) and Santa Clara clay (SCC). Finally, concise conclusions are provided in the last section.

MODEL FORMULATION
AA1-CLAY is an elastoplastic critical state-based constitutive model with an inclined YS to represent the anisotropic composition of clays in their natural state. This model can be considered as the first version and basis of a constitutive model family for natural and reconstitute clays. The proposed constitutive model benefits from a flexible 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. Furthermore, in this model, a non-associated flow rule is adopted by introducing a rotated and sheared MCC type ellipse as the PPS. An innovative RH rule is also introduced, which predicts the evolution of anisotropy more realistically by using a transitional function between plastic volumetric and deviatoric strain increments. Moreover, an equilibrium state of anisotropy is defined through the RH rule to interrupt the YS from excessive rotation when the stress state reaches the CSL. Finally, a new expression for the initial inclination of the YS is presented based on the employed RH rule. In the following, specific features and components of AA1-CLAY are discussed. It must be noted that all stress-related parameters used in this paper represent effective stresses, and the prime sign is therefore neglected. Moreover, the model formulation is presented in the triaxial or − space for the sake of simplicity, and the general stress formulation is provided in Appendix A.

Yield surface
Using expressions with more degrees of freedom to produce a wide range of different shapes (i.e., from a simple ellipse form to a complex inclined teardrop or bullet-shaped surfaces) can provide a higher level of accuracy in capturing the experimental yield stress points, and also enhance the model capabilities for more realistic stress path and stress-strain predictions. For instance, following experimental observations from undrained triaxial tests confirming that failure envelops of the heavily overconsolidated clays lie on the Hvorslev surface 4 , teardrop-shaped yield surfaces were proposed to substitute with the Hvorslev surface in the dry side. 12,18,47 This type of YS expressions also provided the opportunity to define unified models for clays and sands 18,19 ; since constitutive models of sands mostly use a thin open or close wedge as their YS because shearing and variation of stress ratio are known as the primary cause of plastic strain in non-loose granular soils. 48 Additionally, the interpretation of triaxial test results on some isotropically consolidated clays shows that models with a bullet-shaped YS reproduce more accurate simulations. 49 Compared to established simple anisotropic clay models SANICLAY 37 and S-CLAY1, 36 the YS in the proposed model has two additional parameters, and , which enable the model to reproduce a wide range of shapes from a simple ellipse to tear or bullet-shaped YSs and also leaf-shaped YSs with two sharp tips. The YS is expressed in the triaxial stress space as where is a stress-ratio type variable that defines the inclination of YS to represent the plastic anisotropy, 0 is preconsolidation pressure and represents the size of the YS, is the shape factor of the YS, and and are two parameters controlling the shape of the YS. Figure 1 illustrates the YS shape variations with different values of and . As this figure shows, parameter changes the overall shape of the YS while parameter mainly contributes to the adjustment of the YS curvature, especially in the dry side. The possibility of calibrating these two parameters results in enhanced pliability of the model in capturing yield points and subsequently improved undrained soil behavior simulations in both normally consolidated and overconsolidated states. By setting = 1 and = 0 the YS of the model reduces to those of the former constitutive models, as is illustrated in Figure 2. One can substitute shape factor, , with critical state stress ratio, , to obtain the S-CLAY1 model ( = ) or leave it as a constant model shape parameter to replicate SANICLAY constitutive model ( ≠ ).

Plastic potential and flow rule
In addition to determining the relative magnitude of plastic strain increments, the PPS also specifies the ultimate strength.
Although some studies suggest that using identical YS and PPS is a feasible assumption for constitutive models with an oriented YS 19,36,50 , it was found that employing associated flow rule leads to unrealistic simulations. 12,18,[37][38][39]51 Yu 18 showed that models using associated flow rule fail to predict the undrained response of normally consolidated undisturbed clays with softening and pre-failure peak strengths. Taking this point into account, several recently developed models employ the non-associated flow rule in their formulation. 16,17,[37][38][39]52 AA1-CLAY uses the same PPS shape of those in S-CLAY1 and SANICALY models, which is an oriented and sheared MCC type surface and it is defined in triaxial stress space with the following expression where and represent the inclination and the size of the PPS, respectively. It is noteworthy to mention that in the proposed model, the non-associated flow rule is in order, unless when and shape parameters are set to 0 and 1, respectively, and = . The proposed model does not employ any independent RH rule for PPS, and it assumes an identical inclination for YS and PPS ( = ). To set up the non-associated flow rule and determine the PPS gradient at the current stress state, it is necessary to find the size of the PPS by solving Equation (2). As the PPS is defined, now the plastic volumetric and deviatoric strain increments can be determined readily using the following expressions where = ∕ is stress ratio, is loading index, and ⟨ ⟩ represent Macaulay brackets.

Hardening rule
The new AA1-Clay model adopts isotropic and RH rules to represent any changes in size and inclination of the YS due to plastic volumetric and shear strains.

Isotropic hardening
For the expansion and contraction of the YS, the model adopts a volumetric hardening rule with the same general formulation of that of the MCC model where the changes in the size of the YS are related to the plastic volumetric strain increments, , as in which is the specific volume, and are the slopes of the normal compression and swelling lines in the (ln − ) space, respectively.

Rotational hardening
As stated above, the proposed model belongs to the class of anisotropic models, where anisotropy is introduced through the initial inclination of YS, and the evolution of YS during the development of plastic strains is controlled by a RH rule. In the literature, different combinations of the plastic strain increment components are considered to contribute to the soil anisotropy. While Banerjee and Yousif 32 , Dafalias 21 , Whittle and Kavvadas 34 , Ling et al. 35 , Dafalias et al. 37 , Jiang et al. 38 presented RH rules as a function of plastic volumetric strain increment, Chen and Yang 17 introduced a RH rule by incorporating only plastic deviatoric strain increment. Involving both components of plastic strain increment in the evolution of anisotropy is another option that can be found in the works of Pestana and Whittle 19 , Wheeler et al. 36 , Nieto Leal et al. 53 by assuming separate contribution of plastic strain components in the evolution of anisotropy. Moreover, Yang et al. 39 merged the actions of plastic strain components through the root mean square of plastic volumetric and deviatoric strain increments. One of the most crucial aspects of a RH rule is the definition of an equilibrium state that attracts the YS during nonhydrostatic loading and plays the role of a limiting value which controls the excessive rotation of the YS. 44 Under the application of a constant-stress-ratio loading, a specific configuration of soil fabric is formed in its structure that induces a preferred inter-particular state. One realization of this fabric configuration, that also controls the soil's shear strength and stress-strain behavior, is a specific inclination of the YS in the stress space that is termed as equilibrium state of anisotropy. Dafalias 21 considered a fraction of as the bounding value by means of a model constant (Equation 12 in the reference). Although this formulation provides adequate results along with its simplicity, it cannot restrict the rotation of the YS when > . As a remedy, Dafalias and Taiebat 51 redefined parameter as a function of through an exponential function to curb the excessive rotation the YS, and showed that the modified version provides more conformity with experimental results. To address the hypothesis of isotropic fabric orientation during critical state failure, Dafalias, and Taiebat 44 improved their last premise and proposed an expression for the equilibrium state for the range of | | ≥ with a declining trend for bounding value to zero as the quantitative representation of the isotropic state. Wheeler et al. 36 defined equilibrium state of fabric anisotropy implicitly through their RH rule which can be obtained by solvinġ= 0 (explained in the sequel). The equilibrium value varies between , according to the values of the model constant ( Figure 5 in the reference). It should be noted that similar to Dafalias 21 , Wheeler et al. 36 did not anticipate any measures to limit the equilibrium state of anisotropy when > which leads to the excessive rotation of the YS. 51 Yang et al. 39 equipped their proposed RH rules for both YS and PPS with two distinct equilibrium states. The integration of the proposed RH by Chen and Yang 17 led to an exponential function of plastic deviatoric strain increment. Although no explicit bounding value was defined in their model, it was shown in their work that the accumulative anisotropy reaches to zero as the plastic deviatoric strain increases infinitely on reaching the critical state.
In this work, a novel RH rule is proposed to describe the rotation of both YS and PPS. The RH formulation represents the development and erasure of anisotropy due to both plastic volumetric and shear strain increments. In the triaxial stress space, this hardening rule has the form in which controls the absolute rate of anisotropy evolution toward its equilibrium state ( ), and represent the bounding values which the YS tends to incline along with due to pure plastic volumetric and deviatoric strains under constant stress ratio loading condition when = , and , , and are model constants that adjust the equilibrium state at different stress ratios. Moreover, ∕ 0 is implemented to slow down the rate of YS rotation with increasing overconsolidation ratio ( ). As it is shown in Equation (5a), the proposed RH rule employs both plastic volumetric and deviatoric strain increments. It controls their participation through a transitional function ( ) which considers a gradual change in the contribution of governing plastic strain increments from volumetric to deviatoric by increasing stress ratio ( ). This function is defined in such a way that it satisfies the dominance of plastic volumetric strain at low values of stress ratio, while it signifies the role of plastic deviatoric strain simultaneously with the increment of stress ratio.
Moreover, a new expression for the equilibrium state of anisotropy is defined via Equation (5b). The first term in the square bracket characterizes the bounding value when ≤ , while the second exponential term is considered to restrain the rotation of YS when ≥ . By specifying and as two extreme equilibrium states and considering the significance of dominant plastic strain increments by hyperbolic function , is determined as a function of ∕ . Figure 3 shows the variation of with regard to parameters , , and . Parameters and prescribe the equilibrium state of anisotropy between and with regard to the governing plastic strain component when ≤ ( Figures 3A and 3B). On the other hand, parameter controls the variation of equilibrium state for loading beyond the critical state (see Figure 3C). While 0 < < 0.9 may lead to excessive rotation of YS as a consequence of the development of anisotropy, > 0.9 results in the erasure of anisotropy and decreases the inclination of the YS.
It should be noted that in some anisotropic models such as S-CLAY1 36 and MIT-S1 19 , Macaulay brackets are applied to the plastic volumetric strain increments. This is because in these models, yielding on the dry side of the critical state ( > and < 0) can result in excessive rotation of the YS and divergence from the equilibrium state at an accelerating pace. However, as discussed above, for the proposed model the second exponential term in Equation (5b) prohibits any deviation from the equilibrium state. Therefore, in the new model, contributions from the actual values of plastic volumetric strain increments (whether positive or negative) can be considered for fabric evolution and application of Macaulay brackets is not required.
As is demonstrated in Figure 4, the previously proposed equilibrium states of anisotropy can be reproduced by choosing relevant values of RH parameters, as are shown in Table 1. Also, by assuming = , the variation of equilibrium anisotropy state transforms to a straight line as in Nieto Leal et al. 53 Moreover, Figure 5 shows the comparison of actual and simulated values of the equilibrium state of anisotropy for Otaniemi clay. 36 With calibrated values of = 0.4 and

Discussion of the parameter
Controlling the pace of anisotropy evolution via a constant parameter is an integral part of almost every anisotropic constitutive model (e.g., in Wheeler et al 36 , in Dafalias et al. 37 and in Chen and Yang 17 ). This parameter directly affects the  model simulations by controlling the rate of YS inclination. Moreover, it has an implicit contribution to alter the overall response of the model, which has not been clearly elaborated in the literature. Besides controlling the absolute pace of anisotropy evolution, parameter also specifies the participation rate of isotropic and RH rules in the variation of YS. This means indirectly affects the isotropic hardening. Figure 6 shows the effect of parameter on the size of YS under 0 loading condition, when YS is oriented along with the equilibrium state . As the figure shows, low values of signify the role of isotropic hardening in representing the plastic strain F I G U R E 5 Anisotropy equilibrium state values of Otaniemi clay for constant stress ratio loading (data from Wheeler et al. 34 ) F I G U R E 6 Variation of yield surface size with under 0 loading condition F I G U R E 7 Variation of yield surface inclination with under 0 loading condition increments and lead to more expanded YS. On the other hand, when tends to infinity, YS almost remains unchanged in size.
Although it was found that the convergence of YS inclination to the equilibrium state might be guaranteed just for relatively high values of and there is an implicit limiting value which interrupts YS from rotating, 54 no hidden limiting value was observed in the AA1-CLAY, at least in the case of 0 loading condition YS reaches its equilibrium state for any desired > 0. But low values of might lead to unrealistically large sizes of YS. Figure 7 shows the variation of YS inclination with under 0 loading condition. As this figure shows, parameter can dictate how YS changes its inclination to lay in the direction of the prescribed equilibrium state. However, it is noteworthy that in all cases, reaching the prescribed equilibrium state is satisfied.

Plastic modulus
To establish a complete stress-strain relationship, the consistency condition should be satisfied over the YS, which guarantees that the stress state remains on the YS during the plastic flow (̇= 0). By determining the loading index from routine plasticity procedure and following a standard differentiation of the YS with regard to its state variables and substituting the hardening rules and loading index into it, the plastic modulus is obtained as It should be noted that the plastic strain components associated with the YS evolution are obtained using the Euler backward implicit integration scheme, which is described in Appendix B.

Uniqueness of CSL
From the mathematical point of view, the uniqueness of CSL is preserved if it can be proven that the ratio of mean effective stress at CSL ( ) to mean effective stress at NCL ( 0 ) is independent of Lode angle. 51 By rewriting YS equation into the stress-ratio form and setting the values of , and to the corresponding values at the critical state ( , , and , respectively), the following relationship can be obtained By considering a unique NCL in − ln space, Equation (8) shows that for constant and , there is at least a constant 0 ∕ , and by considering = , it will only depend on . Therefore, the CSL in − ln space can be specified easily for any value of 0 and the constant ratio of 0 ∕ . As the left-hand side of Equation (8) is independent of Lode angle, the uniqueness of CSL remains valid for any loading mode process in the -plane. Moreover, reaching a certain value of anisotropy equilibrium state at CSL is necessary. It means utilizing a suitable RH rule should be considered as a crucial part of the constitutive model to maintain a unique CSL. For a RH rule that relates the rate of YS rotation to only plastic volumetric strain increments (Equation 3a), reaching CSL yields ∕ = 0, and since YS does not evolve any further, it cannot be conclusively proven that an identical anisotropy equilibrium state is reached. 37 This is why several recently developed advanced constitutive models employ plastic deviatoric strain increment in their formulation to simply eliminate the ∕ part.
To schematically illustrate the uniqueness of the CSL for the simple case of triaxial loading, two undrained anisotropically consolidated compression and extension loading were simulated. As is shown in Figure 8, the evolution of YS during both compression and extension triaxial simulations leads to identical mean effective stress at CSL ( ) and YS size values ( 0 ). The sickle-shaped stress path near CSL at the extension side proves that the model continues to evolve to reach the anisotropy equilibrium state and satisfies the uniqueness of CSL.

simulation
PPS plays a decisive role in the simulation of 0 loading in isotropic constitutive models. The basic OCC and MCC models fail to predict the 0 loading precisely for normally consolidated conditions. As a remedy, using PPS with a more flexible shape has become a smart option that yields rather satisfactory results. 49,55 In the case of anisotropic constitutive models, in addition to PPS, RH rules are responsible for simulating 0 loading. The capability of the proposed RH rule guarantees that the 0 value can be adequately simulated. Figure 9 shows how the variation of parameters , and can reproduce different under the application of one-dimensional loading.

Inherent anisotropy
Before achieving a satisfactory model simulation, it is essential to properly define the initial state of soil element. In addition to preconsolidation pressure ( 0 ) and initial void ratio ( ), determination of the initial inclination of the YS (i.e., α 0 ) as a measure of inherent soil fabric anisotropy, is of great importance. Due to the sedimentation process of natural clays and also through the preparation of reconstituted samples under 0 condition, it can be assumed that the fabric of soil element alters similarly with the anisotropy equilibrium state ( 0 → ), where each 0 corresponds to a specific stress ratio 0 related to 0 via AA1-CLAY assumes the domain between -line and as the permissible range of 0 and defines it as a simple function of 0 by where for simplicity is defined as the mean of and

MODEL PARAMETERS AND THEIR CALIBRATION
In this section, the input model constants of the proposed model are discussed, and then the model response is elaborated through a step-by-step parameter calibration process using experimental results of TSC. 56 The calibration procedure of each parameter is followed by a more detailed sensitivity analysis in order to demonstrate the influence of a broader range of parameter values on the model simulations.
The model parameters and state variables required for the AA1-CLAY model can be categorized into three groups: (1) isotropic parameters including 0 , 0 (initial void ratio), , , and (Poisson's ratio) which are similar to those of MCC model, (2) YS shape parameters including , , and , and (3) anisotropy parameters including 0 , , , , , and . The evaluation of the isotropic parameters is straightforward, and they can be determined from the results of standard geotechnical tests as is the routine for MCC-based critical state models (see, e.g., Wood 57 ). The other parameters can be obtained by using the following calibration procedure.

Model calibration
The calibration procedure involves sequential tuning of YS shape and RH rule parameters to achieve the most suitable prediction of soil responses. This process starts with optimizing YS shape parameters ( , , and ). Although it is reasonable to fine-tune the shape parameters in order to capture the yield stress points more accurately, due to known uncertainties associated with determining the definitive yield stress points, 49,58,59 the authors recommend fitting the triaxial undrained isotropic and anisotropic stress paths using the variation of shape parameters. However, for the sake of completeness, the result of the best calibrated YS shape parameters can also be compared with the yield points.
The optimization of RH rule parameters ( , , , , and ) should be carried out in the second step of the parameter calibration process. For the best results, it is recommended to calibrate these parameters against one-dimensional loading tests. Nonetheless, using undrained triaxial test results can be another way to calibrate suitable values of RH rule parameters.
The calibration procedure is conducted over the experimental test data of normally consolidated TSC. The specimens were tested at 0 = 1.01 with maximum past consolidation stress of 640 kPa. The values of conventional critical state parameters are taken from Ling et al. 60 Measurements of effective stress at the axial strain of = 10% resulted in friction angle ( ) of 28 • and 33 • in compression and extension, respectively. Given the dependence of critical state ratio to (i.e., = 6 sin ∕(3 − sin ) ), was estimated as 1.1 for = 28 • . Comparison of experimental data with the model simulations ( Figure 10) suggests that = 1.3 (by considering = 33 • ) can lead to more compatible predictions and therefore, calibration of parameters is carried out by assuming = 1.3. It should be noted that a similar value for is used in the simulation of the undrained isotropic triaxial test of TSC in Hung et al. 61 4.1.1 YS shape parameters calibration Figure 11 shows the model responses with different values of YS shape parameters. By increasing the value of , the proposed model shows more distinctive vertical stress path associated with a peak strength while, in the end, all simulations approximately converge to a specific stress state ( Figure 11A). On the other hand, the variation of affects the model responses at high strain levels, and for small strain levels the model yields almost identical stress paths. It is clear from Figure 11B that higher values of produce more softening-like behaviour. The value of the shape parameter has the most significant influence on the model simulations and it can change the overall shape of the stress path. Defining this parameter independent from , brings functional capabilities to the model. The lower values of this parameter can alter the nature of model simulation from purely hardening mode to hardening-softening and purely softening ( Figure 11C). The yield loci for TSC were determined by using triaxial drained tests. Prior to shearing the samples, in-situ stress condition was reproduced, and samples were allowed to consolidate. 56 The yield data were normalized with respect to the in situ vertical stresses to eliminate the effects of different sampling depths. Figure 12 shows a practically good agreement between the experimental yield stress points, and the model YS using the best-fitted parameters in Figure 11.
To clarify the contribution of YS shape parameters on the response of AA1-CLAY, the manner in which these parameters transform the YS is also studied for the isotropic condition ( = 0). For this propose, the YS can be rewritten in the following form In which the first square bracket forms the MCC YS and the second one stands as the shape modifier. Therefore, evaluating the second square bracket defines the role of each shape parameters in adjusting the YS and their corresponding simulations in stress and strain spaces.

Parameter
According to Equation (12), parameter takes part in the transformation of YS via the last parentheses in the second square brackets, called from now on Given > 1, contracts the YS at low values of the stress ( → 0) and expands it for normally consolidated and lightly overconsolidated stress states ( → 0 ). Thus, alters the MCC YS to a teardrop-shaped with a flat cap by increasing . To better understand how affects the model simulations, let us assume that the YS and PPS are associated with the normally consolidated stress state of Figure 13. For the sake of simplicity, the isotropic condition is evaluated, but the results can be readily extended to the anisotropically consolidated soil. It is evident that in the early stages of the loading, both normals to the YS and PPS are nearly coincided and perpendicular to the stress increment. Therefore, the model experiences a nearly neutral loading situation and stress state moves along the YS ( → 0). This condition continues until the norms deviate from each other, and after that, the stress path leans to reach CSL. On the other hand, < 1 reproduces the bullet-shaped YS with an acute cap. Unlike the YSs with a bulky cap, the normals to the bullet-shaped YS start to diverge from the normal of the PPS at the beginning of the loading process, and the stress path heads to reach the CSL. Figure 14 shows the sensitivity of model simulations to the values of . The variations can be evaluated in line with the YS shapes in Figure 1A.
For the calibration purposes, since has the most striking effect on the model response from early stages of the shearing and influences the stress path from low strain levels, this parameter can be calibrated by fitting the undrained triaxial stress paths in the low strain levels.

Parameter
Parameter participates in controlling the shape of YS through the power function = ( ∕ 0 ) . is a monotonically increasing ( > 0) or decreasing ( < 0) function tending to unity when → 0 . Therefore, the variation of does not make any considerable changes in YS shape when → 0 and mainly reshapes the YS for low values of (as illustrated in Figure 1B). Considering the flow rule and consistency condition, one can realize why the stress path remains unchanged with the variation of in the early stages of shearing. The effect of this parameter on both YS and stress path becomes significant when the stress state moves away from 0 . While increasing leads to a narrower teardrop-shape YS with maximum deviatoric stress greater than that in conjugation with CSL, the stress path shows a softening behavior. Figure 15 compares how different values of can affect the response of AA1-CLAY. It is apparent that when → 1, AA1-CLAY shows a softening stress path.
The evaluation of the sensitivity analysis suggests that can be calibrated using the experimental results from undrained triaxial tests at high strain levels.

Parameter N
The term ( ∕ ) 2 as a constant ratio magnifies (when > ) or compresses (when < ) the YS. By changing the volume of the YS, this coefficient adds a desirable amount of flexibility to the model, and consequently, various forms of model response can be obtained. Considering an anisotropic state of the soil ( ) and the ratio ∕ , AA1-CLAY may reproduce various modes of simulations for normally consolidated samples from purely hardening to purely softening. Figure 16 shows a range of model simulations for different values of ∕ . The significant effect of this parameter on model responses is evidently noticeable.
After trying to calibrate parameters and to capture the overall trend of experimental test data, the model response can be fine-tuned by adjusting . It should be kept in mind that changing may lead to re-calibration of the other two shape parameters. Due to low sensitivity of model response to the variations of RH rule parameters in simulating the undrained triaxial test data (as will be discussed in the sequel), it is recommended to calibrate the RH rule parameters using the one-dimensional loading data. Moreover, since the initial inclination of the YS ( 0 ) is defined as a function of RH rule parameters (Equations 10 and 11), yield stress points can also be used to approximate or further validate these parameters' values. Calibration procedure of each RH rule parameter is discussed in the following.  Figure 17 come from the various initial anisotropy state of the soil ( 0 ), it is evident that the variation of these parameters does not have a significant influence on the triaxial test predictions. To clarify this, one can consider Figure 18 which shows the effect of parameters and on the model responses by neglecting their contribution in the changing of 0 . This figure reveals that, contrary to simulations of the 0 condition shown in Figure 9, predictions of undrained triaxial test results are less sensitive to the values of and . Therefore, unless there are no valid one-dimensional loading test results, one can ignore the flexibilities which are brought into the model by these parameters and consider them as constant values by assuming = 5 and = 2. As will be shown, the model verifications indicate that this set of parameters can lead to adequate results for various types of clay. Figure 19 shows a sensitivity analysis over a broader range of RH rule parameters' values. This analysis also proves the inconsiderable effect of these parameters on the simulation of the undrained triaxial test results.

Parameters , , and
Parameter participates in controlling the rotation of the YS for > . Thus, this parameter is applicable for specific loading conditions that continue beyond the CSL. Although using this parameter is optional, it is recommended to set ≥ 0.9 to restrict the YS from excessive rotation. Figure 19 shows how diminishes the sickle-shaped stress path by restraining the YS rotation.

Parameter
This parameter directly controls the equilibrium state of anisotropy at CSL. Since defines the intersection point of the YS and CSL at failure, it has a similar effect as shape parameter on the simulation of the undrained triaxial test data. By considering constant and values, both 0 loading and the undrained triaxial tests can be used to calibrate . Moreover, similar to and , participates in the determination of initial inclination of the YS. Therefore, yield stress points can also be used to calibrate or further validate the value of . It should be mentioned that the equilibrium inclination at CSL ( ) must be less than to have a valid YS (( 2 − 2 2 ) > 0). Figure 20 shows the simulation of undrained triaxial test data using various values. Although the effect of various 0 is concealed in these predictions ( Figure 20A), this parameter can also change the model responses by assuming a constant initial anisotropy state ( Figure 20B). Moreover, Figure 21 shows the sensitivity of model predictions against . To ignore the effect of various inclination of the YS on the model response, constant 0 is considered in these simulations.

Parameter
The remaining RH rule parameter, , controls the absolute pace of YS rotation toward its equilibrium state. Although this parameter controls how YS rotates during shearings, Figure 22 confirms that it has an insignificant effect on the F I G U R E 1 6 Sensitivity analysis for shape parameter in (A) triaxial stress and (B) stress-strain spaces overall prediction of the undrained triaxial test data. As a recommendation, values between 0.8 and 1.2 can be considered for this parameter. Table 2 summarizes the suitable tests to be used in the calibration procedure of YS and RH rule parameters. As discussed in Section 4.1, for each parameter two recommended and alternative tests are suggested. Moreover, for use in the practice with limited or no tests available, permissible parameter values are suggested.

MODEL VERIFICATION
In this section, using available experimental test results from the literature, the performance of the proposed model is studied. First, the ability of the model's YS in producing appropriate fits for yield stress points is investigated. Then, AA1-CLAY is validated using undrained triaxial test data for four different clays, including KC 62 , LCT 63 , BBC 64 , and SCC 65 . Prior to simulations, and following the procedures described in the previous section, the model parameters are first calibrated. Table 3 summarizes the model parameter values for each soil type. It is noteworthy that for all soil types, the RH rule parameters and are considered constant and equal to 5 and 2, respectively (as indicated in the previous section). Moreover, the calibration of is followed by determining the initial inclination of the YS via Equations 10 and 11.

Yield surface
As discussed earlier, AA1-CLAY employs a versatile YS equation that can capture yield points more accurately in comparison with other constitutive models with an elliptic YS. One of the primary motivations for introducing such YS is the known problems associated with capturing the experimental yield points. For the purpose of demonstration, the proposed YS is tuned to fit some of the well-documented clays in the literature (see Figure 23). 22,49 To this end, when there has not been enough data available to determine the initial inclination of the YS, the anisotropy parameter is altered until the best fit is achieved (as suggested by Wheeler et al. 36 ). Based on presented plots in Figure 23 it is reasonable to conclude that in comparison with simple anisotropic models, for example, 36-37 AA1-CLAY provides a more realistic and complete representation of the experimental yield points. Based on the evaluation of experimental data and yield stress points, Sivasithamparam and Castro 49 argued that the value of the YS shape parameter may not be constant but may depend on the degree of anisotropy. By considering this point, the YS for the isotropic condition tends to be bullet-shaped, while for anisotropic condition, an inclined ellipse or teardrop-shaped surface best represents the experimental observations. In addition to Figures 23E and 23F, triaxial simulations in the next section also confirm this hypothesis. In these experiments, the preconsolidation stress, 0 , varied between 357 and 408 kPa for isotropically consolidated tests, and 204 kPa for anisotropically consolidated tests. The experimentally observed soil responses from these tests are used here to validate the prediction capability of the developed model. Figure 24 compares the AA1-CLAY simulations with the undrained triaxial compression test results of KC. In this figure, the solid lines show the predicted results for compression tests, and the scattered symbols represent the test data. It can be seen in the figure that the model predictions are notably consistent with the experimental data in both invariant stress ( Figures 24A, 24C, 24E, 24G) and stress-strain spaces ( Figures 24B, 24D, 24F, 24H), especially for the case of normally consolidated samples. It is important to note that simulations for isotropically and anisotropically consolidated samples are obtained by using different YS shape parameter values ( Table 3). The proposed model shows good agreement with the experimental data for different consolidation conditions which supports the hypothesis of different YS shape parameter values for isotropic and anisotropic conditions, as discussed in the previous section. Considering = 0.1 and = 1.3 for isotropically consolidated sample leads to a bullet-shaped YS, while = 1.4 and = 0.4 combination results in a relatively F I G U R E 1 9 Sensitivity analysis for rotational hardening rule parameters for parameters (A, B) , (C, D) , and (E, F) bulky teardrop-shaped YS. Moreover, the excellent agreement between the model response and the experimental data for = 1.2 also confirms the suitability of calibrated YS shape parameter values (see the corresponding stress path and crossing of the YS in Figure 24A).
For the NC sample at 0 = 0.57 ( Figure 24G), the model fails to produce reasonable predictions using the calibrated YS shape parameter values for the anisotropic sample. Sivasithamparam and Castro 49 related the similarly observed mismatch to the overconsolidation of the sample due to ageing and creep effects, and they captured the best fit by considering an of 1.1. Although this assumption can be considered as a possible solution to improve the simulation quality, in this work a more realistic model simulation is achieved by only adjusting the shape parameter value = 3.6 to signify the  Figure 24G (by the dashed line). It should be noted that the initial inclination of the YS (i.e., 0 ), which is derived using Equations 10 and 11 for different values of 0 , are in a good agreement with those obtained through simulation of the 0 stress paths in Sivasithamparam and Castro. 49

Lower Cromer till (LCT)
Gens 63 performed undrained compression and extension triaxial tests on samples of LCT which were consolidated isotropically and anisotropically. The low plasticity glacial till samples were obtained from Norfolk coast in the UK where the deposit had been part of the North Sea Drift. It is a sandy clay with a liquid limit of 25%, plastic limit of 13% and a plasticity index of 12%. To resolve the variability of this natural material, Gens 63 performed his experiments on reconstituted samples which were taken from uniform specimens. Figure 25 shows the 0 loading simulation which is used to calibrate the values of RH rule parameters and to be used in triaxial test simulations. As shown in the figure, by considering constant and , the value of = 0.23 results in a very good simulation of the experimental data with almost identical stress ratio . Although has no significant effect on the triaxial response of the model, it is used to predict the one-dimensional 0 loading. Figure 26 shows the AA1-CLAY predictions of undrained isotropic and anisotropic shearing of LCT samples corresponding to experimental tests at different 0 conditions. Figure 26A shows that the stress paths predicted by the model for all values of 0 are qualitatively and quantitatively consistent with the experimental data. The model provides good  agreement with the test data in the stress-strain space, as well (see Figure 26B). Similar to KC, the isotropically consolidated sample is simulated with different YS shape parameters while all anisotropic samples are predicted with a uniform set of shape parameters. In this case, the isotropic sample is modeled using a lean bullet-shaped YS by setting = 0.2 and = 0.8, which is different from the tear-shaped YS of the anisotropically consolidated samples. Moreover, the simulations reveal the ability of AA1-CLAY in capturing the softening behavior which can be clearly observed, for example in the cases of 0 = 0.5 and 0 = 0.4 (see Figure 26). Although for the case of 0 = 0.4 the proposed model does not capture the apparent hardening behaviour preceding the softening response, but still for the pure softening part the simulation shows a good qualitative agreement with the stress-strain data (see Figure 26B). In order to better demonstrate the advantages of the proposed model, the simulations are also compared with the results of S-CLAY1 36 model. The RH parameters of S-CLAY1 36 are taken from Yang et al. 39 in which simulation of the LCT behaviour under different loading conditions were used for parameters calibrations. As can be seen in the figure, AA1-CLAY clearly outperforms S-CLAY1 36 in prediction of the stress paths and stress-strain responses at different 0 conditions during undrained triaxial shearing.

Boston blue clay (BBC)
BBC is a marine clay type that exists in a layer of the complex Boston soil deposit. This soil has a saturated unit weight of 1826-1922 kg/m 3 and its dry unit weight is between 1281 and 1393 kg/m 3 . The liquid limit, plasticity index and liquidity index of BBC are 41%, 21%, and 0.8%, respectively. A series of undrained compression and extension triaxial tests on isotropically and anisotropically normally consolidated samples of BBC were performed by Ladd and Varallyay 64 . The specimens were tested at 0 = 0.84-0.89 with 0 ranging from 273 to 785 kPa. The model parameter values, determined for simulations, are summarized in Table 3, among which the values of conventional critical state parameters are taken from Ling et al. 60 Figure 27 shows the AA1-CLAY predictions of normalized stress paths and corresponding stress-strain curves for isotropically and anisotropically consolidated samples. Figure 27A shows that the stress path predictions are in good agreement with experimental measurements for both isotropic and anisotropic samples. Unlike previous clay types, the proposed model predicts the behavior of BBC using an identical set of YS shape parameters for both isotropically and  Figure 27B).

5.2.4
Santa Clara clay (SCC) Venda Oliveira and Lemos 65 carried out a series of isotropic and anisotropic undrained triaxial tests over samples from medium plasticity clay deposits of the Santa Clara dam area in Portugal. They also performed drained triaxial tests in order to determine the yield stress points. To this end, the isotropic sample was first normally consolidated to an effective stress of 200 kPa and then unloaded down to an effective stress of 100 kPa (i.e., = 2). From this state, four drained constant stress ratio tests were conducted ( = 1.0, 2.5, 3.0, 5.0) to capture different yield points. In the case of determining the anisotropic YS, the reconstituted samples were consolidated under 0 = 0.47 condition, where 0 was obtained from the one-dimensional tests (zero radial strain) using a computer-controlled triaxial device. In the next step, the samples were unloaded to = 1.25. Similar to the isotropic condition, four constant stress ratio tests were carried out ( = -1.5, -0.5, 1.0, 3.0). Figure 28 demonstrates the undrained triaxial compression stress path of SCC for isotropic and anisotropic tests associated with their corresponding yield points. Evaluation of the isotropic triaxial test (see Figure 28A) suggests that a bulletshaped surface can best represent the YS of the isotropically consolidated samples. Using = 0.1 and = 0.4 generates a F I G U R E 2 8 Comparison of data and simulation for undrained triaxial tests on Santa Clara clay (data from Venda Oliveira and Lemos 65 ); (A) isotropic, and (B) anisotropic conditions bullet-shaped YS which reproduces admissible responses, both for the yield points and the stress path in isotropic condition. On the other hand, the anisotropic sample is simulated using a bulky teardrop-shaped YS, which offers a reasonable fit to the yield point data and leads to a realistic stress path direction (see Figure 28B).

CONCLUSION
A new adaptive anisotropic constitutive model, namely AA1-CLAY, has been developed within CSSM framework. This model can be considered as the first version and basis of a constitutive model family for natural and reconstitute clays. A flexible YS is incorporated into the model to reproduce a wide range of shapes, from simple ellipse to more advanced teardrop-shaped or bullet-shaped surfaces using two additional shape parameters. The proposed YS enables the model to capture the experimental yield points in both dry and wet sides with good accuracy. The model formulation is developed using a non-associated flow rule and an inclined elliptic plastic potential function. An added advantage of AA1-CLAY model is that based on the values of YS shape parameters, it can reduce into simpler classical constitutive models. Moreover, this model uses a new RH rule to control the rotations of yield and plastic potential surfaces. Based on experimental evidence, a transitional function is implemented into the proposed RH rule to determine the governing plastic strains and equilibrium state (which the YS tends to reach during shearing) more accurately for different constant stress ratios. In addition to the determination of the anisotropy evolution before failure, by choosing suitable model constants, the development and erasure of fabric anisotropy at critical state can also be simulated with the proposed RH rule. The numerical implementation of the model is developed via an implicit Euler Backward integration scheme. The proposed model, AA1-CLAY, has seven new parameters compared to the MCC model. A detailed sensitivity and calibration procedure has been domostrated, which can be considered as a useful guide to determine model parameter values. The model was also employed for prediction of the responses of five different isotropically and anisotropically consolidated clayey soils at different s. The comparison of model predictions against experimental data illustrated the excellent capabilities of the model in capturing different clays responses. To summarize the originality of the model, in addition to a novel yield function that has the flexibility to produce non-elliptical YSs that can capture yield points with high accuracy, it has a versatile RH rule that encapsulates the important aspects of clay anisotropy. The combination of these two novelties provides remarkable capability to outperform existing elliptical constitutive models in predicting soil peak strength, 1D compression and general stress-strain behavior both in small and large strain domains.

A C K N O W L E D G M E N T S
The financial support of the RFCS project MINRESCUE (Contract RFCS-RPJ-899518) funded by the European Commission is gratefully acknowledged.

D ATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are available from the corresponding author upon reasonable request. Yield surface By using the definition of deviatoric stress tensor ( ) and deviatoric fabric tensor ( ), the YS can be generalised as follows where, as stated in section 2.2, = . Similar to the * , * is implemented in the PPS formulation to consider the effect of Lode angle on the variation of between its corresponding values in compression ( ) and extension ( ). According to the relationship suggested by Sheng et al. 66 , the dependency of the slope of CSL to the Lode angle can be formulated as * = *  With PPS defined, now the plastic strain increments can be determined readily using the following expressions

O R C I D
where is the so-called loading index and is the plastic modulus that can be obtained from the consistency conditioṅ= 0 where and are plastic volumetric and deviatoric strains, respectively. It should be noted that the derivatives of YS and PPS with respect to and are involved with the derivatives of * and * , as these parameters are a function of stress and anisotropy state of the clay.

Hardening rules
The isotropic hardening rule is expressed identically in both invariant and general stress spaces. However, the rotational hardening rule can be rewritten in the general stress space as follows

APPENDIX B: NUMERICAL INTEGRATION
AA1-CLAY is integrated using an easy-to-implement and efficient implicit scheme called the cutting-plane algorithm. 67 In this method, the strain increments are split into elastic and plastic components. This procedure starts with the integration of the elastic predictor part, while the plastic hardening parameters remain unchanged. The trial stress increment is determined using the elastic stiffness tensor The resulting stresses are used as initial variables for plastic correction step. The plastic strain components associated with the YS evolution can be obtained using the plastic multiplier which is resulted from linearized YS formulation where the superscript ( ) represents the number of iterations before triggering the convergence criterion. It means that the iteration continues until the prescribed convergence criterion for the iterative stress integration with respect to the YS is satisfied (i.e., < 10 −7 ). Under this condition, the stress increment is calculated as It should be noted that in order to obtain stable solutions, if necessary, the size of the strain increment can be controlled through an automatic sub-stepping algorithm (see for example Sheng et al. 66 ).