AVISA: anisotropic visco-ISA model and its performance at cyclic loading

In this work, a constitutive model able to capture the strain rate dependency, small strain effects and the inherent anisotropy is proposed considering the influence of the overconsolidation ratio (OCR). Small strain effects are captured by using an extended ISA plasticity formulation (Fuentes and Triantafyllidis in Int J Numer Anal Methods Geomech 39(11):1235–1254, 2015). The strain rate dependency is reproduced by incorporating a third strain rate mechanism (in addition to the elastic and hypoplastic strain rate). A loading surface has been incorporated to define a three-dimensional (3D) overconsolidation ratio and to account for its effects on the simulations. Experimental investigations using Kaolin Clay and Lower Rhine Clay with horizontal bedding plane have shown that under undrained cycles of small strain amplitudes (<10-4\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$<10^{-4}$$\end{document}), the effective stress path in the p–q space is significantly inclined towards the left upper corner of the p-q\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p - q$$\end{document} plane. Consequently, a transversely (hypo)elastic stiffness has been successfully formulated to capture this behaviour. The performance of the model has been inspected by simulating the database of approximately 50 cyclic undrained triaxial (CUT) tests on low-plasticity Kaolin Clay (Wichtmann and Triantafyllidis) considering different deviatoric stress amplitudes, initial stress ratios, displacement rate, overconsolidation ratio and cutting direction. Furthermore, 4 CUT tests conducted on high-plasticity Lower Rhine Clay were simulated, whereby the influence of the displacement rate, as well as the deviatoric stress amplitude, has been analysed. The simulations showed a good congruence with the experimental observations.


Introduction
The prediction of the soils irreversible displacements is indispensible for the performance-based design of structures subjected to cyclic loading. For an accurate prediction of the settlement and the bearing capacity of an installed structure as, for example, a vibro-replacement with stone columns, the cyclic behaviour of the underground consisting of clay as well as its time-dependent behaviour is of high importance. The majority of analytical settlement calculation methods pertain to primary settlement only, and thus, the same value of the settlement improvement factor tends to be applied to both primary and creep settlements [28,3].
In the literature, soil anisotropy is distinguished as inherent anisotropy (as already introduced herein) or stressinduced anisotropy. The stress-induced anisotropy is produced by the recent loading history and thus results as a consequence of the applied stresses or strains. On the other side, most natural clays show anisotropic behaviour due to their mode of deposition and the elongated shape of the particles [3][4][5][6]. This anisotropy, known as inherent anisotropy, along with the time-dependent phenomena of soft soils, has practical consequences for, that is, the passive lateral thrust on piles, which has been the target of many investigations in the last decades. By creeping slopes or construction of an embankment on such formations, passive time-dependent lateral pressure may occur to the pile shaft. This may lead to deformations of the pile foundation and of the superstructure. Besides these effects, the behaviour of clays depends not only on the material state, but also on the strain amplitude. Researchers [6][7][8][9][10][11][12][13][14][15] have found that clays behave elastically only under very small strain amplitudes je 1 j\0:01 %. Under medium strain amplitudes 0:01 %\je 1 j\1 %, the so-called small strain effects, namely the stiffness increase due to reversal loading and the reduction in the plastic strain rate, take place. Finally, when the soil is sheared under very large deformations, it tends asymptotically to the critical state, which may result to a failure. Thus, in order to simulate each of these effects, the usage of an ''appropriate'' model is recommended, e.g. [14][15][16][17][18][19][20][21][22][23][24][25][26][27]. Yet, each of these effects is important for the investigations of the soft soil behaviour in urban excavations with shoring walls subjected to a low number of cycles. This cyclic loading may be attributed to the excavation shoring and backfilling processes or in case of offshore foundations with a large number of cyclic loading due to wave and wind actions or for the long-term design of stone columns as soil replacement subjected to a large number of cycles carrying traffic loading. In this paper, a single model capturing all these effects is proposed.
The paper examines the anisotropic material response of clays resulting in transverse isotropy. Only one additional parameter is required to describe this type of inherent elastic anisotropy [32]. Among all strain amplitudes, the strain rate dependency of clays is experimentally investigated. These effects are incorporated into the proposed model by an extension and modification of the (hypo)elastic stiffness tensor and the consideration of an additional strain mechanism, whose intensity depends on the material viscosity as depicted in [9]. The model is proposed under the platform of the ISA plasticity, which has shown good simulation results for of cyclic loading and captures well the small strain effects for sands. These advantages of the ISA plasticity have been extended and adjusted to the mechanical behaviour of clays.
The structure of this article is as follows: At the beginning, the formulation of the proposed model is explained. Then, some details about the numerical implementation and the integration scheme are given. The performance of the model is finally illustrated in the case of undrained triaxial cyclic simulations of two different clays: the medium plasticity Kaolin Clay and the high-plasticity Lower Rhine Clay. The experiments for the Kaolin Clay are borrowed from [36], whereas the experiments for the Lower Rhine Clay are performed by the authors as will be reported in a subsequent paper. The simulations are carefully analysed regarding the strain rate effects, different overconsolidation ratios, various initial stress ratios and mean pressures and different sample cutting directions. The samples were subjected to various deviatoric amplitudes.
The proposed model is able to describe the material behaviour of viscous and non-viscous clays under cyclic as well as monotonic loading, capturing also the present inherent anisotropy of clays. The model covers a wide range of strain amplitudes without restriction on the overconsolidation ratio (OCR). Normal consolidated as well as slightly and heavily over consolidated clays behaviour can be simulated. Simulations of tests with different initial stress ratios show the good performance of the model in capturing the induced anisotropy of Kaolin Clay. From the point of constitutive modelling, this is very interesting because the proposed model is able to simulate each effect without changing the model from non-viscous to viscous clays while keeping the same set of material parameters.

