Comparison of two small-strain concepts: ISA and intergranular strain applied to barodesy

The intergranular strain concept (IGS) and intergranular strain anisotropy formulation (ISA) are state of the art extensions to describe small-strain effects. The main conceptional difference between ISA and IGS is the purely elastic strain range introduced by ISA. In addition, the ISA formulation used in this article includes an additional state variable in order to reduce accumulation effects for cyclic loading with a larger number of repetitive cycles. Barodesy is enhanced here with ISA to improve its small-strain predictions. The performance of this new model is compared with barodesy enhanced with IGS. It turned out that the small-strain extensions do not negatively influence predictions under monotonic loading. Differences between ISA and ISG are only remarkable for very small-strain cycles and even there they are negligible for certain parameter values. Supplementary Information The online version contains supplementary material available at 10.1007/s11440-022-01454-3.


Introduction
It is essential to understand and constitutively describe the mechanical behaviour of fine-grained soils under cyclic loading for many applications such as offshore wind turbines subjected to wave cycles and wind cycles, geotechnical structures subjected to repetitive traffic and highspeed train loads, or even earthquake loading.
In finite element simulations, it is important to use a constitutive model that includes stress-strain relations that are appropriate for the given problem. In order to accurately address the small-strain stiffness and reproduce the soil behaviour under cyclic loading [36] using hypoplastic models, Niemunis & Herle in 1997 [24,26] proposed the concept of the intergranular strain (IGS). This extension aimed to remedy the hypoplastic saw-tooth-like paths (ratcheting), which yields an excessive accumulation of deformation predicted for small stress cycles. Even though the IGS approach considerably reduces ratcheting, it does not vanish completely. Furthermore, the strain amplitude at cyclic mobility during undrained cyclic shearing grows too slow. [17] formulated the hypoplastic model with IGS for sand under the framework of anisotropic critical state theory (ACST) introducing an evolving fabric tensor as a remedy for these shortcomings when dealing with granular materials. In this course, [39] addressed the description of inherent anisotropy of granular soils in hypoplastic models using the ACST. [14] included strength anisotropy into the hypoplastic model [19] by rotating the asymptotic state boundary surface.
Fuentes and Triantafyllidis [7] reformulated the intergranular strain concept in 2015 and named it ISA-model, often denoted also as ISA-plasticity. This new approach introduces an elastoplastic intergranular strain, and hence, it includes a yield surface in the intergranular strain space. It was first introduced for granular materials using projections in the stress space for the mechanical model. for sands [27] was determined by the coupling between the ISA concept of the intergranular strain and the well-known hypoplastic model for granular materials developed by von Wolffersdorff [38]. The latter one loses the ability to account for cyclic mobility effects, and hence, liquefaction phenomena are not reproduced. In 2019 [10], ISA has been extended by introducing an additional state variable that permits the detection of cyclic mobility paths. For fine-grained soils, a first model has been developed in 2015 [8] and extended in 2017 [9,29] by introducing a time and rate dependency according to [33,34]. In both models for fine-grained soils, the intergranular strain formulation of sand [27] is used. Tafili and Triantafyllidis [28,31] enhanced the ISA-plasticity to capture the specific behaviour of soft soils and introduced also the inherent anisotropy effects [30,32] into the mechanical model formulation (AVISA model). In 2020, the ISA-plasticity has been coupled with the hypoplastic model for finegrained soils developed by Mašín [11]. This model also accounts for the inherent anisotropy of fine-grained soils, but does not introduce rate-and time-dependent effects.
Barodesy is a constitutive model, which is attractive for engineering purposes due to the reduced number of parameters and also the simplicity of their determination. A fundamental drawback of the model is, however, the inability of describing the behaviour of soils under cyclic loading. Therefore, in this article, barodesy for clay [22] is extended with ISA according to [27], hereafter denoted as ISA-B. Comparison is made with the IGS [26] extension for barodesy [1], hereafter denoted as IGS-B. We decided to use the most common versions of ISA [27] and IGS [26] to investigate conceptual differences. The main difference between ISA and IGS is that only ISA includes a purely elastic strain range. In addition, the ISA version used in this article includes an additional state variable in order to reduce accumulation effects under cyclic loading with a large number of repetitive cycles.
In this article, the similarities and differences between ISA-B and IGS-B are discussed. In addition, the models are compared with (cyclic) experiments of Kaolin and Lower Rhine clay.

