Small Strains in Soil Constitutive Modeling

This paper reviews the state of the art of soil behavior in the range of small strains and its constitutive modeling, which is an important issue when predicting displacements under serviceability conditions. The factors that control nonlinear, hysteretic and dependent on recent history soil behavior are described. Likewise, concepts of soil constitutive modeling are explored in detail and two criteria are explained and used to classify the analyzed models: (1) a first criterion based on the concept of tensorial zones; and (2) a second criterion based on the elements that defines the hysteretic soil behavior, including reversal criteria, memory rules and the effects of reversals on soil degradation and on soil stiffness recovery. The fundamentals of the formulation of the analyzed models are provided, as well as their scope of application, advantages and disadvantages.


Introduction: The Kinematic Nature of Soil Stiffness
The stiffness of a body, understood as its resistance to deformation under applied forces, depends on its shape, its boundary conditions and the stiffness of its constitutive materials [1].
From the experiments carried out by Hardin and Drnevich [2], Simpson et al. [3], Jardine [4] and Smith et al. [5], Jardine [6] developed the idea previously stated by Skinner [7] and supported by a large number of researchers [8][9][10][11][12][13], to explain soil stiffness behavior depending on the range of stress and strain to which it is subjected, within the framework of the theories of Kinematic Yield Surfaces (KYS) and plasticity. Jardine [6] differentiated four zones around any point in stress space, which move and change their shape according to the stress paths followed. In each of these zones, Jardine identified a characteristic behavior of the soil, which was related to the nature of the strains that take place in them (Fig. 1).