ISA plasticity for the small strain stiffness
In order to accurately describe the behaviour of granular soils under cyclic loading, in particular for a better prediction of the small strain stiffness and the shear modulus degradation curve, Niemunis and Herle [23] proposed the concept of the intergranular strain (IS). In 2015, Fuentes and Triantafyllidis [10] reformulated the intergranular strain and named it the ISA model, hereafter denoted as ISA plasticity, introducing an elastic locus in the intergranular strain space. Table 1 provides a summary of the models developed under the ISA plasticity framework. Thereby, the models in [8] and [9] are proposed to capture the specific behaviour of clays. In [9], the time dependency and rate dependency are introduced, yet the intergranular strain formulation of sand [10] is used. In the following lines, the ISA plasticity formulation for sand [9, is summarized and its extensions to capture the specific behaviour of soft soils are described and marked in red colour.
The main feature of the ISA plasticity formulation for sand is the incorporation of the elastic locus as a strain-type yield surface. It describes a hypersphere with diameter R and depends on the intergranular strain tensor h and the kinematic hardening tensor c: The evolution of the intergranular strain tensor is described in accordance with elastoplastic formulations introducing the consistency parameter _ k H and the flow direction N: whereby _ k H fulfils the consistency condition _ F H ¼ 0 and is derived to: For the plastic condition F H ¼ 0, an associated flow rule N ¼ oF H =oh and a kinematic hardening c of the yield surface are introduced as a transition state until the bounding surface (BS) is reached: The BS has been formerly proposed in [10] according to the bounding surface plasticity characteristics [4] representing it similarly to the yield surface as a hypersphere in the IS space, but with twice the diameter of the yield locus F Hb ¼ khk À R ¼ 0. Following Eq. (4), we propose a changed size of the bounding surface through the new parameter d. Hence, for sand this parameter can be set to d ¼ 1 yielding the formulation of [10]. For clays, the experimental experience gained in this work has shown that d tends to d ¼ 2. The introduction of a new size of the bounding surface has required the adjustment of the following equations compared to [10]. The kinematic hardening tensor represents the centre of the yield surface and reads: In order to fulfil the bounding constraint khk R, the hardening mechanism should evolve towards the image tensor at the BS (projection in direction of N) denoted by c b . Thereby, the hardening function is proposed to read: whereby b ¼ f ðb 0 ; rÞ is here proposed as a variable that controls the hardening rate as will be defined in the subsequent section. For isotropic conditions, it takes the value of the material parameter bðq ¼ 0Þ ¼ b 0 (Fig. 1).

Mechanical model
According to the ISA plasticity, for material states inside the IS yield surface kDek\R the response of the model is visco-elastic; thus, _ k H ¼ 0 and the intergranular strain evolves equally to the applied strain. The (visco-elastic) constitutive equation for these states reads: Once the yield surface is reached F H ¼ 0, the tensor c, describing the centre of the yield surface, evolves towards the image tensor c b located at the bounding surface according to the proposed kinematic hardening function (see Eq. 6). Upon this path, the model describes a smooth transition from the visco-elastic response (Eq. 7) to a fully mobilized response with y h ¼ 1 at the bounding condition The scalar factors y h ð0\y h \1Þ and m ð1\m\m R Þ have been proposed in [10] to interpolate between the elastic state F H \0, whereby m ¼ m R and y h ¼ 0, and the fully mobilized state, whereby m ¼ y h ¼ 1. In order to incorporate the behaviour of clays, the functions proposed in [10] have been extended to read: bounding surface yield surface Acta Geotechnica (2020) 15:2395-2413 2397 whereby h b ¼ d Á R N is the image tensor of the IS at the bounding surface. The exponent v accounts indirectly for the number of constitutive cycles the soil has experienced [25] v ¼ v 0 þ e a ðv max À v 0 ÞÞ. Laboratory tests dealing with cyclic behaviour of clays with different numbers of cycles [35] have shown that the plastic accumulation rate reduces for increasing number of consecutive cycles before reaching the critical state. This effect can be described through a new state variable e a [25], which distinguishes whether the soil is performing a few or several consecutive cycles. Its evolution rate is controlled through a material parameter C a and reads: Thus, for a few cycles e a % 0 ) v ¼ v 0 and the IS exponent proposed by [10] is recovered. e a ! 1 is reached when several consecutive cycles have been experienced and the exponent then increases to v ! v max . v 0 and v max are material parameters that can be calibrated for a few and several number of cycles, respectively. The increase in the elastic stiffness at unloading is achieved through the scalar function m ¼ m R þ ð1 À m R Þ y h , whereby m R is a material parameter introduced first in [20,23].
If the stress loops are performed at the critical state, even though the strain rate may be small jDej\R, the material response is fully mobilized. The scalar value b controls the strain amplitude required to cease the small strain effect caused by a reversal loading as depicted also in [10]. Although the former ISA models assumed b as a constant, it has been noted that the strain amplitude required to cease the small strain effects reduces for increasing stress ratios. Therefore, the present model proposes to consider b as a function of the stress ratio, according to the following relation: with the stress ratio r ¼ r Ã =p and its image at the critical state r c .
The formulation of the general model might at first sight appear high, but under some specific assumptions it reduces to some model families that possibly are simpler to understand, as outlined in Table 2. For example, if the yield surface radius goes towards R ! 0, then the ISA plasticity is disregarded and we obtain the constitutive model formulation under mobilized states.