Barodesy
The constitutive model barodesy differs fundamentally from conventional elastoplastic models: it does not distinguish between elastic and plastic deformations and does not require expressions such as a yield function, a plastic potential or a flow rule. The stress rate is formulated as a tensorial function of the current stress, stretching and the void ratio, i.e. T ¼ hðT; D; eÞ. Thus, barodesy has certain similarities with hypoplasticity [19,38].
Basic concepts from Critical State Soil Mechanics are included in the mathematical formulation of barodesy [22]: (i) Critical stress states of barodesy virtually coincide with the Matsuoka-Nakai failure criterion [6]. (ii) A stress-dependent critical void ratio e c (CSL according to [18]) in the p 0 -e plot enables to distinguish between normally to slightly overconsolidated soil (e [ e c ) and highly overconsolidated soil (e\e c ). (iii) The isotropic normal compression line (NCL) according to [4] defines normally consolidated states.

Intergranular strain concept (IGS)
The intergranular strain concept (IGS) was introduced by Niemunis and Herle [26] to improve the description of soil behaviour under cyclic loading with hypoplastic models and to capture the effects of small-strain stiffness. The short-term deformation history is stored as an additional tensorial state variable, the so-called intergranular strain h. After a change in the deformation direction, the stiffness of the hypoelastic (linear) part of the hypoplastic model is used, however, increased by a scalar factor. For a constant stretching after the change in the deformation direction, the ISA-HP Clay (2020)

Clay
Inclusion of inherent anisotropy [11] IGS performs a reduction of the increased stiffness and an interpolation between the linear part and the full hypoplastic model. Because of its mathematical formulation, the original IGS by Niemunis and Herle [26] is not directly applicable to non-hypoplastic models. To overcome this, Bode et al. [1] proposed a formulation of the IGS to allow for an application to barodesy. Details of this formulation and also of the general IGS formulation can be found in [1, 2, 24, 26, e.g.].

ISA-plasticity for the small-strain stiffness
In the following lines, the ISA-plasticity formulation proposed by Poblete et al. [27] is briefly summarized using the notation depicted in Sect. A. Some explanations about the state variables and parameters responsible for the reproduction of cumulative soil behaviour are provided. For further details, we refer to the original works published in [7,9,27].
The main feature of the ISA-plasticity formulation is the incorporation of the elastic locus, see Fig. 1c, as a straintype yield surface. It describes a hypersphere with diameter R as a material parameter 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 for _ F H ¼ 0 and the associ- Once the yield surface is reached, its centre evolves towards the bounding surface (BS) governed by the kinematic hardening state variable c and its evolution equation whereby b is a material parameter controlling the hardening rate.
According to the ISA-plasticity, the BS is reached after applying stretching in a constant direction D 0 . In order to reduce the plastic strains between this ''fully mobilized'' state and the elastic condition (e.g. due to an unloading process), the scalar function y h has been introduced whereby h b ¼ R N is the image tensor of h at the bounding surface.
According to physical experiments, e.g. [28,35,37], the cumulative behaviour of the soil reduces for increasing number of consecutive cycles before reaching the critical state. To describe this effect, Poblete et al. [27] introduced the so-called internal variable for cyclic history 0 e a 1 to the former ISA-model [7], which indicates whether the soil has been objected to a few or several consecutive cycles. For reconstituted samples, this variable can be initialized to e a;0 ¼ 0, because the sample has not been previously subjected to cyclic loading. As this work deals with the constitutive description of reconstituted soil, e a;0 was initialized to zero in all simulations. For in situ conditions, the initialization is to some extent controversial because it often cannot be established whether the soil has been previously subjected to cyclic loading. Therefore, we recommend initializing the variable to zero while simulating the entire history of previous loading. The evolution rate of the internal variable, which will then save the information about the altering of the soil subjected to consecutive cycles, is controlled through a material parameter C a and reads For a few number of cycles or for monotonic loading y h ! 1, hence e a vanishes, otherwise, for a large number of consecutive cycles y h ! 0 ) e a ! 1. This feature is now used to intensify the effect of the intergranular strain at larger number of cycles through the modification of the exponent v in Eq. 5. As y h as well as the plastic accumulation rate decrease with increasing v, the relation v ¼ v 0 þ e a ðv max À v 0 Þ has been introduced in [27]. Hence, the minimum value v ¼ v 0 resulting in maximal plastic strain accumulation is reached for a few number of cycles or for monotonic loading. On the other hand, a large number of cycles leads to the maximum value v ¼ v max , which in turn reduces the plastic accumulation rate. Hence, v 0 and v max are material parameters that can be calibrated for a few and many number of cycles, respectively.
To reproduce the stiffness increase due to changes in the direction of the strain path, the scalar function m ¼ m R þ ð1 À m R Þ y h is included into the ISA-plasticity, where m R is a material parameter. Note the dependence of the stiffness increase factor on v through y h .
A more recent development by [5] applied to clay hypoplasticity [19] contains a similar improvement for IGS to predict cyclic loading as introduced for ISA by [27].