Zone I
Externally limited by Y 1 (Fig. 1a), soil behavior in this area is linear reversible [6,[14][15][16][17][18][19][20][21]. Some researchers indicate that a purely linear elastic behavior does not occur in any strain range, although it is a sufficient approximation from a practical point of view [1]. Likewise, Hueckel and Nova [22] affirm that irreversible deformations occur in all primary loading processes. Zone I of Jardine can be very small [23][24][25], especially in uncemented and normally consolidated soils [26,27]. Zone II Internally limited by Y 1 and externally by Y 2 (Fig. 1a), soil behavior in this zone is nonlinear reversible, hysteretic and dependent on recent history [3,4,6,15,16,20,22,[28][29][30][31][32][33][34]. Zone III Internally limited by Y 2 and externally by Y 3 (Fig. 1a), soil behavior in this zone shows a first degree of irreversibility (plastic deformations), which depends on the OCR [6,35]. As can be seen in Fig. 1b, if OCR [ 1, plastic strains in this zone tend to increase linearly with the strain, while if OCR ¼ 1, they tend to do it in a nonlinear way. Zone IV Internally limited by Y 3 (Fig. 1a), soil behavior in this zone shows a high degree of irreversibility [6,35]. The plastic strains in this zone always increase in a nonlinear way with the strain (Fig. 1b).
The definition of these zones is useful to understand and to classify the soil response in the context of mechanical constitutive equations, and they will be used in this paper for this purpose.

The Range of Small Strains in Engineering Practice
The application of Soil Mechanics to the new problems faced by civil engineers at the beginning of the twentieth century was focused on the prevention of collapse, but this conception changed in the 70 s due, among other reasons, to: (1) a good knowledge of most soil failure mechanisms; (2) the need to build new structures in dense urban areas and to protect nearby existing structures; (3) the need to build especially sensitive structures, such as nuclear plants; and (4) the advances in numerical analysis tools and in computing power [1]. During the 70 s, it was observed how the soil stiffness obtained in conventional laboratory tests did not coincide with the stiffness obtained in dynamic field tests [8].
The dynamic field tests stiffness was usually an order of magnitude higher than the one of laboratory tests [36]. Later, researchers [15,37,38] realized that these differences were due to stiffness measurement at different strain levels, which allowed unifying the theories that tried to explain, until that moment, the value of the soil shear modulus. These facts suggested that conventional laboratory tests were not adequate to assess natural soil stiffness. In addition to that, it was shown that the low stiffness measured in conventional laboratory tests also came, among others, from deficiencies in the strain measurement techniques and from the disturbance of the samples during their extraction process [9]. During the subsequent years, measurements of soil displacements in various London works, new instrumentation techniques in laboratory tests and advances in numerical analysis methods, helped to improve the understanding of soil behavior in the range of small strains, in which stiffness plays a fundamental role [15,36]. These findings promoted the development of new instrumentation techniques for field and laboratory tests, which Clayton [1], Elhakim [39] and Cudny [40] classify as follows: • Dynamic field tests: Cross-hole [41], Down-hole [42][43][44], Spectral Analysis of Surface Wave SASW [45,46], Continuous Surface Wave System CSWS [47], Suspension logger [48], Seismic cone CPTU [49] and Seismic dilatometer SDMT [50]. • Laboratory tests: Resonant column [2], Hollow cylinder [18,51], triaxial devices with internal strain measurements [9,20] or bender elements and compression and shear seismic wave measuring devices [52,53] in oedometers [18,[54][55][56][57], in direct simple shear apparatus [55], in triaxial apparatus [55,58,59], in resonant column apparatus [60][61][62], in double simple shear apparatus [63] and in unconfined samples immediately after their extraction [64,65].
There are numerous advanced constitutive models capable of simulating soil behavior in the whole range of strains (Zones I, II, III and IV of Jardine [6]), which is essential for a complete analysis of geotechnical problems, although the use of many of these models is reduced to an academic ambit and according to Tamagnini and Viggiani [66] this is due to: (1) its complex formalism; (2) the use of a large number of parameters which are difficult to obtain experimentally; (3) the use of a large number of state variables which are difficult to initialize; and (4) the difficulties associated with the formulation of precise, efficient and robust algorithms for the numerical implementation of equations (although there have been significant advances in recent years regarding this point). However, some of these advanced models have managed to extend to the professional practice, for example, the Hardening Soil Small model (HS-S) [67]. Axial plastic strains e p À Á normalized by the maximum axial total strain (e max ) vs. total shear strain (e s;max ). From [6] Figures 2, 3 and 4 show the need to consider soil behavior in the range of small strains in different geotechnical problems. It can be observed that models incorporating small strain stiffness improve the prediction of soil displacements with respect to the conventional elastoplastic models. and J4, from [68]. b Numerical simulations with a linear elastic model with perfect plasticity and two nonlinear models, Brick and SRD-Brick, from [69] Fig. 3 Soil displacements during an excavation between diaphragm walls in the Taipei Enterprise Center (TNEC), from [70]. a Numerical simulations with the MCC model. b Numerical simulations with the USC model 3 Soil Behavior in the Range of Small Strains

The Range of Small Strains
Jardine [6] identifies soil behavior in the range of small strains with Zones I and II. Zone I usually corresponds to strains less than 10 À6 in sands and 10 À5 in clays, while Zone II does so with strains within the interval 10 À6 to 10 À3 [40]. On the other hand, the range of intermediate and large strains is identified with Zones III and IV defined by Jardine. These zones usually correspond to strains greater than 10 À3 .

Bulk Modulus
Regarding the drained bulk modulus, it is not usual to distinguish between values at small, intermediate or large strains, unlike what happens with the shear stiffness modulus. Duncan et al. [71] proposed the expression (1) to calculate the apparent bulk secant modulus (Fig. 5).
The use of the confinement r 0 3 as a variable in the expressions of the stiffness moduli is a common practice in many hypoelastic and quasi-hypoelastic models. This is because soil parameters of many models are usually adjusted to theoretical curves based on the results obtained in the deviatoric phase of triaxial tests, in which the confinement remains constant, which simplifies the fitting. However, the use of the r 0 3 as a variable in these expressions makes it difficult to generalize these models to multiaxial loading states. The expression (2) for K 0ap s is considered more appropriate.
Duncan et al. [71] proposed an expression analogous to (2) for the drained tangent elastic longitudinal modulus E 0 t;ur p 0 ð Þ, from which it is possible to calculate the drained tangent elastic bulk modulus K 0 t;ur ¼ E 0 t;ur = 3 1 À 2m 0 ur À Á À Á . Likewise, Roscoe and Schofield [73] and Roscoe and Burland [74] proposed the following expression for the elastic tangent bulk modulus in their critical state models (Fig. 6).
Lade and Abelev [76] analyze the hysteretic volumetric soil behavior through the study of the variation of volumetric stiffness during isotropic loading and unloading in sands due to the introduction of small loading cycles (Fig. 7). In the primary loading branch, they observed that the cycles led to a significant increase in volumetric stiffness, which did not happen during the cycles in the unloading branch, in which any differences in the value of such stiffness were hardly observed.

Shear Stiffness Modulus
Unlike the bulk modulus, shear stiffness modulus in soils is more difficult to specify and there are many studies about the parameters controlling it. The parameters considered most relevant are: shear strain (c), mean stress (p 0 ), void ratio (e), plasticity index (PI), overconsolidation ratio (OCR, R 0 ), diagenesis, recent history, loading rate and anisotropy. Fig. 5 Dependence of the apparent secant bulk modulus (B K 0ap s ) on r 0 3 =p a in high plasticity clays from Barind, according to the expression (1), from [72] Fig. 6 Tests in silty clays, from [75]. Isotropic compression, t À 1 ¼ 1 þ e 0 ð Þ e v : • Shear strain (c): There are many evidences indicating that shear stiffness modulus G degrades with shear strain. To generalize the concept of shear strain to a multiaxial state of stress and strain, it is usual to work with the octahedral shear strain. Figure 8 shows the results of diverse experimental tests involving the degradation of the apparent shear modulus with the total shear strain in clays and sands. • Mean stress (p 0 ): Shear modulus depends on confinement. Authors such as Ohde [79], Hardin [80], Janbu [81], Hardin and Richart Jr. [82] or Hardin and Drnevich [2], based on experimental observations, proposed relations as G 0 / p 0 ð Þ m between the value of maximum shear modulus G 0 and mean stress (Fig. 9). The introduction of the Hertzian contact theory in spherical particles for the calculation of G 0 results in m ¼ 0:33 [86,87]. Experiments in sands provide values of m ¼ 0:40 À 0:60 [88][89][90], while in clayey soils, values of m ¼ 0:50 À 1:00 are usually taken, being close to 0:50 for low plasticity clays and to 1:00 for high plasticity clays.
• Void ratio (e): Hardin and Richart [82] observed the dependence of the maximum shear modulus G 0 on the void ratio from experiments with Ottawa sands, and proposed expressions like G 0 /B À e À Á 2 = 1 þ e ð Þ. Authors such as Biarez and Hicher [91] or Lo Presti and Jamiolkowski [92] proposed expressions as G 0 / e Àx , and Bui [93] proposed expressions as G 0 / 1= 1 þ e ð Þ 3 (Fig. 10). • Plasticity Index (PI): It is experimentally observed that higher values of the plasticity index in a soil result in a shift of the degradation curve of the apparent shear modulus towards higher shear strain values (Fig. 11).
Vucetic and Dobry [94] proposed patterns for this dependence. • Overconsolidation ratio (OCR, R 0 ): Hardin and Black [86] proposed correlations like G 0 / OCR ð Þk, based on experimental observations between the maximum shear modulus G 0 and the overconsolidation ratio OCR. On the other hand, Houlsby and Wroth [95], based on the works of Hardin and Black [86] and Atkinson and Little [96], proposed correlations as G 0 / R 0 ð Þk between the maximum shear modulus G 0 and the overconsolidation ratio R o (Fig. 12).
• Diagenesis: The diagenesis is the physicochemical process by which a sediment is transformed into a sedimentary rock. This process gradually alters the stiffness of the soil. Some of the main processes that alter soil stiffness are cementation [98] and aging [99][100][101][102][103][104], understood as the alteration of mechanical properties of the soil resulting from secondary compression under a constant external loading. Trhlíková et al. [105] proposed a relation between the maximum shear stiffness modulus G 0 and the structure as Fig. 7 Stiffening of the bulk modulus in sands of Nevada during small loading and unloading cycles in the primary loading branch and during the unloading branch [76], extracted from [67] Fig. 8 a Degradation of the apparent shear modulus in clay tests with different plasticity index, from [77]. b Degradation of the apparent shear modulus in sand tests, from [78] G 0 / s Ã =s Ã f l . On the other hand, Anderson and Stokoe [106] proposed a relation between G 0 and the aging like G 0 t ð Þ / G 0 t p À Á 1 þ N G;1 log t=t p À Á À Á . • Recent history: Atkinson et al. [30] defined the concept of recent history as that corresponding to recent stress or strain path in relation with the previous one, from which it is differentiated by a change in its direction (reversal) or by an extended period of rest. Soil behavior after a reversal and before a subsequent monotonous strain suggests a gradual adaptation of its internal state until it exclusively depends on r 0 and e [107]. In such state, that it is inside the SOM region [108], proportional strain paths lead to proportional stress paths. It is for the latter that the influence of the internal state of the soil over its behavior can only be revealed during small strains after reversals [107]. The studies of soil recent history require the use of tests with different stress or strain paths. Sayao [109] show these tests as indicated below: Fig. 9 a Influence of p 0 on the value of G 0 [83], extracted from [84]. b Influence of p 0 in the degradation of G, from [85] Fig. 10 Results of experiments with resonant column on well graded sands and reconstituted clays, from [93] • Tests without rotation of the principal stresses ( _ a r ¼ 0Þ.
Volumetric strain reflects the effect of isotropic loadings that tend to increase the value of the contact forces between particles, while deviatoric strain, usually generated by a deviatoric loading, modifies the direction of such forces and affects the stiffness of the soil. Some notable works in which soil recent history has been studied in relation to deviatoric loadings and deviatoric strains are those carried out in triaxial tests by Richardson [29], in hollow cylinder tests by Sayao [109], in biaxial tests by Topolnicki et al. [110] and in true triaxial tests by Sture et al. [111], Fig. 13. Nevertheless, there are few tests in which, in addition, suitable measurement methods have been implemented for the range of small strains, among which stand out the triaxial tests with local strain measurement [29] (despite the limitations they present in terms of the possible stress and strain paths), as well as the hollow cylinder tests [112,113].
The work of Atkinson et al. [30] shows the results of a set of triaxial tests on reconstituted overconsolidated samples (OCR ¼ 2) with London clay (Fig. 14). All of them were initially brought to the same stress state (O in Fig. 14a) using different stress paths (PO and QO in Fig. 14a). After a rest period of 3 h, they were subjected to a deviatoric phase of Dq ¼ 90kPa with Dp 0 ¼ 0 (OA in Fig. 14a). As it can be seen in Fig. 14b, the stiffness of the soil during the OA path depends on the angle of such path with respect to the previous paths (PQ and QO). The wider the angle between paths, the higher is the value of the shear stiffness G 0 at the beginning of the new path.
Clayton and Heymann [20] conducted tests on natural samples of London clay, applying the stress paths depicted in Fig. 15a. After completing the AB phase, they allowed the samples to rest for a period of 6-12 days before starting the triaxial extension (BE) and compression (BC) phases, thus allowing the soil yielding, which did not happen in the tests of Atkinson et al. [30], in which only a rest period of 3 h was left. Thus, Clayton and Heymann observed how plastic strains significantly reduced the effect of soil recent history on the G 0 value (Richardson [29] had observed a similar behavior in his experiments). Nevertheless, despite this attenuation of the path rotation effect on the G 0 value, Fig. 15b shows how, even allowing soil yielding, the stress history still has some effect on the shape of the degradation curve.  Gasparre [75] studied the effects of plastic deformations on the value of soil stiffness (through the control of resting times of the sample before applying new loads), as well as the effects of the loading magnitude prior to the rotation of stress paths, extending the research of Atkinson et al. [30] and Clayton and Heymann [20]. Gasparre et al. [34] conducted a set of tests on natural samples of London clay (Fig. 16), in which the following was observed: 1. When the stress state in recent history remains within the contour Y 2 (Zone I or II according to Jardine), -If plastic deformations are allowed, the effect of stress path rotation on the shear stiffness of the soil is reduced, Fig. 16a (as observed by Clayton and Heymann [20]). Fig. 13 a Triaxial test: stress paths in q=p 0 À e s space with different angles tan h q=p 0 À Á ¼ Dq=Dp 0 in reconstituted London clay, with OCR ¼ 2 and p 0 ¼ 200 kPa, from [29]. b Hollow cylinder test: stress and strain paths in R ¼ r 0 1 =r 0 2 À c max space, where a ¼ a r is the rotation angle of the principal stress, in Ottawa sand with p 0 ¼ 300 kPa and D R ¼ 36%, from [109]. c Biaxial test: stress and strain paths in the respective deviatoric planes in reconstituted Karlsruhe clay, from [110]. d True triaxial test: strain path in the deviatoric plane and deviatoric stress ffiffiffiffiffiffiffiffi ffi space in a Leighton Buzzard sand with D R ¼ 72%, from [111] -If plastic strains are not allowed, the same behavior observed by Atkinson et al. [30] is obtained, that is, a clear dependence of the G 0 value on the stress path rotation, Fig. 16b rotation, regardless of whether plastic strains are allowed or not before starting the stress path after the rotation conducted, Fig. 16c.
• Strain rate and inertia effects: There are numerous experimental tests on plastic soils that show the dependence of their stiffness on the strain rate [36], an effect attributed to their viscosity and plasticity. In sands, generally, this effect is very small or non-existent [114]. On the other hand, this effect is negligible in the range of small strains c\0:001% [115], and increases its relevance for larger strains 0:01%\c\0:1% [116], as it can be seen in Fig. 17. Yong and Japp [117] defined the strain rate shear modulus parameter (Fig. 18) as a G ¼ DG=D log _ c ð Þ ð Þ.
• Anisotropy: Multiple correlations have been proposed to calculate the soil maximum shear modulus value, several of which can be found in the work of Obrzud and Truty [119]. Below there are some expressions for the calculation of G 0 ij ð Þ in a plane i according to direction j, based on the relations discussed above, which also introduce the effect of soil anisotropy: Hardin and Black [89]: Hardin and Blandford [14]: Rampello et al. [120]: Pennington [121]: Pennington [121]:

Considerations about the Hysteretic Behavior
A fundamental aspect of soil behavior in Zone II of Jardine [6], along with the nonlinearity and dependence on recent history, is the hysteresis (Fig. 19).

Fig. 16
Degradation curves in undrained triaxial tests in London clay (G eq;tan G ap t ¼ _ q=3 _ e q ), from [34]. a Within Y 2 (Zones I or II), allowing plastic strains. b Within Y 2 (Zones I or II) without allowing plastic strains. c Overcoming and displacing Y 2 (Zones III or IV), allowing or not plastic strains Fig. 17 The effect of strain rate on stiffness during the deviatoric phase of an undrained triaxial test on an intact sample of London clay, from [116] Soils are formed by particles and, therefore, will experience energy dissipation during loading cycles, which will result in a hysteretic behavior [124]. Nevertheless, the micromechanical mechanism through which the soil dissipates energy during such cycles is not clear and Cudny [40] points out to two possible explanations: (1) the dissipated energy is the result of a process of local yield and friction in the contacts between particles, which are subjected to normal and shear forces, in which case, the energy absorbed by the soil would be a function of the deformation amplitude [125]; and (2) the dissipated energy is the result of a viscous behavior due to the presence of fluid in the soil pores [126].
As Hueckel and Nova [22] indicate in their fundamental work on hysteretic behavior in soils, hysteretic cycles are characterized by having different apparent stiffness after rotations in stress or strain paths. In addition to that, if a perfect hysteretic behavior is considered, the deformations would only be recoverable if the cycles starting from the same reversal point were closed. Nevertheless, according to Hueckel and Nova, there are no perfect hysteretic cycles, generally observing: (1) certain permanent strain after the closing of the cycles, which depends on the number of these; (2) viscous effects or cyclic hardening/softening; and (3) dependence of the stiffness in unloadings/reloadings on the magnitude of the accumulated plastic strain (elastoplastic coupling). Strain rate shear modulus parameter [118], extracted from [67]. a Dependence on strain rate and strain amplitude. b Dependence on the PI value Fig. 19 a Loading cycles in triaxial tests in a loose Lubiatowo sand, from [122]. b Loading cycles in undrained triaxial tests in a Speswhite kaolinite clay with p 0 ¼ 300 kPa, from [123] Small Strains in Soil Constitutive Modeling 3233 As Gudehus [107] points out, when reproducing the hysteretic soil behavior, the state of a representative soil element is not sufficiently characterized by the void ratio e and the effective stress tensor r 0 . In these cases, it is necessary to define hidden state variables v hist that cannot be observed macroscopically and represent the spatial fluctuation of the chains of internal forces between soil particles (called force-roughness effect). It is also possible to attribute to these variables the nonlinearity and the soil behavior dependence on its recent history. As previously commented, soil behavior after a reversal and before a subsequent monotonous strain path suggests a gradual adaptation of its internal state (characterized by the mentioned hidden state variables) until it depends exclusively on the stress tensor and the void ratio, state in which proportional strain paths lead to proportional stress paths. It is for this last reason that the influence of the internal state of the soil on its behavior can only be revealed during small strains after reversals. Furthermore, the state variables v hist , within the set of state variables v, are, in general, of two types: a (back stress) normally in elastoplastic or viscoelastoplastic models (these models tend to underestimate the hysteretic effects of the soil) and d (internal strain) normally in hypoplastic or viscohypoplastic models (these models tend to overestimate the hysteretic effects of the soil) [107]. Although a and d can be formally treated as strains and stresses respectively, they cannot be interpreted physically as such [107], which is due to the fact that the internal variables a and d are obtained, respectively, from r 0 and e (a is obtained from stress through the establishment of an elastic center a ¼ r 0 when a reversal takes place, and _ d is proportional to _ e and is calculated from it).

Constitutive Modeling
As proposed by Popper [127], theoretical models can capture part of reality if they are logically consistent and the hypotheses they use are not refuted by observations. Nowadays, there are several theoretical frameworks that try to explain distinct aspects of soil behavior. Within these theoretical frameworks, several constitutive soil models that consider the mechanical behavior of the soil in the range of small strains have been developed over the last decades.
The mechanical soil behavior must be explained by the interaction of the distinct phases that constitute the soil along with the internal and external actions on it. Despite the discontinuous nature of soil, most of the models used in the professional practice consider soil as a continuous medium. Although there has been important progress in models that considers soil as a discontinuous medium [128][129][130][131], the scope of this paper has been limited to continuous ones.
According to Tobita [132], any type of soil micromechanical behavior (slide, rotation, deformation or breakage of the particles or aggregates) leads to a nonlinear incremental type of macromechanical behavior, therefore, any constitutive model that considers soil as a continuous medium should be able to reproduce this behavior.
Mechanical constitutive modeling of soils can be framed in Continuum Mechanics problems, which involve a set of well-known general conservation and balance equations and a particular constitutive equation, which for simple materials [133], is expressed in (9).
The consideration of a functional F , and not a function, is due to the irreversible soil behavior that should be reproduced. In this type of behavior, knowing the state of deformation e t ð Þ at a given instant t does not allow to know the stress state at that instant and vice versa [134]. Furthermore, not any functional F will represent a valid constitutive relation [133].
As Owen and Williams [135] show, in the case of nonviscous type materials, which are the kind of materials analyzed in this paper, it is common to resort to incremental type formulations in which the functional F is reduced to a tensorial function of the type Considering that the tensorial function G is homogeneous grade one in the term of _ e, allows expressing the constitutive equation as _ r 0 ¼ E 0 t e; r 0 ; v; g ð Þ: _ e. The dependence of the tensor E 0 t on g ¼ _ e= _ e k k leads to the tensorial zone concept defined by Darve [136,137] and Darve and Labanieh [138], which will be very useful to classify the analyzed models. A tensorial zone Z is defined as the part of strain increments space in which the tensorial function G is linear with _ e. This will imply that in a certain tensorial zone Z, the tangent stiffness tensor will be independent from g, and the relation between _ r 0 and _ e will be incrementally linear, that is, t e; r 0 ; v ð Þ,8g 2 Z. Different tensorial zones conform different adjacent hypercones with the vertex in common, and the condition of continuity does not allow the election of any tangent stiffness tensors in two adjacent tensorial zones.

Some Previous Considerations About Constitutive Modeling in the Range of Small Strain
To approximate soil behavior in Zone II, it is necessary to use nonlinear hysteretic models that consider the effect of recent history on stiffness. A common way to simulate this nonlinearity is by using nonlinear elastic stiffness moduli. Such moduli depend exclusively on the elastic strain e e or, equivalently, on the effective stress r 0 . However, the experimental measurement of soil stiffness is normally conducted in tests in which soil is subjected to primary loading processes. During such processes, it is usual to obtain stiffness values lower than those measured during unloading or reloading processes. According to Hueckel and Nova [22], this is due to the fact that irreversible strains occur in every process of primary loading. Based on this, it is possible to define the concept of apparent stiffness moduli of the soil, which are calculated using the total strain and not the elastic ones. The apparent stiffness moduli usually depend on the total strain e, the effective stress r 0 and the state variables v hist .
From a numerical point of view, using the apparent stiffness moduli instead of the elastic ones within an elastoplastic model significantly simplifies the calculation algorithms ( Fig. 20), although it can lead to theoretical inconsistencies.
Tensorial zones [136][137][138] can be bijectively related to the apparent stiffness of the soil and, therefore, to the aforementioned apparent stiffness moduli. In this way, for each tensor zone Z i a value of K 0ap i and G ap i will be obtained. Table 1 summarizes these correspondences.
In any model that uses elastic theory, it should be considered what is pointed out by Zytynski et al. [139] regarding to the choice of elastic moduli. They showed that considering constant values of drained Poisson's ratio in a model implies that this model will not be conservative and the generation of energy in closed loading cycles will take place. In fact, several incrementally multilinear models that consider soil behavior in the range of small strains, such as those of Papadimitriou et al. [140], Wongsaroj [141], SSOM and HS-S [67] or HS-SS of Plaxis, use as elastic parameters the shear modulus G and a constant value of the drained Poisson's ratio. Based on these parameters, such models internally calculate the value of drained bulk modulus K 0 using expression (12). Other incremental multilinear models that consider soil behavior in the range of small strains, such as those of Whittle [142], Al-Tabbaa and Wood [10], Yu [143], Gryczmanski et al. [144] or Gryczmanski and Uliniarz [145], use the drained bulk modulus K 0 and a constant value of the drained Poisson's ratio as elastic parameters. Based on these parameters, such models internally calculate the value of the shear modulus G through expression (12).
From the previous expressions it follows that the fact of adopting a constant value of drained Poisson's ratio implies the proportionality K 0 / G, which generally does not respond to experimental observations in the range of small strains (Fig. 21). The expression of Poisson's ratio as a function of K 0 and G is given by (13), as well as the elastic thermodynamic limitations and the condition m 0 [ 0 that must be satisfied in all soils, based on experimental observations. On the other hand, in the expression (14), _ m is deduced from expression (13).
Experimentally, it is observed how, in the range of small strains, the drained bulk modulus K 0 does tend to stiffen with the volumetric deformation, while the shear modulus G tends to degrade with the octahedral shear strain.
Incrementally nonlinear models #Z i ! 1 i ! 1 Fig. 21 (a) Variation of m 0 with c oct in cohesive soils, from [146]. (b) Variation of m 0 with c oct in sandy soils, from [146]. (c) Variation of m 0 with e oct in different sandy soils, from [147] Considering on the one hand G ¼ G c oct ð Þ with _ G _ c oct \0 (degradation of G with c oct ), and on the other

Classification of Models Using the Elements that Define the Hysteretic Behaviour
To reproduce the hysteretic soil behavior in a constitutive model it is necessary to define hidden state variables v hist [107], within the set of state variables v.
In addition to that, it is necessary to distinguish and define the following concepts in each model, which are directly related to the state variables v hist and define a specific classification criterion for constitutive models: (1) reversal criteria; (2) memory rules; (3) effect of reversals on the variables that control degradation; and (4) effect of reversals on maximum soil stiffness.

Reversal Criteria
The models that consider the hysteretic behavior of the soil use reversal criteria to identify points where changes of direction occur in the stress or strain paths, whose effect induces changes in soil stiffness. These criteria can be divided as follows: Extrinsic Reversal Criteria One or more loading and unloading criteria are defined and added to the model equations. These criteria are usually formulated based on stress [143,144], strain [67,140,142,148,149] or energy/ power [150].

Memory Rules
Memory rules are those that allow models to store information of a certain number of active reversal points, understanding as an active reversal point that which appears in t 0 and can influence the soil behavior for t [ t 0 . Depending on the number of active reversal points from which information is stored, all or part of the recent history of the soil will be considered. Furthermore, the models that store only part of this information may reproduce a certain finite number of symmetric loading cycles without transgressing the 1st Principle of Thermodynamics. It is possible to distinguish three types of models based on the number of reversal points from which information is stored.
Information storage from a single active reversal point These models store information of the last active reversal point and, therefore, only consider the history between such reversal point and the current state, offering important limitations to reproduce hysteretic behavior but requiring little computational memory to store such information. Some models that belong to this group are the following: Simpson et al. [3]; Whittle [142]; Al-Tabbaa and Wood [10]; Bolton et al. [148]; Yu [143]; Puzrin and Burland [152]; Gryczmanski et al. [144]; Pestana and Whittle [149]; Papadimitriou et al. [140].
Information storage of several active reversal points These models store information of a certain number of active reversal points, so they can better reproduce the hysteretic behavior with respect to the previous ones and, therefore, they will have a greater computational memory requirement. Likewise, these models use different typologies of state variables to store information, such as, for example, the situation of yield surfaces in multisurface or bubble models, or state variables in diverse elastoplastic and hypoplastic models. The effect of reversals on these variables depends on the rotation angle of reversal that takes place, and it may be the case that an important reversal erases the effect of previous active reversal points, although these are recent. Some models that belong to this group are the following: Prévost [11,12]; Simpson [31]; Stallebrass and Taylor [123]; Niemunis and Herle [151]; Benz [67]; Cudny and Truty [155].
Information storage of all active reversal points These models can store information of all the active reversal points, improving their capability to reproduce the hysteretic soil behavior, although they will generally require a high computational cost. Nevertheless, for practical purposes, this type of models ends up limiting the number of reversal points from which they store information. The maximum number of reversal points considered will depend on the number of expected loading cycles. Some models that belong to this group are the following: Hueckel and Nova [22]; Niemunis et al. [153]; Schädlich and Schweiger [154]; Castellón and Ledesma [156] and those that fulfill the Generalized Masing Rules [157,158].

Effect of the Reversals on the Variables that Control Stiffness Degradation
The nonlinear models consider shear stiffness degradation with deformation. The mechanisms that control this degradation are identified here with the variables ! i ! 0, which can be considered as part of the variables v hist . These models can be grouped in two categories: The variables ! i that controls the degradation are always fully reinitialized after a reversal These models consider that the variables ! i adopt a value of 0 after a reversal ( Fig. 22) . This implies that in these models the maximum value of stiffness is reached after a reversal, which does not allow simulating the experimental soil behavior described in the work of Atkinson et al. [30] (Fig. 14), Clayton and Heymann [20] (Fig. 15), Gasparre [75] or Gasparre et al. [34] ( Fig. 16). Some models that belong to this group are the following: Hueckel and Nova [22]; Simpson et al. [3]; Whittle [142]; Al-Tabbaa and Wood [10]; Bolton et al. [148]; Yu [143]; Stallebrass and Taylor [123]; Gryczmanski et al. [144]; Puzrin and Burland [152]; Pestana and Whittle [149]; Papadimitriou et al. [140]; Niemunis et al. [153] and those that fulfill the Generalized Masing Rules [157,158].
The variables ! i that controls the degradation are fully or partially reinitialized after a reversal depending on the reversal rotation angle These models consider that the variables ! i reduce their value after a degradation process (Fig. 22) after the reversal R). Depending on the values ! Rþ i different stiffness values are obtained after a reversal. This type of models allows partially simulating the experimental soil behavior described by Atkinson et al. [30] (Fig. 14), Clayton and Heymann [20] (Fig. 15), Gasparre [75] or Gasparre et al. [34] (Fig. 16) Some models that belong to this group are the following: Prévost [11,12]; Simpson [31]; Niemunis and Herle [151]; Benz [67]; Schädlich and Schweiger [154]; Cudny and Truty [155]; Castellón and Ledesma [156].

Effect of the Reversals on Maximum Soil Stiffness
Experimentally, it is observed how the magnitude of shear stiffness recovery depends on the stress/strain reversal rotation angle, as shown in Fig. 23. Considering this aspect of soil behavior, it is possible to classify the models following two different criteria: (1) depending on whether the soil stiffness recovery is continuous or discontinuous with the rotation of stress/strain recent path; and (2) depending on whether such recovery is always total after a reversal or can be partial or total after it. Introducing the second criterion within the first, the models can be classified as follows: Discontinuous recovery of stiffness with the rotation angle of stress/strain recent path This group integrates those models in which the recovery of soil stiffness occurs discontinuously (staggered) according to the rotation angle of the stress/strain recent path. This is because the variables ! i , that control the degradation process of shear stiffness, experience finite jumps in their value for certain values of such rotation (Fig. 23a). These models will be, therefore, incrementally multilinear. In turn, the models that consider a total recovery after a reversal or a recovery that can be partial or total are distinguished in this group (Fig. 23b). Some models with a total recovery after a reversal are the following: Hueckel and Nova [22]; Simpson et al. [3]; Whittle [142]; Al-Tabbaa and Wood [10]; Bolton et al. [148]; Yu [143]; Stallebrass and Taylor [123]; Puzrin and Burland [152]; Gryczmanski et al. [144]; Pestana and Whittle [149]; Papadimitriou et al. [140]; Niemunis et al. [153] and those that fulfill the Generalized Masing Rules [157,158]. Some models with a partial discontinuous recovery according to the rotation angle of the stress/strain after a reversal are the following: Prévost [11,12]; Simpson [31]; Benz [67]; Cudny and Truty [155].
Continuous recovery of stiffness with the rotation angle of stress/strain recent path This group integrates those models in which the recovery of soil stiffness occurs continuously according to the rotation angle of the stress/strain recent path. This is because the variables ! i that control the degradation process of the shear stiffness vary continuously with the values of such rotation (Fig. 23a). These models will therefore be Fig. 22 Effect of reversals on the variables that control degradation, drawn on the graphic of the shear modulus degradation of a London clay, extracted from the work of Atkinson et al. [30] incrementally nonlinear, and the recovery of stiffness is continuous with the rotation of the strain path. Some models that can be considered to have a partial continuous recovery according to the rotation angle of the stress/strain after a reversal are the following: Niemunis and Herle [151]; Schädlich and Schweiger [154]; Castellón and Ledesma [156].

Classification of the Models Using the Tensorial Zone Criterion
The tensorial zone criterion has been used to classify the models that consider soil behavior in the range of small strains (Zones I and II of Jardine [6]), as follows: (A) Incrementally linear models associated to one tensorial zone ( The set of incrementally linear models include the elastic models, although both linear elastic models and nonlinear elastic models can also be formulated non-incrementally, that is, according to r 0 ¼ F e e ð Þ. Finally, it should be noted that the fact of not having the numerical codes of most of analyzed models complicates and limits their study, idea shared by Gudehus [107]. (A) Incrementally linear models The number of tensorial zones #Z considered in the incrementally linear models is reduced to one. These models can reproduce soil nonlinear behavior in Zone II, but not its hysteretic and dependent on recent history behavior.
(A.1) Elastic models Elasticity constitutes one of the fundamental pillars of mechanics of deformable solids.
Two fundamental characteristics of these type of models are: (1) the proportionality between effective stress r 0 and elastic strains e e ; and (2) the possibility of determining the stress state r 0 x; t ð Þ at any point of the continuous medium and at any instant, only from the state of elastic strains e e x; t ð Þ in such point and instant, without the need to know the previous history.
Soil behavior in Zone I can be approximated by a linear elastic model. In addition to that, it is common to find inherent and induced anisotropy in most soils, although it is often complex to separate the influence of each of them on the test results. Several discrepancies in numerical simulations regarding experimental results come precisely from the fact that the effects of anisotropy are not considered. Piriyakul [159] demonstrates the need to consider inherent and induced anisotropy in the range of small strains, and Poulos [160], Simpson [6], Simpson et al. [161], Addenbrooke et al. [68], Zwanenburg [162], Kung et al. [163] and Whittle [36] point out the importance of considering soil anisotropy when elastic models are used. Nevertheless, the difficulty to obtain the parameters associated with anisotropic linear elastic models leads, on multiple occasions, to the use of isotropic linear elastic models.
The elastic stiffness moduli in the range of strains characteristic of Zone I correspond to the dynamic elastic stiffness type that are G ¼ G 0 and K 0 ¼ K 0 0 when considering isotropy. These can be calculated from the soil density q soil and the speed at which the shear waves S (v s ) and compression waves P (v p ) are transmitted through its skeleton, according to the following expressions: (A.1.1) Linear elastic model Linear elastic models use a linear tensorial function F Á ð Þ for the total constitutive equation.
Since the elastic stiffness tensor E 0 is constant ( _ E 0 ¼ 0Þ, the incremental constitutive equation can be expressed as follows: Considering the symmetry of r 0 and e e , and some elastic thermodynamic considerations, the components of E 0 can be reduced from 81 to 21. Anisotropic linear elastic model The most general expression of E 0 for a general anisotropic linear elastic material has 21 independent components. In anisotropic materials, it is common to use E 0 ij , m 0 ij o G ij as elastic parameters. In case the material properties present three planes of symmetry, the material is said to be orthotropic, and the independent components are reduced from 21 to 9. And if the properties of the material have an axis of symmetry, it is said that it has transversal isotropy and the independent components are reduced from 21 to 5. Orthotropic elastic models are rarely used. In contrast, elastic or hypoplastic models that consider soil transversal anisotropy have been used in professional practice to reproduce soil behavior in the range of small strains [14,154,159,164].
Isotropic linear elastic model In case the material properties present three symmetry axes, the material is said to be isotropic and the independent components are reduced from 21 to 2. The constitutive equation of an isotropic linear elastic model can be expressed as follows: On the other hand, considering s ¼ r 0 À r 0 e e k k, it is possible to uncouple the volumetric and deviatoric behavior of the soil in the expression (19), obtaining the expressions r 0 oct ¼ 3K 0 e e oct and s oct ¼ Gc e oct . Linear elastic models: • Are linear models.
• Do not distinguish between stiffness in loading and unloading/reloading and, therefore, cannot reproduce the hysteretic behavior of the soil.
• Do not generate energy in closed stress or strain cycles, H e e r 0 : _ e e ¼ 0.
(A.1.2) Nonlinear elastic models The classification proposed by William [165], based on the structure of the tensorial function F Á ð Þ that relates stress to elastic strain, is used to define the nonlinear elastic models. The hypoelastic models stricto sensu are not elastic models per se. As previously commented, one of the fundamental properties of the elastic models is the possibility to determine r 0 x; t ð Þ only knowing e e x; t ð Þ, without the need of knowing the previous history. On the contrary, in hypoelastic models, the stress state at a given point and instant does depend on such previous history. Despite this, the hypoelastic models have been traditionally considered within the framework of elastic materials, therefore, such classification has been respected in this work.
Algebraic formulation (Cauchy elastic models and pseudoelastic models) Cauchy elastic models Cauchy elastic models use a nonlinear tensorial function F Á ð Þ for the constitutive equation of the material.
A priori, and given the arbitrariness in the functions / i , this type of models can generate energy during the application of cyclic loading, thus transgressing the thermodynamic principles.
The linearized constitutive incremental equation is expressed as follows: Pseudoelastic models Pseudoelastic models are a particular case of Cauchy elastic models. In this case, the nonlinearity in F Á ð Þ is incorporated through the secant elastic stiffness tensor. As with linear elastic models, if an isotropic material is considered, the secant elastic stiffness tensor E 0 s depends on two parameters, being able to be expressed according to the drained secant elastic bulk modulus K 0 s and to the secant elastic shear modulus G s , which leads to the known K À G nonlinear models formulation.
In the same way it was deduced for the constitutive equation of an isotropic linear elastic material, it is also possible in this case to uncouple the volumetric and deviatoric behavior from expression (22) À Á e e oct and s oct ¼ G s c e oct À Á c e oct . In case of considering transversal anisotropy, it is necessary to introduce the secant coupling modulus J s into the previous formulation, resulting r 0 oct ¼ 3K 0 s e e oct þ J s c e oct and s oct ¼ J s e e oct þ G s c e oct . The linearized constitutive incremental equation of the pseudoelastic models is the following: The above expression can be rewritten as t can be written as follows: In the previous expression K 0 t and G t are the tangent moduli which, considering À Á , can be expressed as: Based on the above and considering that e e : _ e e ¼ 3=4c e oct _ c e oct , the constitutive incremental Eq. (23) becomes expression (27), which is possible to express uncoupled as _ Cauchy elastic models and pseudoelastic models: • Are nonlinear models.
• Do not distinguish between stiffness in loading and unloading/reloading and, therefore, cannot reproduce the hysteretic behavior of the soil. • Do not consider the effect of recent history on soil stiffness. • Constitute simple formulations from a conceptual and mathematical point of view. • The parameters of the models are, in general, easy to obtain from simple tests. • In stress/strain closed cycles, they conserve stress and • Depending on the selection of soil parameters, they can generate energy in stress/strain closed cycles, thus transgressing the thermodynamic principles, H e e r 0 : _ e e 6 ¼ 0. If a Cauchy elastic model or a K À G conservative, since, under these conditions, the elastic energy W el is independent from the followed path, complying with W el ¼ 0 in a closed cycle, where Integral formulation (hyperelastic models) The hyperelastic models define an energy potential from which the expression of F Á ð Þ is derived. This potential is called elastic strain energy W ¼ W e e ð Þ. In the case of isotropic materials, W e e ð Þ depends on the invariants of r 0 and e e .
From expression (28), it is verified that this type of models complies with the 1st Principle of Thermodynamics, thus not allowing energy to be generated in stress/strain closed cycles (expression 29).
Hereunder, the most general expression of F Á ð Þ for isotropic materials is derived. For that purpose, the strain tensor invariant momentsÎ e e i ¼ 1=i ð Þtr e e ð Þ i À Á are defined. Using the chain rule and the notation W i e e ð Þ ¼ oW e e ð Þ=oÎ e e i , the expression (30) is obtained: As it can be seen, expression (30) has the same form as expression (20). Nevertheless, the coefficients W i e e ð Þ, unlike / i e e ð Þ, satisfy the theorem of Schwartz, which implies the integrability conditions o 2 W e e ð Þ=oÎ e e i oÎ e e j ¼ o 2 W e e ð Þ=Î e e j oÎ e e i , from which it follows that the tangent elastic stiffness and compliance tensors have greater symmetry in hyperelastic models and it can be shown that this is a necessary condition to satisfy the 1 st Principle of Thermodynamics. These conditions do not have to be complied in the Cauchy elastic models or in pseudoelastic models.
The constitutive incremental equation in the hyperelastic models is expressed as follows: Several authors developed hyperelastic models that consider the relation G ¼ G p 0 ð Þ [166][167][168][169][170][171]. Such models, in the case of isotropic materials, adopt the following uncoupled form Hyperelastic models: • Are nonlinear models.
• Do not distinguish between stiffness in loading and unloading/reloading and, therefore, cannot reproduce the hysteretic behavior of the soil. • Do not consider the effect of recent history on soil stiffness. • Introduce a great number of parameters, which are often difficult to obtain from simple tests. • In stress/strain closed cycles, they conserve stress and • Satisfy the theoretical requirements of continuity, stability, uniqueness and the 1 st Principle of Thermodynamics, H e e r 0 : _ e e ¼ 0.
The hyperelasticity framework has allowed the subsequent development of hyperelastic-plastic models. This is the case of the model of Likitlersuang and Houlsby [172], from which it is possible to deduce multisurface models with kinematic hardening, or the model of Zhang et al. [173], that is able to consider, among others, soil cementation, stress-induced anisotropy, cyclic shear behaviour or shear modulus degradation under small strain conditions due to the fact that energy dissipations and irreversible strains are allowed within all strain ranges.
Differential formulation (hypoelastic models stricto sensu) Incrementally linear hypoelastic models stricto sensu consider the existence of tensorial functions that relate the stress increments with total strains, total stresses and strain increments. Such models do not directly define the tensorial function F Á ð Þ, but this results from integrating the incremental constitutive Eq. (32), leading to expression (33).
where E 0 t r 0 ð Þ, in its most general form, depends on 12 functions C Ã i , which in turn depend on the three invariants of r 0 , as Unlike what happens with the algebraic formulation and the integral formulation exposed in the previous sections, the recent history in the constitutive law of the material is considered in the hypoelastic formulation. This is why hypoelastic models stricto sensu are not elastic models per se.
• Do not distinguish between stiffness in loading and unloading/reloading and, therefore, cannot reproduce the hysteretic behavior of the soil. • The extension of these models led to the appearance of quasi-hypoelastic models, framed within the incrementally multilinear models (see Sect. B.1.2). Quasihypoelastic models do allow the hysteretic behavior of the soil to be reproduced, since they introduce different stiffness depending on whether the processes are of loading or unloading type, and define the reversal criteria, the memory rules and the effect that reversals have on the variables which control the degradation of stiffness and on the maximum soil stiffness.
• Depending on the selection of soil parameters, they can generate energy in stress/strain closed cycles, thus transgressing the thermodynamic principles, H e e r 0 : _ e e 6 ¼ 0.
• Simulations with these types of models should be limited to the stress and strain paths corresponding to the tests with which their parameters were obtained.
(B) Incrementally multilinear models The number of tensorial zones #Z considered in the incrementally multilinear models is a finite number greater than one. These models can reproduce nonlinear, hysteretic and dependent on recent history soil behavior, characteristic of Zone II.
(B.1) Hysteretic models Hysteretic models emerged with the aim of simulating hysteretic soil behavior in the range of small strains originally in dynamic problems. Assuming that the hysteretic behavior of the soil was perfect, deformations would be recoverable only if the stress/strain cycles that started from the same reversal point were closed. This behavior differs from the reversible behavior of elastic models and from the irreversible behavior of elastoplastic models, as can be seen in Fig. 25.
Nevertheless, reality does not show a perfect hysteretic soil behavior. Indeed, Hueckel and Nova [22] point to the elements described in Sect. 3.3 that separate the real soil behavior from such idealization (Fig. 26).
(B.1.1) Paraelastic models A particular type of hysteresis is the so-called paraelastic hysteresis or paraelasticity, a theory that was developed by Hueckel and Nova [22], being one of the fundamentals in this field.
The term paraelastic refers precisely to the fact that strains reversibility is conditioned to the closure of the cycles (Fig. 27). The paraelastic models consider that soil stiffness degrades as a function of variables that depend on the paraelastic strains e pe ¼ e e þ e mp , sum of the elastic strains e e and the microplastic strains e mp , the latter reversible only under the condition of the corresponding stress/ strain cycle closure. The cycles of this kind of models can overlap one another using simple memory rules with two information levels, 1st for the constitutive equation and 2nd for recent history. Table 2 analyzes the stiffness formulation, the treatment of hysteretic behavior and general aspects of some paraelastic hysteretic models.
The clear, robust and intuitive way in which the paraelastic theory defines the hysteretic behavior of the soil, in addition to its versatility to be combined with other theoretical frameworks, endows it with a high potential. So far, interesting models have been developed using concepts Fig. 27 Idealization of the hysteretic behavior of paraelastic models, from [22]. a One simple cycle. b One minor cycle overlapped to a major cycle. c Closure of minor overlapped cycles from paraelastic theory combined with elastoplastic, bounding plasticity, multilaminate or hypoplastic models [140,153,154,156,[174][175][176][177].
(B.1.2) Quasi-hypoelastic hysteretic models Among the incrementally multilinear models stand out the ones here called quasi-hypoelastic hysteretic models, which have been normally used to simulate soil behavior in the range of small strains. A widespread methodology to build quasi-hypoelastic hysteretic models is to incorporate the elements that allow reproducing the hysteretic behavior of the soil into quasi-hypoelastic models, for example, through the Generalized Masing Rules [157,158]. Quasihypoelastic models simulate soil behavior during monotonous loadings and resemble hypoelastic models, although instead of elastic stiffness moduli they use apparent stiffness moduli which act on total strains (and no just on the elastic ones). Following the classification by Sellers [124], quasi-hypoelastic models can be classified into two types depending on how their parameters are defined: • Variable parameter models [178][179][180] are those that propose a relation between stress and strain based on nonlinear functions (hyperbolic, elliptical, logarithmic, spline, etc.). Variable parameter models can generally be reformulated as variable moduli models. • Variable moduli models [181][182][183][184] are those that propose nonlinear functions to describe apparent volumetric and shear stiffness moduli in a constitutive equation with an algebraic structure equal to that of an elastic model.
The general constitutive equation of the variable moduli models usually adopts the following expression: The following decomposition of secant stiffness modulus is usually considered based on experimental observations [66]: The popularity of quasi-hypoelastic hysteretic models come from its relative simplicity and easy incorporation into finite element numerical codes for solving practical engineering problems [66], as can be seen in the work of Jardine and Potts [185] and St. John et al. [186]. Some disadvantages in relation to the quasi-hypoelastic hysteretic models are pointed out [66,182,187]: • Lack of continuity during neutral loading processes.
are calibrated from tests with specific stress or strain paths that can be very different from paths followed by the soil in practical applications.
1. Reversal criterion:intrinsic,  • The constitutive equation implies that stress and strain increases are coaxial, which does not always correspond to experimental observations [188]. Nevertheless, it should be considered that in the tests used to calibrate these models this coaxiality will generally be complied. • To reproduce the hysteretic behavior of the soil, loading/unloading criteria that can lead to numerical instabilities are introduced. • The dependence of G ap s;0 p 0 ð Þ implies that the shear modulus depends on the volumetric strain G ¼ G c oct ; e oct ; v hist À Á which, in case of not considering the coupling modulus J t , may not be thermodynamically consistent.

Generalized Masing rules
A common and widespread way to build hysteretic models from models that simulate the nonlinear behavior of soil during monotonous loading is through the application of the original and extended Masing Rules [157,158], jointly denominated as Generalized Masing Rules, to which the corresponding loading/unloading criteria must be added. Such rules are formulated in a one-dimensional stress or strain context, although the extension to multiaxial states can be made with an appropriate change in the variables. For the one-dimensional case, the Generalized Masing Rules consider that the state variables v hist are identified with the total strain e. However, despite these variables v hist can be formally treated as strains, they cannot be physically interpreted as such (see Sect. 3.3). The original Masing rules Nr. 1 and Nr. 2 [157] for symmetric loading, along with the extended Masing rules Nr. 3 and Nr. 4 [158] are the following: 1. Rule Nr. 1: For the primary loading, the constitutive law adopts the expression r 0 ¼ F e ð Þ. In general, this relation considers the total strain e. This is because in the tests used to calibrate the nonlinear function F Á ð Þ, the reversible and irreversible strain components are not separated. 2. Rule Nr. 2: In loading or unloading, after a reversal point R, the constitutive law considers a two-scale factor, as compared to the constitutive law for primary loading, taking as reference the local origin given by the reversal point r 0R ; e R ð Þ, as r 0 À r 0R ¼ 2F e À e R ð Þ=2 ð Þ . In general, this relation also considers the total strain e, although in elastoplastic models, this will coincide with the elastic strain, as this law is activated within elastic domain. Rule Nr. 2 can also be applied through different strategies [158,189].
3. Rule Nr. 3: If the unloading or reloading curves intersect the initial loading curve, they resume such curve. 4. Rule Nr. 4: If the unloading or reloading curves intersect a previous unloading or reloading curve, they resume such curve.
Some quasi-hypoelastic hysteretic models Some models that are classified as quasi-hypoelastic hysteretic models are described and analyzed hereafter. They all are variable parameter models, except for the SSOM model [67] and the HQH model [156,175], which are variable moduli models. Nevertheless, practically all of them can be reformulated as variable moduli models, which allows obtaining the equivalent apparent stiffness moduli in each case. In this type of models, the state variables that control the value of the stiffness moduli are identified with those proposed by Gudehus [107] to describe the hysteretic behavior of the soil.
In all the models studied in this section, excluding the SSOM and HQH models, it is necessary to define ad hoc components that allow reproducing hysteretic soil behavior for example, through the Generalized Masing Rules. Tables 3 and 4 analyze the small strain stiffness formulation and general aspects of some quasi-hypoelastic hysteretic models.
The SSOM [67] and HQH [156] models have great potential to be used in professional practice, given their ability to integrate multiple aspects of the soil in the range of small strains, in addition to their capacity to be combined with plastic models which let extend its applicability to the range of medium and large strains.
Both models: (1) are based on the apparent secant shear modulus G ap s degradation curve provided by the Dos Santos and Correia model [231], which is limited by a minimum value of the tangent shear modulus G ur ; (2) consider the recent history of deviatoric strains; (3) use history tensors that act as Simpson block models [31]; (4) consider similar thermodynamic-type corrections; and (5) can be combined with multiple plastic models leading to advanced elastoplastic models.
When both models are used by themselves, they use the Hashiguchi [189] criterion to differentiate primary load from unloadings/reloadings, but when they are combined with an elastoplastic model, they consider, on the one hand, a factor n ¼ 2 in the expression of the apparent secant shear modulus G ap s and, on the other hand, a factor h i that modifies the hardening laws to avoid overlapping the mechanisms of these models and the plastic models with which they are combined, that try to explain the reduction of soil stiffness during the primary loading.
The SSOM model describes the elastic part of the wellknown HS-SS model, implemented in Plaxis and based on Based on Ramberg and Osgood [190] model according to Idriss et al. [191] and Hara [192]. Numerical examples can be found in Ishihara [193] Kondner and Zelasco [194] E 0ap s ¼ 1 aþbe1 Hyperbolic model initially developed by Kondner [195]. Parameter setting can be found in Duncan and Wrong [196] Hansen Hyperbolic model based on previous models [194,195]. It introduces the factor R f . Duncan et al. [71] proposed an expression for K 0ap s and E 0 t;ur . Kulhawy et al. [198] proposed an expression for m 0 t : The model has been successfully applied in various boundary value problems [71,199,200] Desai [179] Spline functions Spline functions [201] are used to set experimental data with nonlinear curves that represent constitutive relations. Some applications can be found in Desai [179] and Yniesta et al. [202]. Despite the precision in curve fitting, the computational cost is high Hyperbolic model to reproduce soil behavior subjected to static and dynamic loadings under small and large strains Jardine et al. [8] 3 G ap Generalization to the multiaxial case of the model of Jardine et al. [203] Tatsuoka and Shibuya [115] Y ¼ Hyperbolic model for small, medium and large deformations. It is applicable to sands under compression, constant confinement and plane strain conditions. Applications can be found in Chiaro et al. [205], Tatsuoka and Shibuya [206], Siddiquee [207] and Siddique et al. [208,209] Matasovic and Vucetic [210] G ap s ¼ G0 Based on Kondner and Zelasko model [194] and in the work of Matasovic [211].
Comparison with the hyperbolic model is done in the work of Stewart et al. [212] Ishibashi Based on previous works such as those of Gunn [223] or Bolton et al. [224]. The model provides a good approximation in the simulation of pressuremeter tests [222] Lee and Salgado Subsequently developed by Lee and Salgado [226,227], it is a generalization to the multiaxial case of the model of Fahey and Carter [219]. Some applications can be found in the work of Lee et al. [228] Small Strains in Soil Constitutive Modeling 3247 the HS-S model [67] (see Sect. B.2.1). However, despite its great advantages, the SSOM model presents some limitations that the HQH model solves, with the following con- Atkinson model is applicable to drained or to undrained loading, with appropriate values for the parameters. For typical values of stiffness and degree of non-linearity for soil the value r 0 is generally in the range 0:1 to 0:5 [32] Darendeli [229] G ap Based on the model of Hardin and Drnevich [2]. The model is applied to both cohesive soils and granular soils [229,230] Dos Santos and Correia [231] Based on the model of Hardin and Drnevich [2]. The value of the parameter a has been selected through a least square adjustment    , taking K values of 0:40 À 0:58 for sands and of 0:50 for clays [156] (K % m T =m R , where m T and m R are the parameters that control stiffness in the Niemunis and Herle model [151]).
(B.2) Advanced models The advanced models are those incrementally multilinear models that are capable of reproducing soil behavior in the entire strain range (Zones I, II, III and IV of Jardine [6]).
(B.2.1) Classic elastoplastic models Classic elastoplastic models that consider soil behavior in the range of small strains add, to the elastic regime of the classic elastoplastic formulation, the nonlinear, hysteretic and dependent on history behavior of the soil through the incorporation of new state variables v el;hist to the constitutive equation within yield surfaces _ r 0 ¼ G e; r 0 ; v el;hist ; _ e À Á . This is normally attained by adopting the structure of a quasi-hypoelastic hysteretic model to describe the elastic formulation of the elastoplastic model. In this case, to calculate the stress, the apparent stiffness moduli are operated with the elastic strains and not with the total ones. According to Jardine [6] and Jardine et al. [8], this semiempirical methodology allows attaining very precise predictions thanks to the wide range of strains considered, although its use entails a loss of generality and a lack of consistency of the theoretical framework on which it is based. Among others, these assumptions lead to the following: (1) it is not correct to apply the elastic relations to the apparent stiffness moduli; and (2) the mechanisms of the quasi-hypoelastic hysteretic models and the elastoplastic models that try to explain the reduction of soil stiffness during the primary loading processes, in comparison with the stiffness in unloading or reloading processes, are being overlapped. The basic equations used in classic elastoplastic models, generalized to the case of k ¼ 1. . .q yield surfaces are widely known.
The number of tensorial zones #Z that present these models will depend on the number of yield surfaces and their intersections. In general, it will be complied that #Z ¼ 2 in models with a single yield surface and #Z [ 2 in models with multiple yield surfaces. Table 5 analyzes small strain stiffness formulation, hysteretic behavior treatment and general aspects of some elastoplastic models that consider soil behavior in the range of small strains.
The widespread use of the elastoplastic formulation and its balance between complexity and usability, gives the models that are based on it the greatest potential to be used in the practical applications, among all those analyzed in this paper, of which stand out the HS-SS model (based on the HS-S model [67]), and the EPHYSS model [156]. The incorporation of formulations that can reproduce soil behavior in the range of small strains within models used in professional practice (as it was achieved with the HS-SS model in Plaxis), would allow professional engineers to familiarize themselves with the characteristics and complexity of non-linear, hysteretic and dependent on recent history soil behavior in the range of small strains. Once this first goal has been achieved, introducing more complex models into commercial software could be less traumatic.
The HS-SS model is one of the few models widely used in professional practice that considers soil behavior in the range of small strains and it uses known and relatively easy to obtain/estimate parameters. The HS-SS model is based on the HS-S model which adds, to the elastic formulation of the SSOM, the plastic formulation of the HS model [233] with two modifications: (1) it replaces the dilatancy criterion of Rowe [238] with that of Li and Dafalias [239] to describe the contractive behavior of the soil; and (2) modifies the hardening laws and introduces the factor The HS-S model significantly improves the approach to soil behavior provided by the HS model, although it has certain limitations [67]. Some applications can be found in [240][241][242].
The HS-SS model presents some small differences with respect to the HS-S model on which it is based. These differences appear in the expressions both models use to calculate the dilatancy angle and in the expression of the Cap-surface hardening law, where the increase of the state variable depends on the confinement in the HS-S model and on the state variable itself in the HS-SS model.
Despite its great advantages, some aspects of the HS-SS model could be improved, and some inconsistencies whose effects can have a considerable influence on the numerical simulations of boundary value problems, as these are cumulative, should be corrected [156,175,[243][244][245]. This is achieved with the EPHYSS model. The EPHYSS model adds to the elastic formulation of the HQH model, the plastic formulation of the HS MOD model [156], which is the same as that used in the HS-S model and, like this one, also uses known and relatively easy to obtain/estimate parameters. EPHYSS model improves the approach to the soil behavior in the range of small strains with respect to the HS-SS model (see Sect. B.1.2) and corrects the numerical inconsistencies detected in it. Simulations with EPHYSS model of different element tests and several boundary value problems can be found in Castellón [175].
Although unsaturated soil models are not analyzed in this paper, it should be pointed out that recently some advanced elastoplastic models for unsaturated soils capable of reproducing soil behavior in the range of small strains have been proposed, such as, for example, the model of Ng et al. [246].
(B.2.2) Multisurface models The multisurface models with kinematic hardening were initially developed by Iwan [247] and Mróz [13], based on previous works, as those of Duwez [248] and Besseling The plastic part is described with the MCC model 3. Effect of reversals on degradation variables: total reinitialization after reversals 4. Effect of reversals on soil stiffness: discontinuous stiffness recovery In the elastic zone, it is assumed G ap t ¼ G 0 for e q \10 À5 , and G ap t ¼ G ap t e q À Á for 10 À5 \e q \10 À2 . The elastoplastic formulation of the MCC is assumed when e q [ 10 À2 . Some applications of the SDMCC model can be found in Osman and Bolton [232] HS-S [67] Same elastic formulation as SSOM. The plastic part is described with HS model [233] Idem  [249]. Subsequently, Prévost [11,12] extended these models and applied them to Soil Mechanics. These models consider various nested yield surfaces (f m 00 ð Þ ¼ 0; m 00 ! 2), generally in the stress space (Fig. 28a). Each of these surfaces is activated at the moment when the stress state reaches its boundary, leading to a finite decrement in the global plastic hardening modulus, which, at the same time, generates greater plastic strains (Fig. 28b).
Multisurface models can reproduce a discontinuous variation of elastoplastic stiffness. When j nested yield surfaces are activated, and depending on the followed stress path, they are displaced keeping them tangent at the contact points. The number, form, size and situation of the nested surfaces constitute the state variables v hist . The approximation to the soil behavior that these models provide depends on the number of nested yield surfaces considered, since a greater number of surfaces implies a softer variation in stiffness. However, if this number is high, using these models can be computationally expensive. Based on this limitation, bounding surface plasticity models and bubble models were created.
The multisurface models have two tensorial zones #Z [66]. Table 6 analyzes small strain stiffness formulation, hysteretic behavior treatment and general aspects of some multisurface models.
Puzrin and Burland [152], after studying and evaluating the models of Kondner [195], Hardin and Drnevich [2], Griffiths and Prévost [253], Tatsuoka and Shibuya [115] and Jardine et al. [203], and based on the works of Puzrin and Burland [254], developed a model based on the theoretical framework of Jardine [6], distinguishing, in the stress space, three zones where the soil has different behaviors around any local origin constituted by a point of the stress path. Addenbrooke et al. [68] proposed the generalization of the model of Puzrin and Burland to the multiaxial case and Puzrin et al. [255] extended the original model, introducing the plastic behavior in a thermodynamically consistent formulation. However, Puzrin and Burland [152] model requires unusual parameters in the professional practice, which may be a handicap for its use. On the other hand, this model has certain limitations to reproduce the hysteretic behavior of the soil, since: (1) it only stores information of the last reversion point; (2) after a reversion the state variables that control degradation can only be fully reinitialized; and (3) stiffness recovery is discontinuous with the rotation of the recent path.
The approach proposed by Seyedan and Solowski [252] seems more promising for future practical applications, due to its ability to introduce soil behavior in the range of small strains in classical models. Seyedan and Solowski developed a methodology to introduce shear stiffness nonlinearity to classic elastoplastic models that have a constant elastic shear modulus. A new Small Strain Yield Surface is defined to reproduce the small strain shear stiffness nonlinearity following the expression from Dos Santos and Correia [231]. This methodology increases the total number of parameters of the original model by only two, and does not affect its original plastic formulation. Seyedan and Solowski apply this methodology to the MCC and MC models. State variable a q replaces e p q in the calculation of the hardening law to include errors caused by G e MCC . Once G reduces its value to G e MCC , it is maintained within an elastic behavior until stress state reaches MCC Yield Surface or a reversal occurs. The movement of the Small Strain Yield Surface does not produce any volumetric plastic strain and, therefore, does not affect the MCC Yield Surface, but the deviatoric plastic strain of the MCC Yield Surface does affect the Small Strain Yield Surface. In Seyedan and Solowski [256] it is showed how F ss is well Fig. 28 a Schematic representation of a multisurface model. b Experimental curve q À e 1 (left) and curve q À e 1 resulting from multisurface model (right). From [250] Small Strains in Soil Constitutive Modeling 3251  [11,12] Linear elastic model within f m 00 ð Þ ¼ 0, Puzrin and Burland [152] Within LER (linear elasticity): Within SSR and outside LER (nonlinear elasticity):

.3) Brick models
Most of the elastoplastic models for soils are formulated in the stress space, although some have been formulated in the strain space [31,67,257,258]. Yoder and Iwan [259] and Iwan and Chelvakumar [260] highlight some of the advantages of the models that are formulated in the strain space. Many of them have not prospered, but the brick model of Simpson [31] has led to a theoretical framework that has subsequently been successfully applied to various models. Among them are the model of Clarke [69], the model of Ellison [261] or the HS-Brick Model [155].
To explain how the brick models work, Simpson [31] uses an analogy in which he considers that soil behavior in the range of small strains resembles a man pulling a set of bricks attached to him by inextensible strings. The movement of the man represents the imposed strains path and the movement of the bricks represents the plastic strains. When any of the strings is not tensioned, the part of the plastic strain corresponding to such string is not generated. But when such string is tensioned, as the man is at a distance equal to that of the length of the string and in its direction of stretching, the respective plastic strains are activated, with the consequent decrease in the elastoplastic stiffness of the soil (Fig. 29). The number and situation of the bricks constitute the state variables v hist .
To explain the advantages of working in the strain space, Ellison et al. [262] state that, despite the connection between stress and strain, the changes in soil behavior are determined by micromechanical changes in its fabric, which is reflected at the macroscopic level as an accumulation of strains. Furthermore, they affirm that anisotropy due to the orientation of soil particles changes when large strains occur, even if the stress is the same at the beginning and at the end of these strains. Finally, they remember the well-known statement of Burland (1967) cited by Simpson [31]: ''stress is a philosophical concept-deformation is the physical reality''. Ellison et al. conclude that the models formulated in the strain space are more suitable to explain the mechanisms that govern soil behavior. This approach allows explaining, in a relatively simple way, the observations made by Atkinson et al. [30], Clayton and Heymann [20], Gasparre [75] and Gasparre et al. [34]. From a practical point of view, the models formulated in the strain space can be easily implemented in numerical models without the need to isolate, in each calculation step, the strain as an independent variable [260], just as the models formulated in the stress space should do. In addition, as Yoder and Iwan [259] demonstrate, the models formulated in the strain space do not have the problem of surface intersection, nor some of the instability problems that the models formulated in the stress space do.
Brick models consider as many tensorial zones as possible combinations of active and non-active bricks exist. If n bricks are considered, the maximum number of tensorial zones can be calculated as the sum of the possible combinations with repetition of 1. . .n 00 bricks, plus the case in which no brick is active, that is, n 00 != k! n 00 À k ð Þ! ð Þ . Table 7 analyzes stiffness formulation, hysteretic behavior treatment and general aspects of some brick models.
The original model of Simpson [31] is formulated in the strain space, considering a plane strain state. A finite number of bricks is required to define the degradation curve. This model introduces the effect of soil overconsolidation through the parameter b mod . A higher proportion of volumetric and shear elastic strains is assumed for higher mean stress. It reproduces well the undrained trajectories, and worse the drained ones, unless the effect of dilatancy is introduced. It does not reproduce the critical state, therefore it should be used in stiff clays or its formulation should be extended. The Simpson model [31] formulation that reproduce the behavior of the soil in the range of small strains can be used with relative ease in classic or advanced models (for example, in the HS-Brick Model [155], which improves the HS-SS model). The ability of Simpson model to be combined with other models provides it with great potential. Moreover, the Simpson model establishes a theoretical framework that allows to easily understand part of the formulation of other models, as happens with the history tensors of the advanced elastoplastic models HS-SS and EPHYSS (see Sect. B.2.1).
(B.2.4) Bounding surface plasticity models Bounding surface plasticity models were initially developed by Dafalias [263], Krieg [264], Dafalias and Popov [265] and Dafalias and Herrmann [266] with the aim of attaining a continuous evolution of the elastoplastic stiffness, which is extended to the interior of the yield surface in the range of small strains in these models, and of reducing, therefore, the computational cost of the multisurface models. Bounding surface plasticity models generally use two surfaces: an exterior surface, named bounding surface, and an interior surface, called current loading surface (Fig. 30). Between both surfaces, unlike the multisurface models, which consider m 00 surfaces that give place to a discrete variation in stiffness, bounding surface plasticity models consider continuous functions to represent the plastic hardening modulus, obtaining continuous variations of stiffness (which is equivalent to consider m 00 ! 1 in multisurface models). The form, situation and orientation of both surfaces, along with other information stored at the last reversal point, constitute the state variables v hist .
The formulation of these models is based on that of the classic elastoplastic models, although unlike the latter, the bounding surface plasticity models consider plastic strains also during loading processes within the bounding surface (OCR [ 1Þ. Such plastic strains are made to depend on the distance between the current stress state, located on the current loading surface, and the denominated image point, located on the bounding surface plasticity, according to the mapping rules. The number of tensorial zones #Z of these models is reduced to two when there is only one bounding surface and more than two in case there is a greater number of bounding surfaces [66] (similar to the case of classic elastoplastic models with yield surfaces). Table 8 analyzes small strain stiffness formulation, hysteretic behavior treatment and general aspects of some bounding surface plasticity models.
Although these models can reproduce different aspects of the soil in a wide range of strains, including also small strains, they require parameters difficult to obtain experimentally and with which geotechnical engineers are not familiar. During the coming years it is very possible that these models will remain limited to an academic use.
(B.2.5) Bubble models Bubble models were developed to solve some of the limitations posed by the bounding surface plasticity models. The latter allow reproducing the nonlinear plastic isotropic hardening or softening behavior of the soil during loading processes and assume that part of the deformations corresponding to these loading processes are plastic, even in the range of small strains. Conversely, they consider that soil behavior is elastic during unloading processes, which restricts the degree of coupling between volumetric and deviatoric behavior [84] and limits their ability to reproduce the hysteretic soil behavior in the range of small strains. Bubble models are based on the concepts of bounding surface plasticity models and multisurface models and include one or more surfaces, called bubbles (Fig. 31), with kinematic hardening, in addition to the isotropic hardening within the bounding plasticity surface. The number, shape, size and situation of the bubbles constitute the state variables v hist . Likewise, bubble models must include, apart from the mapping rules that control the plastic hardening modulus as the bubbles move, translation rules that control the movement of such bubbles.
The number of tensorial zones #Z of these models, as in multisurface models, is two. Table 9 analyzes small strain stiffness formulation, hysteretic behavior treatment and general aspects of some bubble models. Other bubble models have introduced the effect of soil destructuration [294,295], soil anisotropy [296] and some improvements to the Three Surface Kinematic Hardening model (3-SKH) [297]. The ability of bubble models to reproduce the hysteretic behavior of the soil is potentially greater than that of bounding surface models. However, bubble models use has been very limited in practical applications, and present a great sensitivity to the exponent that appears in the expressions of the plastic modulus, which has a significant effect on the stress-strain curves [298].
(B.2.6) Multilaminated models The concept of multilaminated model was initially formulated by Batdorf and Budiansky [299]. Zienkiewicz and Pande [300] extended the concept to the analysis of joints in rocks. Subsequently, Pande and Sharma [301] and Sadrnejad and Pande [302] extended these models to Soil Mechanics. Multilaminated models are based on the concept that the three-dimensional stress or strain state at a point of the continuous medium (global state) can be obtained as the sum of the two-dimensional stress or strain state on different i planes with different orientations that pass through such point (local state). For that purpose, a set of fundamental equations are established in each plane: The number of planes considered in multilaminated models is reduced to a finite number and, through the formulation adopted in the model, they are able to represent all orientations in space, grouped by sectors. Bazant and Oh [303] and Ehret et al. [304] demonstrate that using 33 planes in the calculation of _ offers acceptable results and an acceptable computational cost. Global strain increments, which correspond to an increment in global stress, are calculated from the integration of local strains into the planes. For the deduction of the expression that relates these magnitudes, the principle of virtual works is used 2. Memory rules: stores information of the last reversal point (e v;rev , E rev ), which also determine the size of the first yield surface 3. Effect of reversals on degradation variables: total reinitialization after reversals 4. Effect of reversals on soil stiffness: discontinuous stiffness recovery Based on the MIT-E1 model [268], and subsequently developed by Whittle [269] and Whittle and Kavvadas [270]. Ganendra [267] and Whittle [271] analyze discrepancies in the previous formulations of the model and propose some alternative formulations. Several applications can be found in the work of Whittle [36,269] and Whittle et al. [272] CASM-c [143,273] Plastic modulus during primary loading: Plastic modulus during unloading: Plastic modulus during reloading: Based on critical state models [73,74] and on sand models that include the concept of state parameter [274,275].
The CASM-c model simulates nonlinear soil behavior and differentiates loading and unloading stiffness so it can reproduce soil hysteretic behavior. Some applicacions can be found in Khong [276] NAHOS [144] Plastic modulus:  [277], Jastrzebska [278,279] and Sternik [280]. Purely elastic soil behavior is reduced to a single point in the stress space, which agree with the observations of Hueckel and Nova [22]. Gryczmanski et al. [144] comment on the need to introduce a nonlinear elastic law to improve the predictions of the cyclic behavior of the soil. Some applications can be found in the work of Jastrzebska and Sternik [281] MIT-S1 [149] 2. Memory rules: stores information of the last reversal point (e v;rev , e rev ), which also determines the size of the first yield surface 3. Effect of reversals on degradation variables: total reinitialization after reversals 4. Effect of reversals on soil stiffness: discontinuous stiffness recovery Based on the works of Whittle [142], Whittle [269], Whittle and Kavvadas [270] and Pestana [282], MIT-S1 model is as a generalization for clays and sands of the MIT-E3 model [142]. It introduces the void ratio as a new independent variable with the aim of reproducing the characteristic contractive/dilating behavior of sands at different confinement pressures. As well, it considers the hysteretic, nonlinear and plastic anisotropic behavior of the soil, and the critical state behavior of clays. Some applications can be found in the work of Whittle [65] and Nikolinakou et al. [283] and assuming dS 0 % DS 0 ¼ sinã ð Þ sinb DãDb where Dã and Db are finite spherical sectors, ð Þ are satisfied. In this type of models, the state variables that control the value of stiffness modulus in each plane, the number of considered planes and variables that store information of reversion points are identified with the state variables v hist .
The number of tensorial zones #Z in multilaminated models will depend on the amount of planes considered. Nevertheless, this number is generally high; therefore, despite not being incrementally nonlinear models, for practical purposes, it can be considered that the stiffness of this type of models does depend continuously on the increments in the total strain tensor _ e. Table 10 analyzes small strain stiffness formulation, hysteretic behavior treatment and general aspects of Schädlich and Schweiger [154] multilaminated model that considers soil behavior in the range of small strains. Schädlich and Schweiger model considers the elastoplastic behavior of the soil using a Cap-Cone type model (based on the model of Wiltafsky [305] and Wiltafsky et al. [306]), which allows describing plastic behavior in each plane. It also considers the following concepts to describe the nonlinear, hysteretic and history dependent behavior of the soil within yield surfaces: (1) inherent and strain-induced cross-anisotropy (using the spectral decomposition of the global elastic compliance matrix [307]); (2) dependence on stiffness with stress; (3) degradation of elastic stiffness with the strain e deg ; and (4) dependence on strains history satisfying the Generalized Masing Rules. Schädlich and Schweiger model is based on the model of Scharinger et al. [308], but it introduces some differences: (1) it uses a linearization for the degradation curve of the shear modulus in each plane instead of the expression proposed by Soga et al. [309]; and (2) it reproduces the effect of loading history using an approach similar to that of Simpson's brick model [31] instead of the linear approximation of the experimental results from Richardson [29].
Multilaminate models can reproduce multiple soil characteristics, including its complex behavior in the range of small strains. To do this, these models use parameters known in professional practice and a clear and not contrived formulation, which places them in an advantageous position compared to other incrementally multilinear Plastic modulus:  [284] and Papadimitriou [285] and has been subsequently developed by Papadimitriou and Bouckovalas [286]. It uses various concepts: (1) plasticity with kinematic hardening; (2) nonlinearity with hysteresis in the range of small strains based on the models of Ramberg and Osgood [190] and Hueckel and Nova [22]; (3) state parameter for sands [275]; and (4) the effect of the soil fabric. It adopts the expression of Hardin [287] for the calculation of G 0 Fig. 31 a Bounding surface and bubble of a bubble model. b Conceptual scheme of the bubble movement of a bubble model. From [84] Small Strains in Soil Constitutive Modeling 3257 models. Even so, its use has been very limited in practical applications, although some commercial software, such as ZSoil, already incorporate basic multilaminate models in their libraries.
(C) Incrementally nonlinear models Incrementally nonlinear models emerged from the generalization of the hypoelastic theory of Truesdell [310] and from the theoretical framework of elastoplasticity [66]. The most important characteristics of incrementally nonlinear models are: (1) the absence of explicit decomposition of strain increments in reversible and irreversible parts; and (2) the continuous nonlinear dependence of the tangent stiffness tensor on the direction of the increment strain tensor. The most general constitutive equation of this type of models is expression (36) [136,311] and it is composed of a linear part L r 0 ; v ð Þ : d and a nonlinear part d À1 Q r 0 ; v ð Þ :: d d ½ : It is considered that incrementally nonlinear models have an infinite number of tensorial zones #Z ! 1. In some cases, nevertheless, as in hypoplastic models with intergranular strain d, the incremental nonlinearity in the range of small strains is linked to the signd : d, not appearing the strain rate tensor d in the resulting tangent stiffness formulation ford : d [ 0 nor ford : d 0. Therefore, for practical purposes, such models can be considered conceptually closer to multilinear models than to incrementally nonlinear models when describing soil behavior in the range of small strains.
These models can reproduce nonlinear, hysteretic and dependent on recent history soil behavior, characteristic of Zone II.
(C.1) Hypoplastic models Hypoplastic models represent a particular case of incrementally nonlinear models. Specifically, hypoplastic theory was initially developed by Kolymbas [312]. Subsequently, two hypoplastic formulations were distinguished: K-hypoplasticity [313] and CLoE hypoplasticity [314,315]. Afterward, these formulations were extended by Gudehus [316], Bauer [317], Wolffersdorff [318], Niemunis and Herle [151] and Herle and Kolymbas [319], among others. The first generation of hypoplastic models Using the envelopes of Gudehus in the Rendulic space, it is possible to decompose the effect of the linear and nonlinear part of the previous constitutive equation (Fig. 32).
Hypoplasticity offers a theoretical framework for describing the nonlinear mechanical behavior of granular materials, which is capable of reproducing soil strains by reordering their solid skeleton. Furthermore, unlike elastoplastic theory, in hypoplastic theory it is not necessary to neither define yield surfaces nor flow rules, since these can be derived from expression (36).
Nevertheless, classic hypoplastic models have some limitations. Among them, the requirement of a large number of parameters that limit their practical use, as well as the fact that they are not able to reproduce correctly soil behavior in the range of small strains during cyclic loadings due to ratcheting effect. Classic hypoplastic models must incorporate new state variables to describe the hysteretic and dependent on recent history soil behavior [107]. To solve this problem, diverse strategies have been proposed, some of them exposed by Tamagnini and Viggiani [66]: • Include the void ratio as a state variable to obtain a unified formulation for soils with different densities and to incorporate the concept of critical state. In this case, ÞI are considered. An expression of the tensors L r 0 ; e ð Þ and N r 0 ; e ð Þ appears in the work of Wolffersdorff [318] or Gudehus [316], including the barotropy and picnotropy factors.
• Introduce a structure tensor to consider the inherent soil anisotropy. • Introduce an internal variable dependent on strains for modelling the cyclic soil behavior. • Combine a hypoplastic model with a paraelastic one (hypoplastic hybrid models). Table 11 analyzes small strain stiffness formulation, hysteretic behavior treatment and general aspects of the hypoplastic model with intergranular strain of Niemunis and Herle [151] that considers soil behavior in the range of small strains.
According to Niemunis and Herle [151], the strain in a granular material is composed of the granular interface strain, which should be understood as a macroscopic measure of the microdeformations of the interface between particles, and the strain caused by the reordering of the grains that constitute the solid skeleton. The intergranular strain tensor d, which aims to characterize the strain of the granular interface in the range of small strains, is defined and partially solves the effect of excessive stress or strain accumulation during cycles. Other remarkable aspect is that this model allows considering different degradation curves based on the rotation of the strain path which provides versatility to the model and adaptability to experimental results. However, it requires parameters that are unusual in the engineering practice. Soil model parameters for different sands can be found in Mohammadi-Haji and Ardakani [320]. Mašín [321] proposed some modifications in this model to solve the effect of g ¼ q=p 0 on G 0 and the effect that the axis p 0 ¼ 0 has on the degradation curve when it is crossed by a stress path. Wong and Mašín [322] introduced the behavior of unsaturated soil to the Niemunis and Herle model, and recent applications can be found in Mohyla et al. [323].
Hypoplastic models have historically been barely used in professional practice due to their formulation and parameters away from the elastoplastic ''standard''. However, considerable efforts are being made to incorporate these models into commercial software. In the case of the Niemunis and Herle model, it is considered especially useful, beyond its ability to reproduce certain aspects of soil behavior, the mathematical definition of the concept of intergranular strain by means of the state variable d that takes the role of internal strain in the sense given by Gudehus [107]. The mathematical structure of d allows extending its use to elastoplastic models [324], taking advantage of its great potential, as is the case of the EPHYSS model [156].
(C.2) Hypoplastic hybrid models Some recent advanced hypoplastic models combine different concepts from previous theoretical frameworks. This section briefly shows the specific case of the combination between hypoplasticity and paraelasticity, proposed by Niemunis and Prada-Sarmiento [176]. Table 12 analyzes small strain stiffness formulation, hysteretic behavior treatment and general aspects of Niemunis and Prada-Sarmiento hypoplastic hybrid model that considers soil behavior in the range of small strains.
The Niemunis and Prada-Sarmiento model improves the hypoplastic model with integranular strain, providing it with greater robustness to reproduce soil hysteretic behavior, however, as an hybrid model, it not only inherits the strengths of the base models, but also the weaknesses, such as the unusual parameters and the lack of familiarity with these theoretical frameworks in the professional practice. Al-Tabbaa and Wood [10] Elastic behavior (f 1 \0):
It is an extension of the model of Al-Tabbaa and Wood [10], in which a new bubble f 1 ¼ 0 is considered to introduce the influence of the soil history. The expression given by Viggiani [290] is adopted for the tangent elastic shear modulus G t . The three surfaces are based on the formulation of the MCC model. Messerklinger [291] identifies some contradictions in the predictions of the model, which are resolved by McDowell and Hau [292] and Masín [293] Small Strains in Soil Constitutive Modeling 3259

Conclusions
Soil behavior in the range of small strains should be considered in the analysis of geotechnical problems, especially when it is necessary to build new structures in dense urban areas and to protect nearby existing structures, or when sensitive constructions are built or affected. Several circumstances during the 70 s of the twentieth century motivated the development of new theoretical frameworks and constitutive models that tried to explain the soil behavior in the range of small strains, characterized by being nonlinear reversible, hysteretic and dependent on recent history.
The theoretical framework of Jardine [6], based on the theories of KYS and plasticity, allows describing the behavior of soil stiffness according to the range of stress and strain to which it is subjected, giving place to four differentiated zones (Zones I to IV). This theoretical framework has been very useful to identify the strategies to reproduce soil behavior in the range of small strains used in the different constitutive models that have been analyzed.
Furthermore, in order to classify these constitutive models, the concept of tensorial zone [136,137] has been used, grouping: (1) within the incrementally linear models with a single tensorial zone, the linear and the nonlinear elastic models (Cauchy elastic models, pseudoelastic models, hyperelastic models and hypoelastic models stricto sensu); (2) within the incrementally multilinear models with several tensorial zones, the hysteretic models (paraelastic and quasi-hypoelastic models) and the advanced models (classic elastoplastic models, multisurface models, brick models, bounding surface plasticity models, bubble models and multilaminated models); and (3) within the incrementally nonlinear models with infinite tensorial zones, the hypoplastic and the hypoplastic hybrid models.
In addition to that, there have been established four elements considered fundamental for the characterization of soil hysteretic behavior and which have led to a specific model classification criterion. These elements are: (1) reversal criterion, which can be extrinsic or intrinsic; (2) memory rules, which allow storing information of one, several or all active reversal points; (3) the effect of reversals on the variables that control the degradation, which can be fully or partially reinitialized after a reversal; and (4) the effect of reversals on maximum soil stiffness, where the stiffness recovery can occur continuously or discontinuously with the rotation angle of the recent stress/ strain path, being total or partial.
Using the previous classification criteria, a total of 54 constitutive models that consider soil behavior in the range of small strains have been analyzed and classified. The study of these models has focused on describing and explaining the constitutive relations that allow the reproduction of nonlinear reversible, hysteretic and dependent on recent history behavior of the soil, as well as their applications, advantages and disadvantages. It should be noted that the fact of not having the numerical codes of most of these models complicates and limits their study, idea shared by Gudehus [107].
The most relevant aspects of the formulation of those models have been presented. All of them are capable of reproducing, to a greater or lesser extent, the behavior of the soil in the range of small strains. Some features that should be highlighted follow: i. The robust, clear and intuitive formulation of paraelastic models, based on the Hueckel and Nova model [22], to reproduce soil hysteretic behavior, in addition to its versatility to be combined with other theoretical frameworks. So far, interesting models have been developed using concepts from paraelastic theory combined with elastoplastic, bounding plasticity, multilaminate or hypoplastic models [140,153,154,156,[174][175][176][177]. ii. The ability of quasi-hypoelastic hysteretic models [67,156] to integrate multiple aspects of the soil in the range of small strains, in addition to their flexibility to be combined with plastic models, give them a great potential to be used in practical applications. However, these models also present some disadvantages [66,182,187].  iii. The balance between complexity and usability of elastoplastic models that gives them the greatest potential to be used in professional practice. Among all those analyzed in this paper, the Plaxis HS-SS model (based on the HS-S model [67]) and the EPHYSS model [156] stand out. The latter improves the HS-SS model and resolves some inconsistencies whose effects can have considerable influence on the numerical simulations of boundary value problems, as they are cumulative [156,175,[243][244][245]. iv. The potential use of brick models, based on the Simpson model [31], in the development of future constitutive models. This is due to their ability to be easily combined with classic or advanced models (such as the HS-Brick Model [155], which improves the HS-SS model), in addition to its clear theoretical framework that allows to easily understand part of the formulation of other models, as happens with the history tensors of the advanced elastoplastic models HS-SS and EPHYSS. v. The ability of multilaminate models, such as the Schädlich and Schweiger [154] model, to reproduce multiple soil characteristics, including its complex behavior in the range of small strains, using parameters known in professional practice and a clear and no contrived formulation. This place them in an advantageous position compared to other incrementally multilinear models, and although their use has been limited in practical applications, some commercial software, such as ZSoil, already incorporate basic multilaminate models in their libraries. vi. The advantages of hypoplastic models with intergranular strain, based on the Niemunis and Herle model [151], to define the soil behavior in the range of small strains through the state variable d (intergranular strain), whose mathematical structure allows its use to be extended to elastoplastic models [324], such as the EPHYSS model [156]. However, hypoplastic models have historically been barely used in professional practice due to their formulation and parameters away from the elastoplastic ''standard''.
In addition to these considerations, and regarding the behavior of the soil in the range of small strains, it is considered that the multisurface, bounding plasticity, bubbles or hypoplastic hybrid models that have been analyzed in this paper do not present great potential for their use in practical applications, due to: (1) their complex   formulations; (2) the requirement of numerous unusual parameters in professional practice; (3) the lack of description of the hysteretic and dependent on recent history soil behavior in some cases; or (4) their limited use in solving boundary value problems. Future improvements of these models, especially in the case of multisurface models [252], might allow them to be considered in practical applications. It should be mentioned, however, that these models can correctly reproduce other soil characteristics, such as irreversible strains, plastic anisotropy, plastic hardening, contractive/dilatant behavior, critical state in clays and the effect of the soil fabric. Finally, the incorporation of formulations that can reproduce soil behavior in the range of small strains within well-known models used in professional practice, would allow professional engineers to familiarize themselves with the characteristics and complexity of non-linear, hysteretic and dependent on recent history soil behavior in the range of small strains. Once this first goal has been achieved, introducing more complex models into commercial software could be less traumatic.

List of symbols
Soil parameter in the expression of Viggiani [290] used for the calculation of G t in the 3-SKH model [123] A Soil parameter in the model of Jardine et al. [207] A Soil parameter in the expressions that relate G 0 with e, OCR and p 0 A i Functions of soil parameters Þin the model of Tatsuoka and Shibuya [115] A G Soil parameter in the SDMCC model [148]  Soil parameter in the MIT-S1 model [149] related to soil compressibility after a reversal Secondary consolidation coefficient Soil parameter that controls nonlinearity during an unloading in the MIT-S1 model [  Soil parameter in the CASM-c model [143] Water pressure p 0 a Value of p 0 in the center of the f 1 of the 3-SKH model [123] p 0 b Value of p 0 in the center of the bubble f 2 of the 3-SKH model [123] p 0 c Size of the major axis of the bounding plasticity surface (consolidation stress) in the NAHOS model [144] and bubble models of Al-Tabbaa and Wood [10] and Stallebrass and Taylor [123] p 0 j Mean stress in the image point on the bounding plasticity surface in the models CASM-c [143] and NAHOS [144]  Reference mean stress in the expressions of G ap s and G t;ur in the HQH and EPHYSS models [156,175] p 0 rev Transformed variable p 0 ¼ P 0;p r 0 ð Þ on the last reversal point in the models MIT-E3 [142] and MIT-S1 [149] p 0 SSR Projection of the mean stress on the SSR p 0 a Values of p 0 in the center of the bubble f 1 in the model of Al-Tabbaa and Wood [10] p 0 0 Mean stress that fulfills G 0 p 0 ¼ p 0 Soil parameter in the model of Atkinson [32] r Soil parameter in the CASM-c model [143] r Soil parameter in the model of Ramberg and Osgood [194] r Soil parameter that controls nonlinearity during an unloading in the MIT-S1 model [149] Soil parameter in the model of Matasovic and Vucetic [210] s Soil parameter in the model of Papadimitriou et al. [140] which adopts values Parameter that represents the soil structure s 0 0 Initial value of s 0 s Ã f Reference value of s Ã for reconstituted soil Soil parameter of the HS-Brick model [155] that represents the string length of b-th brick Value of s 0 in the last reversal point in MIT-E3 model [142] and value of s in the last reversal point in MIT-S1 model [149] s R Deviatoric stress tensor in the reversal point R S 0 Surface of V S Ratio between the sizes of the bubbles f 2 and f 1 in the 3-SKH model [ Soil parameter in the model of Papadimitriou et al. [140] X X 2 X point of the initial configuration y max ¼ s max =G 0 c max y 1 Soil parameter in the model of Ramberg and Osgood [190] Soil parameter in the model of Jardine et al. [203] c Soil parameter in the model of Puzrin and Burland [152]  Value of c t in the last reversal point c 0: 7 Shear strain value in the model of Dos Santos and Correia [231] for which G ap s c 0:7 ð Þ ¼ 0:722G 0 c 1 Threshold shear strain in the model of Papadimitriou et al. [140] from which G t ¼ G 0 = 1 þ 2 1=a 1 À 1 ð Þ ð Þ (minimum value of G t ) is considered c R Shear strain at the last reversal point in the model of Pyke [158] c r Shear strain in the model of Darendeli [229] which complies G ap sĉr ð Þ ¼ 0:50G 0 c b Shear strain invariant of actual relative strain distance between the man (following Simpson's analogy) and b-th brick in the HS-Brick model [ Total strain in 1D e deg ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi e n;dev À e n;dev;centre À Á 2 þ c s À c s;centre À Á 2 þ c t À c t;centre À Á 2 q e f Axial strain at failure e i I-th principal strain, i ¼ 1; 2; 3 e n;dev Total deviatoric strain due to the normal stress to the considered plane e n;dev;centre Value of e n;dev in the last reversal point Deviatoric strain e qr Reference deviatoric shear strain equal to deviatoric shear strain in which sear modulus reduces to 70% of its maximum value e r Reference strain equivalent to E 0 =E 0 max ¼ 0:50 e s ¼ 2=3 e 1 À e 3 ð Þ e th Total threshold strain that separates linear from nonlinear behavior Initial value of e v e v;rev Value of e v in the last reversal point in the models MIT-E3 [142] and MIT-S1 [ Slope of the swelling line in the plane t À ln p 0 ð Þ j Ã Slope of the swelling line in the plane ln e ð Þ À ln p 0 ð Þ j 0 Maximum slope of the swelling line in the plane t À ln p 0 ð Þ k Slope of the primary consolidation line in the plane ln e ð Þ À ln p 0 ð Þ k Lambda coefficient of Lamé k Ã Slope of the primary consolidation line in the plane t À ln p 0 ð Þ _ k i Plastic multiplier associated to g i k m Eigenvalue of the spectral decomposition of r 0 k Si I-th eigenvalue of e in the SSOM and HS-S models [67] and in the HS-SS model implemented in Plaxis K Soil parameter in the EPHYSS model [156] l Mu coefficient of Lamé l Ã Soil parameter in the NAHOS model [144]  Scale factor that controls the shape of the degradation curve of G ap s . A value of n ¼ 1 is considered for the calculation of G ap s during primary loading and n ¼ 2 for the calculation of G ap s during unloading/reloading n 0 Variable that measures the mean stress with respect to the last reversal point in the MIT-S1 model [149] n 0 s Variable that measures the deviatoric stress with respect to the last reversal point in the MIT-S1 model [149] Maximum historic effective vertical stress Effective stress in the reversal point R r Total stress tensor r 0 ¼ r À p w 1 is the effective stress tensor r 0 Effective stress tensor of the image point on the bounding plasticity surface in the NAHOS model [144] r 0 i Effective stress tensor on the i-th plane in the multilaminated models r 0 S Elastic centre in the NAHOS model [144] r 0L Effective stress tensor in the reversal point L DL r 0 ¼ r 0 À r 0L s Shear stress s f ¼ G 0 c r s max Maximum shear stress s oct ¼ ffiffiffiffiffiffiffi ffi 1=3 p s k k s 1 ¼ y 1 G 0 =c max t ¼ 1 þ e # 0 i i ¼ 1; 2 soil parameters in the model of Tatsuoka and Shibuya [115] u 0 Maximum effective friction angle u X; t ð Þ Position of material point X 2 X at time t v Parameter of the model of Niemunis and Herle [151] Soil parameter in the model of Niemunis et al. [153,174]  Soil parameter in the models MIT-E3 [142] and MIT-S1 [149] x s Soil parameter that controls nonlinearity during the application of shear stress in the MIT-S1 model [149] X Domain on which the problem to be solved is defined X Second order symmetric tensor of constant coefficients that complies that C 0 Á X is semi-defined positive [22] #Z Number of tensorial zones   [123] 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/.