Constitutive model under mobilized states
The fully mobilized state is reached when the intergranular strain h lies at the IS bounding surface Hence, the constitutive equation for clays at fully mobilized states reads: This relation is visco-hypoplastic as described in [32]. The evolution equation of the hypoplastic strain rate _ e hp reads: and for the viscous strain rate _ e vis : with the compression index k and the viscosity index I v .
The same direction is used for both plastic strains [32,9,30,31]: The equations governing the reference model [9] are summarized in Appendix 1. In this section, the differences between the reference model and the AVISA model are in particular presented. The new approach for the description of inherent anisotropic effects, as well as the influence of the overconsolidation ratio (OCR), is formulated with a parameter less than the reference model, and the new relation for the degree of nonlinearity Y is described in detail.  Furthermore, it should be noted that experimental studies  have shown that the viscous effects of (highly) overconsolidated clays vanish with a high exponent. For undrained cyclic loading of a normal consolidated sample with constant deviatoric amplitude, this means that the viscous effects decrease with a high exponent with each cycle. In the constitutive Eq. (7) of the AVISA model, the viscous strain is active also for elastic conditions. Yet, Eq. (14) suggests that its magnitude decreases with increasing OCR exponentially. The exponent 1=I v takes high values for 0 I v 0:5 and therefore ensures a significant reduction in the viscous strain rate for increasing OCR. For heavily overconsolidated states OCR [ 2, the viscous strain rate is almost negligible.

Inclusion of inherent anisotropy for clays
A hypoelastic stiffness has been introduced in the AVISA model: with the normalized Kronecker delta 1 and the unit tensor for symmetric tensors The barotropic bulk modulus K of clays linearly increases with increasing mean stress p as shown in [32,9]: whereby Y 0;max ¼ ðk À jÞ=ðk þ jÞ is the degree of nonlinearity deducted for normal consolidated states upon unloading [32,9]. The shear modulus G can be deduced from K and the Poisson ratio m to read: Most natural clays show anisotropic behaviour because of their mode of deposition and the elongated shape of the particles [32]. This anisotropy, known as inherent anisotropy or fabric, resulting from the deposition process tends to induce a horizontal bedding plane in the soil layer and can be regarded as a transverse isotropy [12,32]. For undrained monotonic loading, the transverse isotropy influences the slope of the effective stress path. In case of undrained cyclic loading, the described elastic anisotropy besides the effective stress slope affects the direction of the strain accumulation, whether it occurs in extension or compression. Thus, a performance-based design is not possible without describing the fabric well. However, these effects are not introduced into the reference model [9].
Furthermore, experiments show that under undrained triaxial conditions after a stress reversal first an increase in the mean pressure occurs in some clays [36] in contrast to sand where the reversal starts with the highest contractancy. This can also present a reason why most clays do not reach liquefaction p ¼ 0 kPa. Thus, also in order to capture the dilatancy/contractancy behaviour of soils well, first the accurate description of elasticity and of the inherent anisotropy is required.
For a complete transverse isotropy, five material parameters are required. However, according to the pioneer work of Graham & Houlsby [12] for the special case of anisotropic clays, a simplified transversal isotropic elasticity is deduced, whereby only three parameters are needed: the Young's modulus E ¼ E v , the Poisson ratio m ¼ m h and the anisotropic factor a obeying the relation: The index h stands for the horizontal direction and v for vertical direction. Two of these material constants are already present in the proposed model with E ¼ f ðjÞ and m. Thus, only one additional material parameter a needs to be introduced into the model. Now the question arises how can a stiffness of a model be transformed to a transversal one? After some mathematical manipulations, it becomes obvious that each hypoelastic stiffness can be transformed to a transversal one by scaling it according to the following schema [22]: m s is the unit vector along the sedimentation axis, e.g. if the sedimentation axis is horizontal, then m s ¼ f1; 0; 0g. Hence, only 1 additional parameter termed as anisotropic coefficient a compared to the reference model [9] is required.

Overconsolidation ratio OCR
Whereas in [9] an interpolation function is used for OCR introducing an additional material parameter, a more precise and comprehensive OCR formulation is proposed herein. Following Hvorslev [17], we introduce the equivalent isotropic pressure p ei ¼ expððe i0 À eÞ=kÞ. After incorporating a loading surface for the stress F b ¼ 0 where OCR¼ 1 holds and along the same lines of thoughts as described in [32], the OCR ¼ f ðp; q; e; hÞ is defined for general 3D states by: Hence, OCR is a function of the void ratio e and the stress invariants p ¼ Àtrr=3; q ¼ ffiffiffiffiffiffiffi ffi 3=2 p kr Ã k and the Lode angle h accounting for multidimensional loading paths. Thus, one parameter less is required compared to the reference model [9] as the controversial parameter n ocr introduced in [9] is omitted in Eq. (22).

Degree of nonlinearity Y
The degree of nonlinearity is proposed according to some hypoplastic requirements: • states lying at the loading surface F b ¼ 0 or at the critical state surface F c ¼ 0: experimental observations have shown that at these states the hypoplastic degree of nonlinearity takes the value Y ¼ 1.
• the degree of nonlinearity should provide information about the state of the material: the loading surface provides the boundary where Y ¼ 1 and inside the loading surface Y decreases towards Y ¼ Y 0 which is reached at the isotropic axis. This requirement is firstly termed in this work.
A suitable relation to fulfil these requirements is the following equation: For isotropic states jjrjj ! 0, relation (24) tends to: thus, the first requirement is fulfilled. States lying at the loading surface or at the critical surface are normal consolidated OCR ¼ 1, hence Y ! 1. The incorporation of the loading surface as a boundary for the degree of nonlinearity is obvious from Eq. (24). Both the degree of nonlinearity proposed in Eq. (24) and the reference model's degree of nonlinearity Y refmod [9] are illustrated in Fig. 2.
Despite the OCR independence of the Y refmod , it possesses also numerical issues approaching the isotropic axis as depicted in Fig. 2a with dashed green lines. The orange limit Y ¼ Y 0 holds for both definitions equally.

Inspection of the model's performance compared with laboratory tests
The model's performance is in this section evaluated by simulating experimental results for two different clays: Kaolin Clay from [36] and Lower Rhine Clay (experiments conducted by the authors and not published yet). The lowplasticity Kaolin Clay used for the experiments has a liquid limit of w L ¼ 47:2 %, a plastic limit of I P ¼ 12:2 % and a grain density of q s ¼ 2:590 g/cm 3 . The highly plastic Lower Rhine Clay stems from the excavation site of an opencast lignite mine and has a liquid limit of w L ¼ 56:1 %, a plastic limit of I P ¼ 34:1 % and a grain density of q s ¼ 2:675 g/cm 3 [36]. The permeability of Kaolin and Lower Rhine Clay has been determined by the authors in a permeability test in a cylinder with variable pressure head during the test resulting to k ¼ 1:3 Á 10 À9 m/s and k ¼ 3:6 Á 10 À11 m/s, respectively. Both materials have been tested in disturbed form, i.e. using reconstituted samples. The material parameters used for the simulations are listed in Table 3.  Fig. 2 Definition of the degree of nonlinearity: a contour plot of Y in p À q space, b Y in dependence on the stress ratio g (green = Y from the reference model [9], blue = Y proposed in Eq. (24)) (color figure online)

Simulations with Kaolin Clay
In several test series performed by Wichtmann & Triantafyllidis [36], the stress amplitude q 0 , the displacement rate _ s, the initial mean pressure p 0 , the OCR, the initial stress ratio g 0 ¼ q 0 =p 0 and the cutting direction of the samples have been varied as described in [36,37]. Each of these tests has been simulated with the AVISA model and showed a good agreement with the experimental results (see the following sections). Only one material parameter set listed in Table 3 has been used for all simulations. The notation of the tests with CXX, e.g. C02 is the same as in [36], so that missing information can be depicted from [36].  The tests C01-C08 was performed with different deviatoric stress amplitudes q ampl ¼ f30; 40; 45; 50; 60; 70g kPa. The initial conditions OCR 0 ¼ 1; p 0 ¼ 200 kPa, g 0 ¼ 0; _ s ¼ 0:1 mm/min were the same for all samples. In the simulations, the same number of cycles as documented in [36] has been applied. For the initialization of the OCR, the simulations used the initial void ratio e 0 ¼ e i0 À k lnðp 0 =ð1 kPa Þ OCR Þ.
The simulated effective stress paths (red) and the experimental results (blue) are given in Fig. 3 for tests C02, C04-C08. It is evident that the effective stress path before cyclic mobility is significantly inclined towards the left upper corner of the p À q plane [36]. This observation is reproduced very well with the AVISA model through the anisotropic coefficient of a ¼ 1:8 resulting in E h ¼ 3:24 Á E v . The reference model [9] could not reproduce the described slope of the effective stress path. The inherent anisotropy effects also the axial strain path, whereby the sample first accumulates in extension. Due to the p À q inclination, first the compression triaxial critical state line is reached, and thus, the accumulation of the axial strain changes to compression as demonstrated in Fig. 4. The AVISA model captures well all of these features of the material behaviour, whereas the reference model [9] accumulated only in extension. For a performance-based design of an, for example, flooded opencast mine under earthquake loading, the accurate description of these effects is indispensable.
The tests C05 and C06 have been performed with the same deviatoric amplitude q ampl ¼ 50 kPa but with different initial void ratios e 0 ¼ 1:145 and e 0 ¼ 1:252, respectively. As reported in [36], lower initial void ratios due to higher overconsolidation ratios lead to a higher number of cycles to reach cyclic mobility or the failure criterion of je 1 j ¼ 10 % used in the experiments. The simulations are of course completely in agreement with the described behaviour, presented also in Fig. 5a in terms of accumulated pore water pressure against the number of cycles N and Fig. 5b, whereas the accumulated axial strain against N is illustrated. This behaviour is observed also by other authors [35,14,16]. Moreover, Fig. 3d shows a better agreement for the first cycle between experiment and simulation. Due to the experimental data scatter of the initial void ratio, the first cycle in some simulations does not coincide exactly with the experiment.
The tests pair C03 and C04 were performed with the same amplitude of q ampl ¼ 45 kPa and showed for the higher initial void ratio of e 0 ¼ 1:224 corresponding to the test C03 an increasing number of cycles up to failure in opposition to the above-described effect for tests C05 and C06. The authors addressed this occurrence to the natural scatter of experimental data [36]. The AVISA simulations showed a good agreement with the test C04 as illustrated in Fig. 3b.
The last cycles at cyclic mobility show a prevented dilatancy at loading and a contractant behaviour after an elastic regime is overcome for unloading. This behaviour, resulting also in an axial strain ''hysteresis'' in Fig. 4, is not accounted by the model as for its description the precise description of the energetical correct hyperelasticity of the Kaolin is necessary. The authors are currently working in this direction. Yet it is of high importance to note that at the stress reversal clays first exhibit an increase in the mean pressure in contrast to sand, whereas the reversal starts with the highest contractancy. This can present a reason why most clays do not reach liquefaction p ¼ 0 kPa. Therefore, the question arises if there is an angle for the cutting direction of the sample at which the clay reaches a total degradation of the mean stress similarly to sand. In Sect. 4.1.7, we will notice a reduction in the relaxation of the mean stress compared to vertically cut samples for the same initial and boundary conditions.

Variation of displacement rate _ s
Six tests with variation of displacement rate between _ s ¼ f0:01; 0:05; 0:5; 0:01; 0:05; 0:5g mm/min have been documented in [36] and simulated with the AVISA model.  Isotropically normal consolidated samples with an initial mean stress of p 0 ¼ 200 kPa were tested with a deviatoric amplitude of q ampl ¼ 45 kPa (C09-C12) and q ampl ¼ 50 kPa (C13-C15). Figure 6 shows the results in the effective stress plane, whereas the corresponding deviatoric stress-axial strain relations are, for example for, tests C13 and C14 illustrated in Fig. 7. In order to distinguish between the effects of stress amplitude and displacement rate, either tests  C09-C12 or C13-C15 are compared. Due to the low plasticity of Kaolin, an effect of the displacement rate can be seen rather comparing the lower rates _ s ¼ 0:01 and 0.02 mm/min with a number of cycles to failure N f ¼ 81 and 58 (C09-C12) or N f ¼ 24 and 40 (C13-C15). The increasing number of cycles to failure with decreasing frequency has been observed by other authors too . In agreement with this tendency is an acceleration of the accumulated pore water pressure and of the accumulated axial strain for lower displacement rates as depicted in Fig. 8. The simulations show a good agreement with all effects shown by the experiments. The frequency dependence of the model is established through the viscous strain e vis with the evolution relation given in Eq. (14). The inherent anisotropy resulting in a compressive accumulation of the axial strain and an inclination of the effective stress path to the upper left of the p À q diagram is very well captured by the simulations with the AVISA model. Responsible for its simulation is the transversal (hypo)elasticity with the anisotropic coefficient of a ¼ 1:8. The reference model [9] is not able to capture the inherent anisotropy would result in a vertical effective stress path against the experimental observations.

Variation of initial mean pressure p 0
Six isotropically normal consolidated samples with variation of the initial mean pressure p 0 ¼ f50; 75 100; 125; 150; 250g kPa have been tested with deviatoric amplitudes resulting in a constant amplitudepressure ratio of f ¼ q ampl =p 0 ¼ 0:2 and a displacement rate of _ s ¼ 0:1 mm/min. A comparison between test C21 with p 0 ¼ 250 kPa presented in Fig. 9f with N f ¼ 995 and C5 with p 0 ¼ 200 kPa presented in Fig. 3c with N f ¼ 68 shows the significant influence of the initial mean stress. A decrease of the mean pressure leads to an acceleration of the accumulated pore water pressure and of the residual strain. This behaviour is also captured with the model simulations illustrated in Figs. 9 and 10. No simulations with different initial mean pressures are presented in [9] with the reference model. However, since the radius of the bounding surface is randomly 2 times larger than that of the flow surface in the reference model [9], it is to be expected that laboratory tests with different p 0 cannot be accurately reproduced with the same parameter set. The proposed model introduces the parameter d in Cap. 2 in order to overcome these disadvantages and to achieve a good performance as depicted in Figs. 9 and 10.
The mismatch in the first cycle of the effective stress paths can be attributed to the fluctuation of the experimental initial void ratio.

Variation of initial overconsolidation ratio OCR 0
The influence of the initial overconsolidation ratio was studied in 2 tests with OCR 0 ¼ f1:5; 2g. The initial conditions were chosen to p 0 ¼ 100 kPa, q ampl ¼ 30 kPa and _ s ¼ 0:1 mm/min. The expected dilative response in the first cycle of undrained cyclic loading of overconsolidated samples is recognizable in the negative accumulated pore water pressure at N ¼ 1 in Fig. 11a. Afterwards, the material response is contractive _ p\0 ) _ u acc [ 0. The number of cycles to failure N f ¼ f116; 531g grows significantly with the OCR 0 ¼ f1:5; 2:5g. The rate of e acc 1 and u acc decreases with OCR 0 resulting in a stiffer response of the material for higher OCR 0 . All these effects are properly described by the models simulations, solid lines in Fig. 11.

Variation of initial stress ratio g 0
Normal consolidated samples with initial mean pressure of p ¼ 200 kPa have been anisotropically consolidated to different deviatoric stresses leading to various initial stress ratios in compression g 0 ¼ f0:375; 0:25; 0:125g and in extension g 0 ¼ fÀ0:125; À0:25; À0:5g. Then, cyclic undrained shearing has been performed with q ampl ¼ 30 kPa and displacement rate of _ s ¼ 0:1 mm/min. The effective stress paths are presented in Fig. 12a-c, whereby the slope and the relaxation of the mean pressure are properly reproduced with the model. Both the experiments and the simulations, showed that the failure criterion is reached due to excessive axial strain accumulation and not due to large strain amplitudes as observed for isotropically consolidated samples (compare Figs. 4 with  13). Furthermore, the accumulation of the pore water pressure almost stabilized after a certain number of cycles, and consequently, a complete relaxation of p is not expected.

(b) (a)
The simulation in Fig. 12c shows a partially flow of the mean pressure when approaching the extension critical state line. This behaviour is not evidenced by these experiments, but it has been documented in some other works [18,33].

Variation of stress amplitude q ampl at anisotropic consolidated samples
The performance of the model is tested at anisotropic initial stress ratios also for different stress amplitudes q ampl ¼ f20; 25; 30; 35; 40g kPa (notation in [36] C33, C34, C25, C35 and C36, respectively). The stress-strain relation is exemplary illustrated in Fig. 13 for C35, whereby the axial strain amplitude is observed to remain constant during the cycles resulting in a constant but considerable rate of the irreversible axial strain. All these effects are captured by the simulation with the AVISA model. For other deviatoric amplitude was observed a similar behaviour as depicted in Fig. 14b. Thereby, the accumulated axial strain is plotted versus the number of cycles. Obviously, the number of cycles to failure decreases with increasing amplitude. The existence of a lower limit of the number of cycles to failure, for example N f ¼ 19, is indicated, as for both tests C35 and C36 19 cycles are passed through until the failure criterion of je 1 j ¼ 10 % is reached. The accumulated pore water pressure stabilizes after a certain number of cycles, and it indicates the existence of a threshold value of u acc ¼ 0:5 p 0 for each test. The simulations are in very good agreement with the experimental evidence presented by [36].

Samples cut out in horizontal direction
The inherent anisotropy incorporated in the model requires the calibration of the material parameter a termed as anisotropic coefficient and the definition of the unit vector along the sedimentation axis m s (not to be mixed up with the tensorial flow rule denoted with m). The calibration of a is in detail described in [32] and once more summarized in Appendix 2. The vector m depends on the sedimentation axis of the material, and we define for a vertical cutting direction m s ¼ f1; 0; 0g and consequently for a horizontal cutting direction m s ¼ f0; 0; 1g. In [36], despite samples cut out in vertical direction, three tests 2 Â g 0 ¼ 0 (C40 and C41) and 1 Â g 0 ¼ 0:5 (C42) on samples cut out in horizontal direction have been performed as well.
The experimental results and simulations with the AVISA model are presented in Figs. 15 and 16. The lower initial void ratio of the test C41 implies a lightly overconsolidated initial state with OCR ¼ exp½ð1:76 À 1:07Þ= 0:13=200 ¼ 1:48. As a result, a stiffer response of the material is evident in the p À q plane as illustrated in Fig. 15 and a decay of the accumulation rate of the axial strain can be observed in Fig. 16b. These observations are of course mirrored in the accumulated pore water pressure too. Moreover, the inclination of the effective stress path is opposite to the inclination of samples cut out in vertical direction (compare Figs. 3 with 15). The final effective stress loop of horizontal samples is practically mirrored along the p-axis. Due to the higher elasticity, the rates of the accumulation of pore water pressure as well as of the axial strain rate are much lower for the horizontally cut samples compared to the vertical ones (compare Figs. 5 with 16). Consequently, for a given amplitude or another boundary condition the horizontal samples can withstand a much higher number of cycles to failure than the vertical ones. All these observations are very well reproduced with the AVISA model. Furthermore, the final relaxation of the mean pressure at the failure criterion of je 1 j ¼ 10 % is lower for horizontal samples. This occurrence gives a hint to the answer of the question if a specific angle for the cutting direction exists, at which a total relaxation near to liquefaction p ! 0 would be experienced.
Of course, the reference model [9] is not able to reproduce a distinction between different cutting directions, because of the absence of inherent anisotropy in its mathematical formulation.

Simulations with Lower Rhine Clay
The experiments with the high-plasticity Lower Rhine Clay have been performed by the authors and Wichtmann in IBF, KIT and will be published once a complete database is finished. In the following, some results of undrained cyclic triaxial tests with different deviatoric amplitudes as well as with variation of the displacement rate are presented. By testing a specified stress amplitude q ampl , the loading direction was changed once the defined deviatoric stress q was reached. The cyclic loading was stopped when a certain value of axial strain j e 1 j¼ 10 % was reached, as done also for the Kaolin tests [36]. The displacement rates were chosen low enough to guarantee a homogeneous field of pore pressure within the samples. The material parameters used for all the simulations are depicted in Table 3. Table 4 provides a short overview of the experiments performed on Lower Rhine Clay with the specifications of initial conditions in terms of: mean pressure p 0 , overconsolidation ratio OCR 0 and deviatoric stress q 0 . The deviatoric amplitude q ampl , the displacement rate _ s and the number of cycles to failure N f which the experiments were subjected to are also documented in Table 4. Figure 17 shows the results for the test RBT-4 with a deviatoric amplitude of q ampl ¼60 kPa and N f % 1:5 until the failure criterion of je 1 j ¼ 10 % is reached. For the same amplitude but with a 2.5 times faster displacement rate, the number of cycles to reach an axial strain of 10 % decreases to approx. 0.5; see Fig. 18. The decay of the number of cycles to failure with increasing frequency of the same amplitude has been observed by other authors too and in accordance with this tendency is an acceleration of the accumulated pore water pressure to observe. Thus, the accumulated axial strain increases for lower displacement rates. The simulations show a good agreement with all effects shown by the experiments. The frequency dependence of the model is established through the viscous strain. Figures 19 and 20 present the experimental and simulated behaviour of Lower Rhine Clays for lower deviatoric amplitudes corresponding to q ampl ¼ 45 and 40 kPa. As expected, the number of cycles leading to failure increases with decreasing amplitude. Among other effects, also the inherent anisotropy resulting in a compressive accumulation of the axial strain is very well captured by the simulations with the AVISA model. No simulations with the reference model are shown in [9] with cyclic tests of Lower Rhine Clay. However, it is obvious that the

Concluding remarks
The present work proposes a constitutive model with some interesting capabilities, such as the simulation of the strain rate and time dependency and the incorporation of small strain and of inherent anisotropic effects, permitting the user to evaluate the influence of each effect on a simulation, without changing the model and keeping the same set of material parameters. The small strain effects of clays have been simulated through an extended ISA plasticity formulation (first version of ISA plasticity has been proposed for sand by Fuentes & Triantafyllidis [10]). The rate dependency is captured by a third strain mechanism (in addition to the total and the hypoplastic strain) that is dependent on the stress and the void ratio. The loading surface has been well defined in order to provide a clean and complete formulation for the overconsolidation ratio. The degree of nonlinearity presented in the reference model [9] turned to be unsuitable as it does not provide information about the preloading state of the material. For this purpose, it has been reformulated depending on the OCR. Some investigations regarding the inherent anisotropy have shown that the effective stress path before cyclic mobility is significantly inclined towards the left upper corner of the p À q plane as documented also in [36] for Kaolin Clay and that the inherent anisotropy of clays can be described through the transversal isotropy as theoretically described also in the pioneer work of Graham & Houlsby [12]. Due to the p À q inclination, first the compression triaxial critical state line is reached, and thus, the accumulation of the axial strain changes from extension to compression. For a performance-based design of an, for example, flooded opencast mine under earthquake loading, the accurate description of these effects is indispensable. This holds not only for the low-plasticity Kaolin Clay as documented in [36] but also for the high-plasticity Lower Rhine Clay evaluated in this work. Consequently, the hypoelastic stiffness has been transformed to a transversal hypoelastic stiffness. The extensions and novelties compared to the reference model [9] are in detail presented in Appendix 1.
Furthermore, the last cycles at cyclic mobility show a prevented dilatancy at loading and a contractant behaviour after an elastic regime is overcome for unloading. It is of high importance to note that at stress reversal both clays investigated in this paper first exhibit an increase in the mean pressure in contrast to sand whereas the reversal starts with the highest contractancy. This can present a reason why most clays do not reach liquefaction p ¼ 0 kPa. Particularly, the fact that the final relaxation of the mean pressure at the failure criterion of je 1 j ¼ 10% is lower for horizontally cut samples than vertically cut ones gives a hint to the answer of the question if a specific angle for the cutting direction exists, at which a total relaxation near to liquefaction p ! 0 would be experienced.
The performance of the model has been inspected by simulating the database with about 50 cyclic undrained triaxial (CUT) tests on the low-plasticity Kaolin Clay published by Wichtmann & Triantafyllidis. The Kaolin database includes variation of the deviatoric stress amplitude, initial stress ratio, displacement rate, overconsolidation ratio and cutting direction. Furthermore, 4 CUT tests performed on the high-plasticity Lower Rhine Clay have been simulated. These experiments have been performed by the authors with the same device as [36], whereby the influences of the displacement rate as well as the deviatoric amplitude have been studied. The simulations showed to be in a very good agreement with the experimental observations.
The research has shown that the loops at the cyclic mobility represent an interaction between inherent anisotropy and dilatancy. The later one is not considered in the present work and will be a matter of future research.
Acknowledgements Open Access funding provided by Projekt DEAL.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/. Appendix 1: Summary of constitutive relations of the reference and the proposed model Table 5 provides a summary of the relations constituting the reference model compared to the constitutive equations of the in this paper proposed AVISA model (in red the reformulations or extensions of the AVISA model).

Appendix 2: Short overview of the calibration of AVISA material parameters
The AVISA model requires the calibration of 15 parameters (one more than the reference model [9]) listed in Table 3. However, the ''controversial'' n OCR parameter termed as viscous exponent is eliminated due to the explicit solution of the loading surface according to the void ratio (see Sect. 3.2). Despite the bounding surface radius d, all parameters have been used in former works [32,9,[23][24][25][26][27][28][29][30][31]25] and will be briefly summarized herein by separating them according to their role within the model: • transversal (hypo)elasticity requires the determination of 4 material parameters, whereby the compression index k can be calibrated upon an compression path with a constant rate of strain and the swelling index j is defined as the slope of the corresponding unloading path. Two remaining parameters can be defined in conjunction with each other, as they both influence the slope of the effective stress path in stress space and the shear modulus. Thus, for the Poisson ratio m h and the anisotropic coefficient a, the effective stress slopes Dq=Dp of two monotonic or cyclic undrained triaxial tests with different cutting directions need to be measured, and then, the hypoelastic formulation can be used: with De v ¼ 0 in undrained tests. Along these lines of thoughts, the authors provided in [32] a Mathematica code presenting a very simple and self-explanatory method for the calibration of m h and a. • the intergranular strain requires the calibration of 7 material parameters whereby the stiffness factor m R , the The parameter R corresponds to the amplitude of the threshold strain which encloses the elastic locus of the material, and when results of the RC test are not available, a value of R ¼ 10 À4 is recommended for the elastic behaviour of clays. The parameter b 0 controls the strain amplitude required to reach the mobilized states. A relation for b 0 was provided in [10]. The minimum v 0 and maximum v max exponent control the shape of the secant shear modulus degradation curve and should be calibrated with tests including a few and a considerably increasing number of cycles, respectively; for further details, see also [25]. The accumulation rate factor C a controls how fast the accumulation is produced during the cycles and can be calibrated through trial and error procedure of an undrained triaxial test. The ratio d of the radius of the bounding surface and the yield surface d should be calibrated through a RC test in the transition between the viscoelastic and the visco-hypoplastic behaviour defining the end of this transition area as depicted in Fig. 21. • the critical and loading surface introduces 1 and 2 additional material parameters, respectively. The critical state slope in triaxial compression M c can be calibrated through a triaxial test with samples sheared up to high strains ensuring that the critical state (e ¼ e c and g ¼ g c ) is reached. The maximum void ratio e i0 corresponds to the maximum void ratio at the reference pressure p ref ¼ 1 kPa at the virgin isotropic consolidated line OCR ¼ 1. For its calibration, the most suitable test is an isotropic compression test with constant rate of strain. In the absence of this test, an oedometer test can be used instead as explained in [32,9]. The loading surface factor f b0 controls the peak strength of the overconsolidated soil. Thus, it can be calibrated using triaxial data of initially overconsolidated samples. • the viscosity requires the calibration of only one parameter denoted as viscosity index I v . It describes the intensity of creep for normal consolidated states and can be determined through a creep stage or consolidation process evaluating the secondary settlement curve.
Other methods as the consideration of two different isotachs can be used as well [20,32,9,8].