Implementation of ISA in barodesy
In this section, the ISA implementation of hypoplasticity [27] will be extended with the ideas used in the IGS implementation of barodesy [1]. The basic formulation of hypoplasticity reads In this formulation, the term linear in D (i.e. L H : D) represents an incrementally elastic response T el . The other term is not linear in D and represents the non-elastic response.
The stress rate of the ISA hypoplastic model is given by with the material stiffness matrix Combining (8) and (9), using D ¼ kDkD 0 and rearranging N H N : Using the decomposition of the objective stress rate (7) in an elastic and a non-elastic part yields This formulation is now ready to be used with barodesy. Barodesy does not allow for a computation of L Baro and N Baro . However, it is easy to compute the incrementally elastic response, which delivers equal stress rates for reversal stretchings and is therefore the odd part of the constitutive model The nonlinear stress rate is then the rest which is the even part of the constitutive relation. This approach is called internal elastic model in [1], where also an alternative and a bit simpler approach by employing a hypoelastic model is proposed, which is used for the computations in this article. This so-called external elastic model approach uses the incrementally elastic stress rate The parameters of the pressure depending elastic stiffness must be calibrated such that the hypoelastic model is consistent with barodesy, see [1] and ''Appendix C''. The nonlinear stress rate is the difference of the stress response of barodesy T m and the elastic stress rate For certain conditions, the formulation of the model simplifies to certain forms, which are summarized in Table 2.

Analysis of the qualitative behaviour of the models
The effect of the different mathematical formulations of ISA and IGS is investigated in this section and reveals some differences and similarities between the predictions of the models. The element tests for these investigations have been carried out using the incremental driver by Niemunis [25] available on www.SoilModels.com [13]. Calibration A default set of clay parameters according to Tables 3  and 4 is used. For the investigations, v was set equal to v 0 ¼ v max and for IGS 1 m R was set to m T , in order to achieve results as comparable as possible. 180 strain path reversal A 180 strain path reversal implies isochoric triaxial reloading subsequent to isochoric triaxial unloading. Figure 1 shows the simulation of a CU test, where a 180 strain path reversal is performed at fully mobilized intergranular strain jjhjj ¼ R. Thus, the intergranular strain tensor is initialized with A and for ISA the kinematic hardening tensor is initialized with c 0 ¼ h 0 =2. The initial stress state is p 0 0 Á I ¼ 200 Á I kPa with an initial overconsolidation ratio 3 (OCR) of 2.5. The stress-strain behaviour of IGS and ISA virtually coincides, see Fig. 1a. The shear stiffness follows from the elastic model and is constant until The elastic response, whose stiffness is governed by m R , is followed by a transition between the elastic model and barodesy. Figure 1c shows the Rendulic plane of the intergranular strain h space and the stretching space D. The shear stiffness is constant as long as jjhjj\R, followed by an interpolation between elasticity and barodesy. When the bounding surface (larger bubble) is reached and the direction of the stretching coincides with the direction of the intergranular strain, the material behaviour is purely barodetic.  Table 3 Parameters for barodesy, default parameter set for clay [19] Table 4 Parameters for the qualitative analyses with IGS and ISA, Kaolin clay the stiffness degradation of IGS and ISA slightly differs. For the simulations in Fig. 2 the closer is the agreement between IGS-B and ISA-B.
Note that also b influences the stiffness decay with ongoing shear strain. The larger b is, the higher is the decay with ongoing shear strain.
A 90 strain path reversal occurs, if undrained triaxial compression is performed after isotropic compression. Figure 3 shows the simulation of a CU test, following a strain path reversal of 90 , starting at fully mobilized intergranular strain with h 0 ¼ À1= The initial stress state is p 0 Á I ¼ 200 Á I kPa again with an initial OCR of 2.5. The stress-strain response of IGS-B and ISA-B is similar, see Fig. 3a, with small differences, better visible in Fig. 3b and c. For IGS: If the angle between h and D is between 0 and 90 , h stays fully mobilized (jjhjj ¼ R) and rotates towards the direction of D [20,26], see Fig. 4c.
For ISA-B: The intergranular strain h moves inside the bounding surface on the yield surface towards the bounding surface. As for IGS-B, the stiffness is obtained through interpolation between the elastic model and barodesy if the angle between h and D is between 0 and 90 .
a Deviatoric stress q-deviatoric strain e q plot, b the secant shear stiffness G-deviatoric strain e q plot, c Rendulic planes of intergranular strain h and stretching D for the state indicated by circles in (a, b) The simulations with IGS-B and ISA-B are similar. a Deviatoric stress q-deviatoric strain e q plot, b the secant shear stiffness G-deviatoric strain e q plot, c Rendulic planes of intergranular strain h and stretching D for the state indicated by circles in (a, b) [ 90 strain path reversal For rotations of the strain path larger than 90 , the behaviour follows from the elastic model until the angle between h and D is smaller than 90 , see Fig. 4. 4 Note that different to the 180 strain path reversal, the interpolation between elasticity and barodesy starts before Initial state h 0 ¼ c 0 ¼ 0 Figure 5 shows the response of an undrained triaxial test, with h 0 ¼ 0 and in addition for ISA-B c 0 ¼ 0. In order to highlight the differences between the models, v is set to 1. ISA-B is initially elastic until h reaches the yield surface (jjhjj ¼ R=2) and thus at e q ¼ ffiffiffiffiffiffiffi ffi 2=3 p Á R=2, see Fig. 5b, c. The deviatoric strain at which ISA starts the interpolation/transition between elasticity and barodesy is marked. The stiffness computed by the IGS decreases from the beginning (e q ¼ 0, h ¼ 0) due to the interpolation between barodesy and elasticity until jjhjj ¼ R.
Simulations with IGS-B and ISA-B. a Deviatoric stress q-deviatoric strain e q plot, b the secant shear stiffness G-deviatoric strain e q plot, c Rendulic planes of intergranular strain h and stretching D for the state indicated by circles in (a, b). ISA behaves purely elastic until e q ¼ ffiffiffiffiffiffiffi ffi 2=3 p Á R=2 (see Supplementary file 3: Animation related to Figure 5) 4 In order to obtain a strain path ! 90 , the following strain path before undrained compression is chosen: Figure 6 shows again the response of an undrained triaxial test again with h 0 ¼ 0 and for ISA- On the contrary, the stiffness computed by IGS results from the interpolation between barodesy and elasticity until the decrease in stiffness for small strain is also very low, so there is hardly any difference between ISA-B and IGS-B. 5

Oedometric compression
As for the undrained triaxial compression tests, the response of IGS-B and ISA-B is similar for oedometric compression. Figure 7 includes an oedometric simulation with loading, un-and reloading of an initial normally consolidated sample with the parameters according to Tables 3 and 4 Summarizing, it can be seen that with a certain parameter set, the response of ISA and IGS is very similar for low cycles or only one strain paths reversal.

Cyclic element tests
It is therefore interesting to further investigate conceptional differences between IGS-B and ISA-B with a focus on a higher number of cycles. Different to IGS, ISA includes a yield surface (small bubble) wherein the behaviour is purely elastic for a certain small-strain range. This fact becomes visible in cyclic experiments with small-strain amplitudes jjhjj\R and is explained in this section.
CU tests Figures 8 and 9 contain simulations of undrained triaxial compression with 100 applied deviatoric strain cycles De q ¼ 7:5 Á 10 À5 \ ffiffiffiffiffiffiffi ffi 2=3 p R. For a better visibility, only the first three and last three cycles are displayed. The initial A and thus the initial volumetric intergranular strain h 0;v ¼ 0 and the deviatoric part h 0;q ¼ ffiffiffiffiffiffiffi ffi 2=3 p R. In order to better show the differences in the response, we start with a small value for v ¼ v 0 ¼ v max ¼ 1 in Fig. 8.
R, h in the simulation with ISA-B stays inside the elastic bubble (Fig. 8b, c). Thus, there is a purely elastic response and no accumulation in deviatoric stress, see Fig. 8a. For IGS-B, the response follows from the elastic model after each loading reversal and is then controlled by transition as soon as h and D point in the same direction (in general for h : D [ 0). In Fig. 8c, the elastic response of IGS () is indicated by the bold line, and the thin line indicates interpolation between elasticity and barodesy, which is responsible for the deviatoric stress q accumulation in Fig. 8a.
Note that for high values of v ¼ 20, also for IGS-B no accumulation of deviatoric stress q is obtained even though Simulations with IGS-B and ISA-B. a Deviatoric stress q-deviatoric strain e q plot, b the secant shear stiffness G-deviatoric strain e q plot, c Rendulic planes of intergranular strain h and stretching D. ISA is described through elasticity until e q ¼ ffiffiffiffiffiffiffi ffi 2=3 p Á R=2. Note that the shear stiffness of IGS-B is higher than that of ISA-B the response of IGS-B is obtained from the elastic model and also from transition between elasticity and barodesy, see Fig. 9. This is also visible from Fig. 6b, where the response for IGS-B (no elastic behaviour for 0\De q \ ffiffiffiffiffiffiffi ffi 2=3 p R) is as stiff as for ISA-B (purely elastic for 0 De q \ ffiffiffiffiffiffiffi ffi 2=3 p R). In Figs. 10, 11 and 12, cyclic undrained triaxial tests with applied strain increments of De q ¼ 10 À4 are investigated applying 100 cycles. The sample is initially normally consolidated, i.e. h 0 ¼ À1= ffiffi ffi 3 p RI and c 0 ¼ h 0 =2. The deviatoric strain increment De q ¼ 10 À4 reaches into the transition zone for both ISA-B and IGS-B. In Fig. 10, v ¼ v 0 ¼ v max is set to 1. It is interesting to note that the obtained accumulation in stress is slightly larger for ISA-B. Still, the results are comparable.
For the simulations in Fig. 11, The obtained accumulation is highly reduced for both models    Fig. 10. Again, the stress accumulation is slightly higher for ISA-B.
If v develops from v 0 ¼ 1 to v max ¼ 20, the results for ISA-B according to Fig. 12 are obtained: There is a higher stress accumulation for the first cycles (compared to Fig. 10) and reduced accumulation for a higher number of cycles (compared to Fig. 11). Note that IGS-B used in this work does not allow an evolution of v. However, very similar results as obtained with ISA-B are expected with IGS according to [5].
Arbitrary paths: It is of further interest to investigate non-conventional strain paths as the star-shaped strain path displayed in Fig. 13a. We apply these star-shaped strain paths with 100 very small-strain cycles inside the elastic ISA-B bubble jjhjj\R, see Fig. 13b. The initial values are A and for ISA-B c 0 ¼ 0. In   Fig. 13, v is set to v 0 ¼ v max ¼ 1. The obtained stress response is displayed in Fig. 13c. As both ISA and IGS include a hypoelastic formulation, there is a stress accumulation 6 , also for ISA where only the elastic model affects the results. The accumulation is certainly higher for IGS as there is transition between barodesy and hypoelasticity (for h : D [ 0) and v is set to a low value of 1. However, the results between ISA-B and IGS-B become very similar when v ¼ v 0 ¼ v max is set to 20, even though IGS-B does not include a purely elastic strain range.
R performing 100 very small-strain cycles The results for ISA and IGS are comparable. a Deviatoric stress q-deviatoric strain e q , b stress accumulation in p 0 -q space

Comparison with laboratory data
In the following, both models are validated through simulations of the laboratory data presented by Wichtmann & Triantafyllidis [37] on Kaolin silt. Due to its low plasticity, the behaviour of this material is only marginally time dependent. Consequently, it has no disadvantage to perform the experiments with a time-independent stress-strain rate relation like the new models ISA-B or IGS-B [1].
The experiments comprise an oedometric test with loading, unloading-reloading cycles, undrained monotonic triaxial tests with variation of initial mean pressure p 0 0 and undrained cyclic triaxial tests with variation of the deviatoric stress amplitude q ampl , initial mean stress p 0 0 , initial overconsolidation ratio OCR 0 and initial stress ratio g 0 . Furthermore, four undrained monotonic triaxial tests with variation of initial mean pressure p 0 0 performed on Lower Rhine Clay [35] are simulated with the new model. The adopted parameters for the simulations are listed in Table 5 for the basic barodesy model and in Table 6 for ISA or IGS extension. Until otherwise specified, black dashed lines represent the laboratory data from [37] in each figure. Orange or blue solid lines represent the simulations with either ISA-B or IGS-B, respectively.
For all simulations, the intergranular strain tensor h has been initialized at fully mobilized state h 0 ¼ À1= ffiffi ffi 3 p RI. The kinematic hardening tensor of ISA has been initialized to c 0 ¼ 1=ð2 ffiffi ffi 3 p ÞRI. The internal variable for cyclic history was initialized considering only monotonic previous loading of the sample to e a;0 ¼ 0.

Oedometric test
A one-dimensional (oedometric) compression test with three unloading-reloading stress cycles is included in Fig. 15. In both, the experiment and the numerical calculations, all processes were stress-controlled. The three parameters N; k Ã and j Ã of the basic barodesy model are calibrated to the values depicted in Table 5 using this experiment.
Both models give comparable results, since the loading is mainly monotonic or occurs with large unloadingreloading loops. This also confirms the qualitative behaviour illustrated in Fig. 7. A very good agreement of the numerical calculations with the laboratory data is evident even at small stresses during loading. The compressibility, hence the change in the void ratio, is very well reproduced. Also, the hysteretic behaviour is described well by the models, yet slightly overestimated. The overestimation arises from the range of small stresses (e.g. r 0

Monotonic triaxial tests
The predictive capabilities of the constitutive model under monotonic behaviour are also checked through the simulations of undrained triaxial tests for different confining pressures. Five undrained triaxial paths of initially normally consolidated and reconstituted samples with p 0 0 ¼ f50; 100; 200; 300; 400g kPa are considered in Fig. 16.
One may check with these simulations the performance of the four material parameters of the basic barodesy model considering that no cyclic loading has been performed. Thus, both ISA-B and IGS-B are expected to provide similar paths for large strains, as is also shown through the qualitative simulations of undrained triaxial tests subsequent to 90 strain path reversal in Fig. 3 for Therefore, the degradation of shear modulus in the transition area occurs faster with IGS. In fact, the effective stress paths, see Fig. 16a, of the two models differ from each other mainly at the beginning. The paths calculated with ISA-B start from the isotropic state vertically upwards (pure shear), whereas those calculated with IGS-B start slightly inclined to the left, which confirms the qualitative behaviour presented in Fig. 3b. In both models, the stiffness and stress response is obtained through interpolation between the elastic model and barodesy, as the angle between h and D is between 0 and 90 . It can also be seen in Fig. 3 that the responses of IGS and ISA differ after a strain path reversal of 90 in sense of shear modulus degradation at transition area. Of course, both models provide the same shear strength at critical state. Apart from the previously described slight differences between the models, they also provide comparable results in the strain-stress space, see Fig. 16b.
Comparing these simulations with the experimental results, some issues arise:   (iii) the stress-strain relationship, see Fig. 16b, is in average well reproduced. However, the initial stiffness of both models overestimates the stiffness of the material. The reason for this also stems from the lack of inclusion of the inherent anisotropy into the basic barodesy model, which can be   Tables 5 and 6 Acta Geotechnica (2022) 17:4333-4358 4345 described by the transversal isotropic elasticity [11,12,31].
These issues motivated to simulate similar experiments of another clay and check the performance of the basic barodesy model for fine-grained soils. In that course, Fig. 17 comprises four undrained triaxial tests of initially normally consolidated and reconstituted Lower Rhine Clay samples with p 0 0 ¼ f50; 100; 200; 400g kPa from [35]. Hereby, the numerical calculations are performed only with the barodesy with ISA model as for monotonic loading both models provide similar results. Both the curvature of the effective stress path and the shear strength are reproduced quite satisfactory, see Fig. 17a. This is also reflected in the stress-strain relationship, see Fig. 17b, which is reproduced almost flawlessly by the model. Hence, as denoted also in [35] the Lower Rhine Clay does not possess a crossisotropic elasticity in contrast to Kaolin.

Cyclic triaxial test
In several test series in [37], the deviatoric stress amplitude q ampl , the initial stress ratio g 0 ¼ q 0 =p 0 0 as well as the initial overconsolidation ratio OCR 0 have been varied on reconstituted Kaolin samples. The cyclic loading was applied with a constant displacement rate. In order to test a certain deviatoric stress amplitude, the loading direction was changed once the specified minimum or maximum values of deviatoric stress were reached (pseudo-stress control). The cyclic loading was stopped when a certain value of axial strain (failure criterion, usually je 1 j ¼ 10%) was reached.
All simulations have been carried out with both ISA-B model and IGS-B model. The adopted parameters for the basic barodesy model are given in Table 5, and those for the small strain stiffness extensions are listed in Table 6. It is important to note that all simulations have been performed with the same parameter set, i.e. an average best fit of all experiments was targeted.
From the mathematical point of view, it is important to capture the fact that contrary to the qualitative calculations in Figs. 8 and 9, where the cyclic loading was within the yield surface of the intergranular strain, the cyclic loading in this section is in the transition region and passes into the fully mobilized region. Furthermore, the cyclic loading in the qualitative simulations presented in Figs. 10, 11 and 12 reaches into the transition area with De q ¼ R.
In Fig. 18, the effective stress paths of the experiments are compared to the numerical calculations. Except for the slope of the path, which is, as expected, not correctly reproduced, the simulations with both models fit very well to the laboratory data. It can be observed that both models show very well the trend of decreasing number of cycles to failure with increasing deviatoric amplitude, which is also evident in Fig. 20. The larger the loading amplitude, the faster the relaxation of the effective mean pressure and thus the critical state is approached faster. A state with zero effective stress is not reached in these tests on Kaolin with isotropic consolidation. However, the lower the deviatoric amplitude, the closer to the origin of the p 0 À q plane are the measured final eight-shaped effective stress loops. This behaviour is only partially reproduced by the modelsonly until the critical state is reached. From this state on, the stress path no longer moves to the left. Note that this behaviour could not be foreseen with the qualitative simulations in Sect. 3 as there the stress state was far away from the critical one.
The stress-strain relationships of these experiments are shown in Fig. 19. Due to the inherent anisotropy, the accumulation of the axial strain occurs in the compression area, where also the failure strain of e 1 ¼ 10% is reached. The models can only reproduce an accumulation of the axial strain in extension, which is also expected from the qualitative simulations in Figs. 10, 11 and 12. Despite of this, all features of the experiments are very well reproduced with the ISA-B model. The IGS-B model overestimates the strain accumulation at each experiment. For C04, it even reaches axial strains of e 1 % À100%. These shortcomings are not present in the ISA-B model due to the ISA internal variable for cyclic history 7 , which was not expected from the qualitative simulations in Sect. 3 as the respective simulations ranged either in the elastic area or in the transition area and the parameters controlling the accumulation rate where simplified to , the accumulation of ISA-B slightly exceeded the one of IGS-B. Only in the simulation of Fig. 12, where the evolution of the parameter v was not turned off (v 0 ¼ 1; v max ¼ 20), the accumulation of the ISA-B model was reduced for a higher number of cycles, suggesting the performance of the model observed in this section. Figure 20b presents the axial strain amplitudes against the number of cycles. Even though ISA-B fits better the experimental behaviour, it also underestimates the axial Fig. 18 Simulation of cyclic triaxial tests performed on Kaolin with variation of deviatoric amplitude q ampl ¼ f30; 45; 60; 70g kPa. Parameters according to Tables 5 and 6 Acta Geotechnica (2022) 17:4333-4358 4347 Fig. 19 Simulation of cyclic triaxial tests performed on Kaolin with variation of deviatoric amplitude q ampl ¼ f30; 45; 60; 70g kPa. Parameters according to Tables 5 and 6 strain amplitudes for the lower deviatoric stress amplitudes. The trend of faster increasing strain amplitude with increasing stress amplitude is well described by both models. The curves of accumulated pore water pressure u acc ðNÞ ¼ uðNÞ À uðN ¼ 0Þ are provided in Fig. 20a. A value of u acc ¼ 200 kPa, which means zero effective stress, has neither been reached in the present tests on Kaolin nor in the numerical calculations. As expected, the accumulation of pore water pressure and thus the relaxation of mean effective stress evolves faster with increasing stress amplitude, which is very well reproduced by the models. Some discrepancies for the last cycles between the simulations with IGS and the experiments for the higher deviatoric stress amplitudes can be observed.
Variation of initial stress ratio g 0 at p 0 0 ¼ 200 kPa Another four experiments are presented in Figs. 21 and 22. Therefore, the initial stress ratio has been varied between g 0 ¼ q=p 0 ¼ f0:25; 0:125; À0:125; À0:25g (tests named chronologically C26-C29) for p 0 0 ¼ 200 kPa. All samples were initially normally consolidated and have been subjected to a cyclic stress amplitude of q ampl ¼ 30 kPa with the displacement rate of _ s ¼ 0:1 mm/min. The considered number of cycles for each test amounted to N ¼ 150.
The effective stress paths are depicted in Fig. 21. Similar to the previous observations for the tests with variation of stress amplitude, the simulations with both models are in a very good agreement with the experiments apart from the slope of the path. Again, a state with zero effective mean pressure is not attained either in compression (tests C26 and C27) or in extension (tests C28 and C29) regime, which is well reproduced by the models.
The stress-strain relationships of the experiments and simulations are presented in Fig. 22. Therefore, a considerable accumulation of permanent axial strain was observed, while the axial strain amplitude remained almost constant in contrast to the previous described tests with isotropic consolidation. The accumulation of strain continued even after the pore pressure accumulation had stopped, hence, even though the effective stress path no longer moved to the left. In contrast to the tests with isotropic consolidation, the failure criterion of je 1 j ¼ 10% was reached due to excessive permanent axial strains, and not due to too large strain amplitudes, which is very well reproduced by the models. Furthermore, the accumulation  Tables 5 and 6 Acta Geotechnica (2022) 17:4333-4358 4351 direction is characterized by the direction of the preceding anisotropic consolidation. Hence, the samples consolidated anisotropically in triaxial compression (g 0 [ 0) accumulate permanent axial strains in compression regime and vice versa. This behaviour is also very well simulated by the constitutive models. It is hereby recalled that the isotropically consolidated samples showed a pronounced accumulation of axial strain in compression direction due to the strong inherent anisotropy, see Fig. 19. Furthermore, especially for the samples with anisotropic consolidation in compression regime, a strong overestimation of the accumulated axial strain can be observed by the simulations with IGS-B. Of course, at the expense of other laboratory tests, one could reproduce these experiments with IGS-B more accurately using a different set of parameters. The parameter set shown in this work has been selected as the mean best fit for all experiments. Nevertheless, the observed excessive strain accumulation presented in Fig. 19 was expected for medium values of the parameter v considering the examined qualitative behaviour of IGS-B in Sect. 3, Figs. 13 and 14. Higher values of v would render in overall a lower strain accumulation, which may be accurate only for experiment C26 or C27. Hence, a single parameter set which would capture accurately different deviatoric stress amplitudes as well as different initial stress ratios is not possible.
The effective stress paths obtained from these three tests as well as the simulations with the two considered models are provided in Fig. 23. The dilative response observed in the experiment during the first cycle is reproduced by the models only for OCR 0 ¼ 2:5 (test C39), see also the negative excess pore water pressure in Fig. 25a. The decreasing pore pressure accumulation with increasing OCR 0 , resulting in a decay of mean pressure relaxation, is very well captured by the models. Also, the decreasing rate of pore pressure accumulation with increasing OCR 0 is accurately described by the numerical calculations, see Fig. 25a. Figure 24 presents the stress-strain relationships, whereby the failure criterion of je 1 j ¼ 10% has been reached only for the initially normally consolidated sample (test C37). In all experiments, IGS-B overestimates the excessive permanent axial strains, whereas ISA-B shows a good agreement with the experiments. The accumulation direction is observed in the experiments to occur in extension regime, which is well captured by both models.
The axial strain amplitudes e ampl 1 versus the number of cycles are illustrated in Fig. 25b. It can be observed that the number of cycles to failure grows considerably with increasing OCR 0 in both experiments and numerical calculations. Furthermore, IGS-B overestimates the axial strain amplitude especially for OCR 0 ¼ 1.

Summary, conclusion and outlook
In this work, barodesy [22] is extended with ISA plasticity [27] to improve the small-strain predictions of the basic model. Comparisons are made with IGS [26] as well as with experiments [37]. The main differences between ISA and IGS comprise the following aspects: (i) only ISA includes a purely elastic strain range, (ii) the ISA version used in this article allows a development from v 0 to v max using the so-called internal variable for cyclic history 0 e a 1 .
As a result of (i), it is obvious that purely elastic behaviour for very small-strain cycles is only obtained for ISA. However, for high values of v, the results of ISA and IGS are very similar, even for very small-strain cyclic behaviour. As the elastic formulation is hypoelastic, an accumulation of stress/strain is still obtained for arbitrary stress/ strain paths inside the elastic ISA bubble. The predictions of the experiments performed on Kaolin for a higher number of cycles reveal a few advantages for the ISA model, which results from the fact (ii), hence the inclusion of a transition state from v 0 to v max . Equivalent results to ISA are also expected for IGS according to [5]. In fact, if the model allows v 0 to develop to v max and introduces the state variable responsible for cyclic mobility, then also for very small strain or stress cycles, the stress or strain accumulation vanishes even without a predefined purely elastic bubble.
Future work will focus on the extension of the new coupled model for inherent anisotropy, which is shown in the quantitative simulations to present a drawback of the constitutive model. When dealing with medium to highly plastic clays, an influence of the loading frequency on the accumulation behaviour as well as relaxation and creep effects is existent. Therefore, the model will be extended to account for rate-and time-dependent effects of fine-grained soils.

Symbols and notations
In this article, the symbolic notation is used for the effective Cauchy stress T and stretching D, but in some cases the more familiar symbol r 0 i instead of T i is used for the principal stresses. For the principal components of stress, compression is defined negative. Tensors are written in bold capital letters (e.g. X). jjXjj ¼ ffiffiffiffiffiffiffiffiffi ffi tr X 2 p is the Euclidean norm of X, and tr X is the sum of the diagonal components Fig. 23 Simulation of cyclic triaxial tests performed on Kaolin with variation of initial overconsolidation ratio OCR 0 ¼ f1:5; 2:0; 2:5g. Parameters according to Tables 5 and 6 Acta Geotechnica (2022) 17:4333-4358 4353 of X. The superscript 0 marks a normalised tensor, i.e. X 0 ¼ X=jjXjj. Stresses are considered as effective ones. T is the co-rotational, objective stress rate. The stretching tensor D is the symmetric part of the velocity gradient.
The void ratio e is the ratio of the volume of the voids V p to the volume of the solids V s . Note that with p 0 ¼ À 1 3 tr T, the mean effective stress is positive for compression. e v ¼ tr e is the volumetric strain. For compressive strain, e i is defined negative. Fig. 24 Simulation of cyclic triaxial tests performed on Kaolin with variation of initial overconsolidation ratio OCR 0 ¼ f1:5; 2:0; 2:5g. Parameters according to Tables 5 and 6 For axisymmetric conditions, often the Rendulic plane is used. For a conventional triaxial compression or oedometric compression test, the axial stress is denoted with r 0 1 and the radial stress is denoted with r 0 2 ð¼ r 0 3 Þ. The associated strains are e 1 and e 2 ¼ e 3 . Other symbols used are q ¼ Àðr 1 À r 2 Þ and e q ¼ À2=3 Á ðe 1 À e 2 Þ. Initial values are labelled with the subscript 0.
Any symmetric second order tensor can be written as vector with the principal values X v ¼ X 1 ; X 2 ; X 3 ½ . We use this to display tensors in figures; however, we do not use the notation X v , as it is implicitly clear that X is shown as a vector in these figures. Bold calligraphic letters denote tensors of 4th order (e.g. M). We use different kinds of tensor operations employing the Einstein summation convention. In particular, the indices follow the lexicographic order: X Y ¼ X ij Y kl , X : Y ¼ X ij Y ij and L : D ¼ L ijkl D kl . We employ the unit tensor of second order I with I ij ¼ d ij and the fourth order tensor I with I ijkl ¼ d ik d jl , using the Kronecker delta d ij .

Equations of barodesy
In this appendix, all equations of barodesy for clay [22] are summarized.
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/.