Development of a new advanced elastoplastic constitutive model that considers soil behavior at small strains. The EPHYSS model

The Elastoplastic Hysteretic Small Strain (EPHYSS) model is an advanced elastoplastic model as a result from the combination of the Hysteretic Quasi‐Hypoelastic (HQH) model that considers strain‐induced anisotropy and can reproduce the nonlinear reversible, hysteretic and dependent on recent history soil behavior, and the Cap‐Cone Hardening Soil Modified (HSMOD) model that can reproduce soil plastic behavior. EPHYSS model uses state variables that define different short and long‐term memory levels which provide it with robustness for the reproduction of soil hysteretic behavior and confers it a great versatility and adaptability to experimental results. It also corrects some inconsistencies of the Hardening Soil with Small Strain Stiffness (HS‐SS) model of Plaxis whose effects can have a considerable influence on the numerical simulations of boundary problems. The performance of EPHYSS and a comparison with the HS‐SS model is presented in some experimental and numerical tests, and in a boundary value problem of a large excavation in Barcelona.

of Jardine, in which this behavior is not reversible, either in a linear or nonlinear way, due to the appearance of plastic strains. 1 There are numerous advanced constitutive models capable of simulating the behavior of the soil in the whole strain range, 2 although the use of many of them is reduced to an academic use. However, some of these advanced models have managed to extend to the professional practice. Among them, stands out the well-known Hardening Soil with Small Strain Stiffness (HS-SS) model, based on the work of Benz, 3 implemented in Plaxis and recently in other numerical software.
After analyzing in depth 54 constitutive models that consider the soil behavior in the range of small strains, there have been identified some aspects that can be improved in the elastic part of the HS-SS model in relation to: (1) nonlinear behavior; (2) hysteretic behavior; (3) deviatoric strains reversals effect on the elastic bulk modulus; (4) consideration of the strain induced anisotropy; and (5) correction of some inconsistencies that have been detected. 2,[4][5][6][7] All these points have motivated the development of the Elastoplastic Hysteretic Small Strain (EPHYSS) model with the aim of being used in the geotechnical engineering professional practice. The entire code of the EPHYSS model, implemented in the User Defined Soil Model (UDSM) modulus of Plaxis may be consulted in Castellón. 4

THEORETICAL FRAMEWORK OF THE EPHYSS MODEL
The EPHYSS model has multiple tensorial zones 8,9 and, therefore, it is an incrementally multilinear model, 10 although its indirect stiffness dependence on deviatoric strain increment when describing soil behavior in the range of small strains allows to consider it, in practice, as an incrementally nonlinear model. Specifically, the EPHYSS model belongs to the subtype of incrementally multilinear models called as advanced models of the elastoplastic type following the classification of Castellón and Ledesma. 2 Such model results from the combination of the Hysteretic Quasi-Hypoelastic (HQH) model, which aims to reproduce the quasi-static reversible behavior of the soil in the range of small strains (Zones I and II of Jardine), and the Hardening Soil Modified (HS MOD ) model, which aims to reproduce the irreversible behavior of the soil in the range of intermediate and large strains (Zones III and IV of Jardine). 4 The HQH model is a quasi-hypoelastic hysteretic model of the variable moduli type. 2,[11][12][13][14] The Hysteretic denomination is due to the capacity of its structure to partially comply the Generalized Masing Rules, 15,16 while the Quasi-Hypoelastic denomination is due to the use of nonlinear apparent stiffness moduli in its constitutive equation.
The reversible behavior of the soil in the HQH model is conditioned to the closure of the strain cycles, in a similar way to what happens in paraelastic models. 2 To define the cycles, HQH model requires information about soil state at the extreme points of those cycles, which constitutes the active reversal points. An active reversal point that belongs to a specific cycle has an effect on the soil stiffness while the cycle remains opened. Once the cycle (or larger cycles that contain it) is closed, the active reversal point disappear and it is erased from HQH state variables. It should be added that the EPHYSS model allows reproducing the behavior of the soil along paths that include common strain cycles in quasi-static geotechnical problems, but it cannot reproduce soil dynamic behavior, since it does not consider factors such as the accumulation of irreversible strain after multiple cycles, cyclic hardening/softening or elastoplastic coupling. The HS MOD model belongs to the models denominated Cap-Cone and has four yield surfaces. This model considers two different isotropic hardening mechanisms: a shear hardening mechanism on the Cone-type surface, which allows to reproduce the plastic strains in hard soils generated by deviatoric loadings, and a compression hardening mechanism on the Cap-type surface, which allows to reproduce the plastic strains generated by oedometric and isotropic loadings. EPHYSS model presents some important limitations that are inherited from the plastic model HS MOD and are described in Section 6. Therefore, it is recommended to not use EPHYSS model for numerical analyses of geotechnical works in which intermediate or large strains play a significant role.

Strain domains of the model
Two strain domains are considered. Domain 1 corresponds to values Δ oct ≤ ur and domain 2 corresponds to values Δ oct > ur , where refers to the last reversal point that conforms the endpoint of the active strain cycle and refers to the deviatoric strain rotation angle corresponding to the active degradation curve.
In domain 1 the reversible component of the strain (whose reversibility is conditioned to the closing of the open strain cycles) follows a behavior that is nonlinear with both deviatoric strain and mean stress, and strain-induced anisotropy is also considered, while in domain 2 the reversible component of the strain (whose reversibility is also conditioned to the closing of the open strain cycles) follows a behavior that is both linear with the deviatoric strain and nonlinear with the mean stress.
The model considers the hysteretic behavior in both domains thanks to its ability to reproduce cycles in which shear strain produces a degradation of the elastic shear stiffness (which happens in domain 1) until it reaches a minimum value (which happens in domain 2). It should be clarified that both domains can appear during primary loading, but it is meaningless to talk about hysteretic behavior until a reversal takes place and at least one cycle appears.
In addition to that, the model considers the recent deviatoric strain history in both domains. It has to be understood as recent history those deviatoric strain paths that modify the model state variables which control elastic shear stiffness. Once strain cycles are closed, these state variables are reinitialized and the recent history erased.
In order to reproduce the soil hysteretic behavior (partially complying with the Generalized Masing Rules), as well as to consider the recent deviatoric strain history, the model uses a set of state variables whose formulation will be introduced in Section 3.5. These variables are the following:

General and incremental constitutive equation in domain 1
The general constitutive equation in domain 1 ( Δ oct ≤ ur ) is the following: The elastic secant bulk modulus ′ in domain 1 adopts the following expression based on a generalization of the expression of Duncan et al., 17 which is widely accepted: The apparent secant shear modulus ap adopts the Expression (4) For the calculation of ap , the structure of the HQH model allows directly applying the Hashiguchi strategy, 19 which consists of considering a value = 1 in primary loading branch and = 2 in unloading or reloading branches. However, that, as well as , depends on the total recent deviatoric strain history, giving place to an area of possible values of ap ( Figure 1B), so that the HQH model considers infinite degradation curves of soil shear stiffness as a function of the rotation in the recent deviatoric strain path, which offers great versatility and adaptability to experimental data. The incremental constitutive equation in domain 1 can be deduced from the Expression (1).

Incremental constitutive equation in domain 2
The incremental constitutive equation in domain 2 ( Δ oct > ur ) is the following: The expression of the elastic tangent bulk modulus ′ in domain 2 is the following: Regarding to the (hypo)elastic tangent shear modulus in domain 2, a minimum value ,ur is considered from the limit strain between both domains ( Δ oct = ur ): The fact that ur coincides with the maximum value of oct allows applying the expression of ap in domain 1 to deduce the value of ur = ur ( ,ur Based on the foregoing, the value of ur is obtained by substituting ap = ,ur and oct = ur , in the Expression (16).

Poisson's ratio
Drained Poisson's ratio will be variable in both strain domains, since it depends on the value of the stiffness moduli according to ′ = ′ ( ′ , ap ) in domain 1 and ′ = ′ ( ′ , ,ur ) in domain 2.
The value of ap in the range of small strains corresponding to domain 1 can lead to values of ′ too small or even negative and, therefore, physically unreal. To avoid this, it is possible to limit the value of the bulk modulus ′corr through forcing drained Poisson's ratio to be maintained above a threshold value ′ min when deviatoric loadings appear. The same type of correction can be applied in domain 2 to ′corr . It is assumed ′ ,min ≈ ′ ,min ≈ ′ min .

HQH model state variables
After a reversal and before a later monotonous strain, the soil behavior suggests a gradual adaptation of its internal state, characterized by hidden state variables, until it depends exclusively on the stress tensor and the void ratio, 20 state in which proportional strain paths lead to proportional stress paths and, therefore, soil internal state is located within the Swept Out Memory (SOM) region. 21  , , , , MEM , MEM and MEM ) to describe the nonlinear reversible, hysteretic and dependent of the recent history soil behavior in Zone II of Jardine, which includes domain 1 and part of domain 2, and one additional state variable ( ap ,min, =1 ) to simulate initial stiffness and stiffness degradation during primary loading when HQH and HS MOD models are combined.

State variable ℎ
The tensor is a history tensor that stores the value of recent total deviatoric strains, understanding as recent those strains that have not been erased yet by a reversal. The state variable acts as a short-term memory variable that can be totally or partially reinitialized after a reversal. The proposed formulation for the history tensor ratėis based on that of the intergranular strain tensor ratė, 22 although significant modifications are introduced.
where cos ( ) =ˆ∶ˆ̇, = oct ∕ ur and is a numerical parameter that controls the rate with whicḣevolves (values ≥ 100 are adopted to havė≈̇in degradation processes after a 90 • ≤ ≤ 180 • deviatoric strain rotation, complying → 0 if < 1 and → 1 if ≈ 1). When the closing of an active strain cycle takes place between the active reversal points − 1 and ( ≥ 2 when a strain cycle exists), which happens when Expression (21) is fulfilled for = 1, 2, the history tensor is updated by modifying its modulus but maintaining its previous direction as indicated in Expression (22). This allows recovering the value of the variable ,( +1) oct corresponding to the cycle that is immediately superior to the one that is being closed, with a reasonable computational cost. acts as a short-term memory variable that can be totally or partially reinitialized after a reversal.
where cos ( * ) = (ˆ∶ˆ̇) * . The value of the maximum shear modulus after stress reversals in different soils has been measured in many studies. [23][24][25][26][27] According to these studies it can be concluded that the ratio between the maximum shear modulus after a 180 • stress rotation ( ) and the maximum shear modulus after a 90 • stress rotation ( that considers stress (nor strain) reversals). Nevertheless, in Niemunis and Herle model the reversals depend on the total strain rotations, thus giving place to a coupling between the variation of volumetric strains and the variation of the shear stiffness, 22,28 but it is possible to show that ifˆΔ ∶ˆ= 0 holds, thenˆΔ ∶ˆ̇= 0, provided that ocṫoct = 0, so, in these cases, for the estimation of from ∕ it will be assumed that the condition ocṫoct = 0 is satisfied, that is, the soil does not experience volumetric strains or these strains are small before or immediately after the 90 • strain rotation. In Table 1  The state variable , which stores the value of the total deviatoric strain tensor at the last reversal point that conforms the endpoint of the active strain cycle, appears, on the one hand, in the expression Δ oct = √ 4∕3 (‖ ‖ − ‖ ‖), which allows to know if the incremental constitutive equation that must be applied is the one corresponding to domain 1 or domain 2, and, on the other hand, in the expression Δ = − , with which it is possible to calculate the productˆΔ ∶ Δt hat identifies the reversals in which the state variables , , , MEM , MEM and MEM must be stored. Conversely, the state variable , , which stores the value of the elastic deviatoric strain tensor at the last reversal point that conforms the endpoint of the active strain cycle, appears in the deviatoric part of the elastic incremental constitutive equation of domain 1,̇= 2 aṗ+ 2̇a p ( − , ). Both state variables act as short-term memory variables that are totally reinitialized after a reversal, and their update occurs according to the following scheme, where ( ) is the number of reversal points that are active at the beginning of the calculation step ( ) → ( + 1): • Whenˆ( ) ∶ Δˆ( +1) ≤ cos( * ,( ) ) andˆΔ ,( ) ∶ Δˆ( +1) ≤ cos( * ,( ) ), a deviatoric strain rotation that generates a new reversal point that will conform the endpoint of a new strain cycle occurs: • When , with = 1, 2, the closing of the strain cycle between the active reversal points ( ) − 1 and ( ) ( ( ) > 2) occurs: when ( ) = 2, ,( +1) = 0 and , ,( +1) = 0.
• When neither the deviatoric strain rotation that generates a new reversal point nor the closing of any strain cycle occur: The first cycle to be closed will be the one corresponding to the smallest active strain cycle. The extremes of such cycle correspond to the last two values of ‖ ‖, ‖ ‖ and ‖ ‖, memorized, respectively, in the state variables MEM , MEM and MEM . Such values are those that have been identified with the ( ) − 1 and ( ) components of such variables in the calculation step ( ) → ( + 1). Moreover, as can be seen, when a cycle closure occurs, the state variables and , are updated by modifying their modulus but maintaining their previous direction. This update allows partially recovering the characteristics of and , corresponding to the cycle that is immediately superior to the one that is being closed, with a reasonable computational cost. • Whenˆ( ) ∶ Δˆ( +1) ≤ cos( * ,( ) ) andˆΔ ,( ) ∶ Δˆ( +1) ≤ cos( * ,( ) ), a deviatoric strain rotation that generates a new reversal point that will conform the endpoint of a new strain cycle occurs (in this case ( +1) = ( ) + 1): • When , with = 1, 2, the closure of the strain cycle between the active reversal points ( ) − 1 and ( ) occurs (in this case ( +1) = ( ) − 2 and the ( ) − th components of these variables, that comply ( ) > ( ) − 1 and correspond to the reversal points that conform the endpoints of the active strain cycles imbricated in the cycle that is being closed, will be erased): • When neither the deviatoric strain rotation that generates a new reversal point nor the closure of any strain cycle occur (in this case ( +1) = ( ) ): The condition = 0 is necessary for the soil state to be located within the SOM region.

State variable
ap ,min, =1 The HQH model does not directly apply the criterion of Hashiguchi 19 for the calculation of ap , despite having an internal structure that allows it. To simulate the initial stiffness of the soil and its degradation in the range of small strains during primary loading, the hardening laws of the plastic model with which the HQH model is combined are modified through the factor ℎ = ℎ ( ap ,min, =1 ). The state variable ap ,min, =1 is updated as follows in each calculation step:

Hysteretic behavior
The HQH model considers the hysteretic behavior in the deviatoric strain component and not in the volumetric one, because the variation of the bulk modulus between elastic isotropic unloadings and reloadings is negligible. 34 This distinguishes the HQH model, provided that ′ > ′ min is met, from several models that consider the dependence ′ = ′ ( , ′ ) or = ( ′ , ′ ), in which the hysteretic behavior that controls induces an hysteretic behavior on ′ or viceversa.
A model that consider the hysteretic behavior of the soil must define the following concepts: (1) reversal criterion; (2) memory rules; (3) reversal effects on the degradation variables; and (4) reversals effect on the maximum soil stiffness. 4 3.6.1 Reversal criterion The reversal criterion of the HQH model is intrinsic, which means that it arises from the own model equations. It is formulated using strains and the reversals take place in a continuous way with the rotation angle of the recent deviatoric strain patĥ∶ˆ̇. The HQH model distinguishes between: • Reversals that affect the value ap through the modification of the state variables

Memory rules
HQH model can memorize information in multiple reversal points and, combined with a plastic model, partially satisfies the Generalized Masing Rules 15-16 to reproduce the hysteretic soil behavior.
• Rule Nr. 1: To be able to reproduce the initial soil stiffness and its degradation in the small strain range during primary loading branch, the hardening rules of the plastic model with which the HQH model is combined are modified through the factor ℎ (see Section 4.3). and all the information of the total and elastic deviatoric strain tensors in the last reversal point that conforms the endpoint of the active strain cycle through the state variables and , , respectively. As well, HQH is capable to store partial information about the strains in all the reversal points that configure endpoints of active strain cycles through the state variables MEM , MEM and MEM .

Reversals effect on the maximum soil stiffness
The reversals effect on the maximum soil stiffness prior to such reversal ( Figure 3).

Stability criterion
To meet the stability criterion of Hill, 35  , so that the value of ap cannot be reduced after a reversal.

EPHYSS MODEL AS A COMBINATION OF HQH AND HS MOD MODELS
The EPHYSS model is an advanced elastoplastic model that considers the soil behavior in the range of small strains.
The elastic part of the EPHYSS model is described with the HQH model, so that: On the other hand, the plastic part of the EPHYSS model is described with the HS MOD model, whose theoretical fundament is exposed in this section.

4.1
Yield surfaces of the EPHYSS model

Cone yield surface and Mohr-Coulomb limit surface mc
The EPHYSS model considers a Cone yield surface that controls the plastic strains of the soil under deviatoric loading. Such surface is formulated in terms of axial strain as = 2 1 − 2 1 − , being the plastic state variable ≈ 2 1 under the hypothesis of hard soils in which = 3 oct ≈ 0. To define the Cone yield surface, it is enough to do it in the stress space sectors corresponding to The EPHYSS model does not explicitly use the ′ref 50 modulus in the formulation of the Cone yield surface, but uses the initial stiffness modulus ′ref , which is an internal parameter of the model that has to be calculated through an internal algorithm in which it is necessary to introduce the input parameter The yield surface can harden until the plastic state variable reaches the limit value determined by the resistance criterion of Mohr-Coulomb, which happens when = .
Facing the numeric implementation of the EPHYSS model, it is interesting to express the Mohr-Coulomb limit surface mc in the same stress space sectors used for the Cone yield surface :

Cap yield surface
The EPHYSS model considers a third yield surface denominated Cap yield surface, which allows to close the elastic region in the direction of the hydrostatic axis and reproduce the plastic volumetric soil behavior. In Expression (48) ′ is an internal parameter of the model that has to be calculated through an internal algorithm in which it is necessary to introduce the input parameter NC 0 (see Section 4.5) and is a plastic state variable that controls the size of .

Tension Cut-off yield surface
Finally, the EPHYSS model considers a fourth yield surface, denominated Tension Cut-off yield surface, which considers the soil tensile strength (− ′ 3 = ′ trac ≥ 0). To define the Tension Cut-off yield surface, it is sufficient to do it in the principal stress space sectors corresponding to

Plastic potentials and flow rules
Dilatancy in critical state models is zero when the soil reaches the critical state. EPHYSS model is not a critical state model but can aproximate critical state phenomenon by using a "Dilatancy cut-off" criterion that consists in forcing the dilatancy to be zero when the void ratio reaches a maximum value (that is, = 0 when e = e max ). This criterion adds a new parameter e max to the model. In the model version presented in this paper, the "Dilatancy Cut-Off" criterion has not been considered.
The definition of the flow rule for the Cone yield surface and the Mohr-Coulomb yield surface is carried out through the following expressions of the plastic potentials in the principal stress space sectors corresponding to Finally, regarding to the Cap yield surface and the Tension Cut-off yield surface, associated plasticity is assumed, considering, therefore, = and ,jkl = ,jkl .

Hardening laws
The EPHYSS model uses two hardening variables: (1) the plastic deviatoric strain that controls the size of the Cone yield surface ; and (2) the preconsolidation stress that controls the size of the Cap yield surface . The hardening laws of the EPHYSS model are the following:̇= The parameter ′ ref in the Expression (59)  and vanishes when ap ,min, =1 = ,ur (see Section 3.5.5). The evolution of ap ,min, =1 follows the degradation curve corresponding to the primary loading branch, whose values are lower than that of the elastic unloading or reloading branches. For that, in the calculation of ℎ the value = 1 is used following the Hashiguchi criterion. 19 Moreover, the plastic modulus depends on ′ref = ′ref ∕(1 − 1 ) according to the Expression (60), where ′ref is an input parameter of the HQH model, and ′ ,ur ∕ ′ is an internal parameter that has to be calculated through an internal algorithm in which it is necessary to introduce the input parameter ′ref oed (see Section 4.5).

EPHYSS model parameter identification
The parameters of the EPHYSS model are listed and described in can be obtained from the results of tests with internal strain measurement in which it is possible to control the rotation of the deviatoric strains, such as, for example, biaxial, true triaxial or hollow cylinder with torsion tests (since these tests are very rare in the professional practice, it has been proposed the following expression

EPHYSS model internal parameters
The EPHYSS model considers three internal parameters ′ref , ′ and ′ ,ur ∕ ′ that explicitly appear in the model formulation. These internal parameters are calculated from the three model input parameters The algorithm that has been implemented to calculate such internal parameters 4,41 considers two numerical tests: first, a drained triaxial test is reproduced to obtain the value of ′ref ; and, subsequently, an oedometric test is reproduced to obtain the value of ′ and ′ ,ur ∕ ′ .

LOCAL INTEGRATION OF THE CONSTITUTIVE EQUATIONS
The standard calculation strategy in elastoplastic models is used for the local integration of the constitutive equations of the EPHYSS model. This strategy uses the trial stress ′(tr) (63), which assumes that the soil behavior is totally elastic. The Expression (63) will be the one used when Δ ,( ) oct > ,( ) ur , that is, in domain 2. However, such expression will be modified when Δ ,( ) oct ≤ ,( ) ur , that is, in domain 1, where ′(tr,NL) is used (64).
After calculating the trial stress ′(tr) or ′(tr,NL) according to the case, and defining the elastic domain as  = { ′ | ( ′ , pl ) < 0, = 1⋯ }, the procedure is as follows: • If ′ ( ) ∈  or ′ ( , ) ∈  according to the case, the following variables will be updated: (1)  will be adopted and the following variables will be updated: (1) the stress tensor ′ ( +1) , according to Expression (67) or Expression (68) as appropriate; (2)  ∕ ′ (67) ) ∕ ′ The so-called Implicit Closest Point Projection Algorithm has been used for the local integration of the equations. This algorithm uses a Backward Euler Elastic Predictor/Return Mapping type integration scheme within the methods of radial return, and belongs to the algorithms class denominated Generalized Midpoint Algorithms within the linear multistep methods. 3,[42][43][44] Furthermore, for the return on the intersection of 2 or 3 yield surfaces, Koiter's Rules have been adopted. 45 As well, iterative algorithms within the elastic part have been incorporated to the general algorithm.
A major problem in multisurface elastoplastic calculations is to stablish a strategy to determinate on which surface or surfaces must be the Return Mapping done, especially when the stress state is next to an intersection zone between several yield surfaces. The strategy for the surface selection considered in the EPHYSS model is based on the proposal of Bonnier. 46

DIFFERENCES AND COMMON ASSUMPTIONS BETWEEN EPHYSS AND HS-SS MODELS
Differences and common assumptions between EPHYSS and HS-SS models are discussed in Castellón 4 and Castellón and Ledesma, 2 and will be summarized in this section.
The HS-SS model implemented in Plaxis is one of the few models widely used in professional practice that considers soil behavior in the range of small strains and uses known and relatively easy to obtain/estimate parameters. The HS-SS model is based on the Hardening Soil Small (HS-S) model 3 which adds, to the elastic formulation of the Small Strain Overlay Model (SSOM), 3 the plastic formulation of the Hardening Soil (HS) model 47,48 adding two modifications: (1) it replaces the dilatancy criterion of Rowe 38 with that of Li and Dafalias 36 to describe the contractive behavior of the soil; and (2) modifies the hardening laws by introducing the factor ℎ . The HS-S model significantly improves the approach to soil behavior provided by the HS model, although it has certain limitations. 3 Some applications can be found in Foster et al., 49 Ramos et al. 50 and Ledesma and Alonso. 51 The HS-SS model presents some small differences with respect to the HS-S model. 2 (5) can be combined with multiple plastic models leading to incrementally multilinear advanced elastoplastic models. When both models are used by themselves, they use the Hashiguchi 19 criterion to differentiate primary load from unloadings/reloadings. However, since the intention is to compare the EPHYSS and HS-SS models, while minimizing the differences caused by their respective plastic formulations, it has been adopted, in the EPHYSS model, the same strategy used in the HS-SS model, which consists of considering, on the one hand, a factor = 2 in the expression of the apparent secant shear modulus ap and, on the other hand, the factor ℎ that modifies the hardening laws to avoid overlapping the mechanisms of SSOM and HQH models and the mechanisms of plastic models with which they are combined, that try to explain the reduction of soil stiffness during the primary loading. Despite its great advantages, the elastic part of HS-SS model (SSOM model) presents some limitations that the elastic part of the EPHYSS model (HQH model) solves, such as: (1) the use of an intrinsic reversal criterion which means that it arises from the own model equations (and not an arbitrary extrinsic one like in the SSOM model, where the reversals depend on the product sign of the eigenvalues oḟwith the components of − in the same direction); (2) the ability to recover the stiffness in a continuous way with the rotation angle of deviatoric strain recent pathsˆ∶ˆ̇, through the use of a history tensor inspired by the intergranular strain tensor of the Niemunis and Herle model 22  (2) EPHYSS model considers a void ratio independent formulation that does not include the concept of critical void ratio, so the evolution of plastic volumetric strain at intermediate or large strains is not properly described; (3) the model is not suited for modeling materials at different void ratios with a single set of parameters; (4) kinematic hardening is not considered; and (5) the cone yield surface is formulated under the hypothesis of hard soils ( ≈ 0) and, therefore, EPHYSS is not appropriate for reproducing deviatoric stress paths in soft soils if plasticity is predominant.
The elastic part of HS-SS and EPHYSS models (SSOM and HQH models respectively) requires the same or equivalent few and simple parameters, except for an additional parameter required by the HQH model, which is the maximum shear modulus after a reversal of the deviatoric strains of 90 • (  Finally, regarding the parameters corresponding to the plastic part of the EPHYSS model, both internal and external are the same than those used in the HS-SS model. In addition to that, in EPHYSS, HS, HS-S or HS-SS models, the plastic modulus depends on ′ref = ′ref ∕(1 − 1 ), where ′ref is an input parameter of the EPHYSS model but it has to be calculated from the elastic relation ′ref = ′ref ( ′ref ,ur , ′ ur ) in the HS, HS-S or HS-SS models.

TESTS FOR THE VALIDATION AND VERIFICATION OF THE EPHYSS MODEL AND ITS COMPARISON WITH THE HS-SS MODEL
Different triaxial, oedometric and biaxial tests, compiled from diverse thesis and papers, have been simulated to carry out a partial verification of the EPHYSS model, as well as a validation and a comparison with the HS-SS model. The simulations have been run with PLAXIS 2D 2015. Table 3 shows EPHYSS and HS-SS parameters for the different soils used in the simulations. The parameters of the EPHYSS model for each material have been obtained from the HS-SS model parameters, that were extracted from the indicated references, using the expressions in Section 6.

Triaxial tests
Numerical simulations of drained and undrained triaxial tests on loose and dense Hostun sand with different confinement stresses have been conducted, as well as drained triaxial tests on reconstituted kaolinite clay with different confinement stresses (Figure 4-7) and triaxial tests with stress rotations in the deviatoric plane on London clay (Figure 8).
In general, the approximation of both models to the experimental measurements of the triaxial tests is quite good. However, the dilatant behavior of the dense sand with − ′ 3 = 100 kPa (Figura 4C) is not correctly reproduced, and the volumetric stiffness in both loose and dense sands are bigger that those obtained in the tests ( Figures 4B and 5B), which is due to the dependence of ′ with and ′ (either in the entire range of strains in the HS-SS model, or when ′ ≤ ′ min in the EPHYSS model). EPHYSS model can improve the volumetric stiffness approximation to the experimental curves by reducing ′ min , although this would lead to important errors in the prediction of horizontal displacements in boundary value problems with significant deviatoric loadings and, therefore, it is no recommended.
In standard traixial tests shown in Figures 4-7 no strain rotations are applied and, therefore, EPHYSS and HS-SS models predict similar soil behavior as both: (1) use the same value of the maximum shear stiffness ( In addition to that, in Figures 4-6 certain differences are observed between the degradation curves of ap provided by the EPHYSS model and the theoretical curve ap according to the Expression (4). These differences are due to: (1) the effect that the factor ℎ introduces on the degradation curves in presence of plastic strains; and (2) the fact that the value of ap is affected by the condition ap = ,ur in both models when oct ≥ ur or HIST ≥ , respectively, which is not considered in the Expression (4). In Figure 7, both the EPHYSS and the HS-SS models present some limitations in the approximation of the experimental curves q−(− 1 ) of undrained triaxial tests on dense Hostun sand. This is because the adopted dilatation formulation in both models does not consider the void ratio of the material as a state variable.
In Figure 8 four samples of reconstituted London clay are consolidated to point and subsequently taken to point (OCR = 2). Once at point , a stress path is applied to each of the samples following different angles ∕ ′ in − (− ′ ) space. It is necessary to remember that in this type of tests the deviatoric strain path rotation angle is = 180 • in all paths, except in the BOX one, where = 0 • . Therefore, the results predicted by EPHYSS and HS-SS models are very similar, which is logical since the only active degradation curve in the EPHYSS model is the one corresponding to the maximum shear stiffness (

Oedometric tests
Simulations of oedometric tests have been conducted on dense and loose Hostun sand with four load cycles (Figure 9). The approximation of both models results to experimental measurements is good. All the reversals in these tests have a

Biaxial tests
Numerical simulations of drained biaxial tests on loose and dense Hostun sand with different confinement stresses have been conducted (Figure 10). The approximation of both model results to the experimental measurements is good, except in the range of strains in which a process of localization and shear band formation appears. Since no strain rotations take place in these tests, the EPHYSS and HS-SS model predict similar results. The slight differences between them are the same that have been explained in Section 7.1 for triaxial tests. Furthermore, numerical simulations of biaxial tests with rotations in the recent strain path have been performed on Hochstetten and Ticino sand ( Figure 11). Following the scheme of Figure 11A Figure 11B shows the degradation curves of the apparent tangent shear modulus ap in the simulations conducted with the following models on Hochstetten sand: (1) hypoplastic with intergranular strain according to Benz 3,22 ; (2) Figure 11B). This is due to the fact that the value of the parameter corresponding to the maximum shear modulus for strain rotations of 90 • is common in both models, but this parameter does not exist in the HS-SS model, which reproduces the increase of the stiffness corresponding to such rotation only by the reduction of the state variable HIST . As well, the jump that can be seen in EPHYSS results for = = 90 • is due to the fact that belonging to one domain or another is evaluated by the strain Δ oct , while the stiffness depends on the state variable oct . Finally, a biaxial test with strain reversal on Ticino sand has been carried out (Figure 11 produce the same soil stiffness recovery that appears in the other three models. This is due to the high value of the strain − xx = − yy = Π before the strain rotation is produced, which, in turn, gives place to high levels of the HS-SS history tensor components ij before the reversal and, therefore, HIST ≫ . In these cases, after the strain rotation, the yy component of the tensor is reinitialized, but not the xx component, which keeps the previous value of such rotation and is high enough to provide values of HIST which still are superior than , therefore, small or no stiffness recovery is produced. The latter does not occur in any of the other three models (hypoplastic with intergranular strain, multilaminated for small strains or EPHYSS), in which the elastic stiffness recovery is less dependent on the value of the accumulated strain Π due to the consideration of more complex history tensors or other state variables ( MEM , MEM and MEM in the case of the EPHYSS model) that allows memorizing more information of the recent history apart from that stored in the history tensor. And additional advantage of the EPHYSS model is the use of very similar parameters to those of the HS-SS model, which facilitates its use in the geotechnical professional practice.

Numerical tests
A group of numerical simulations of oedometric ( Figure 12) and triaxial ( Figure 13) tests have been conducted with the aim of demonstrating that the formulation of the EPHYSS model allows to correct the inconsistencies detected in the HS-SS model. [4][5][6][7] These simulations have been run with both PLAXIS 2D 2015 and PLAXIS 2D 2018. Parameters corresponding to Hostun loose sand (Table 3)   In oedometric tests, the HS-SS model shows the same inconsistencies related to variations in stiffness during consolidation phases, nil phases and phases with small unloading/reloading, which give place to deviations in the oedometric curves, both in the elastoplastic and in the elastic branch ( Figure 12A). In all these cases, the reinitialization of all the components of the history tensor has been detected. These effects can have a very relevant influence on the results of numerical simulations, since they are cumulative.
Regarding to the deviatoric phases of the numeric triaxial tests conducted using the HS-SS model, no stiffness variations are observed when nil phases are introduced during the primary loading ( Figure 13A). Stiffness variations are observed, due to the reinitialization of all the components of the history tensor , when a nil phase is introduced during the unloading branch ( Figure 13A) and when a small unloading/reloading is introduced both in the primary loading and in the unloading branch ( Figure 13A).
Furthermore, in the EPHYSS model no variations of the stiffness that give place to the deviation of the oedometric or triaxial curves are observed (Figures 12B and 13B). This is because EPHYSS model, despite experiencing reversals similar to those observed in the HS-SS model, is capable of correcting its effect in the subsequent calculation phases thanks to the introduction of the state variables MEM , MEM and MEM which accumulate enough information from the recent history TA B L E 4 Geotechnical units 53,54  Blue-gray marly clays Pl2 of the soil to recover the historical stiffness corresponding to a specific branch of previous loading/unloading/reloading, despite the possible cycles embedded in it. This implies an improvement of the compliance level of the Generalized Masing Rules Nr. 3 and Nr. 4.

Study of a large urban excavation in Barcelona
An example of application involving a large excavation in an urban environment is presented. It refers to the future railway station of La Sagrera in Barcelona, which will constitute an important intermodal hub in the local, national and European public transport network. The construction of the station involved an excavation 650 m long, 35 to 80 m wide and 20 m deep. The geology of the area was described in the construction project. 55 Additional geotechnical site investigations were carried out when defining the groundwater drainage of the excavation. 56 The site is located in the Barcelona plain, which consist of Quaternary formations that overlie a substrate mainly formed by Paleozoic and Pliocene series. Table 4 shows the geotechnical units considered from the available information from top to bottom.
For the hydrogeological characterization of the subsoil of the future station, the transmissivities of the different layers of the soil were calibrated from: (1) the information available in the historical records; (2) the information of the different geotechnical site investigations carried out during and after the project phase; and (3) the results obtained in the analysis of the pumping tests performed in March 2011. 56 Given the nature of the soil concerned, anisotropy of permeability ( ) has been considered by adopting values of vertical = 0.1 horizontal in all materials.
The magnitude of the excavation has required a large site investigation. Multiple field and laboratory tests were conducted, including: soil identification, pressuremeter tests, SPT, pumping tests, CD, CU and UU triaxial tests, direct shear tests, oedometric tests and chemical analysis of the soil and the groundwater. Furthermore, within this investigation 4 high-quality block samples were taken from the excavation bottom in order to conduct resonant column tests using different confinement stresses (100 kPa, 200 kPa and 300 kPa) to determine soil parameters in the range of small strains. 4 Table 5 provides soils specific weight, permeabilities and parameter values of the EPHYSS and HS-SS models for each of the geotechnical units. All the materials are normally consolidated, so that OCR = 1.
As indicated in Section 6, it is important to highlight the fact that ′ and ap vary, respectively, with (− ′ ) 1 Figure 14 shows in the work plant the situation of the transversal cross-sections that contain each of the extensometers and the longitudinal cross-section that includes all of them.
Based on the geometry of the excavation and the selected phases, a total of five 2D numerical models have been made, corresponding, respectively, to the four transversal cross-sections and the longitudinal one, which refer to the profiles shown in Figure 14. The geometry of the models analyzed, in each of the calculation phases, is shown in Figures 15 and 16. Earth heaps during the work have been simulated as equivalent loadings. The water table is at the reference level + 2,5 m, and it was not reached in any of the analyzed excavation phases.
Two types of analysis have been performed: (1) drained analysis which considers drained conditions for materials corresponding to all geotechnical units; (2) undrained-consolidated analysis which considers: (a) in the materials corresponding to the less permeable geotechnical units (Ra, Qa1-Qa2, PQ2 and Pl2) and in each excavation phase, a first subphase with undrained conditions, followed by a second consolidation subphase, whose duration corresponds to the period between the respective phases of excavation in which water pressures have been almost dissipated (after the application of the undrained loading in Phase D, no subsequent consolidation phase is applied); and (b) in the materials corresponding to the most permeable geotechnical units (Qa3, PQ1 and Pl1), drained conditions have been considered in the two subphases of each calculation phase described in point (a). In the undrained loading subphases it has been considered ≈ 1 and a maximum allowed suction of = 100 kN∕ 2 in the materials located under the water table and ≈ 0 in materials located above the water table in order to simplify the problem. It should be noted that these two types of analysis will generally provide different ground deformation profiles because different effective stress paths are followed in each of them. However, in the range of small strain, the state of the soil is far from yield surfaces and the resulting profiles in both analyses will be closer as differences will only come from elastic non-linearities, which are similar in both models if no big reversals take place. Deformation profiles would be identical if linear elasticity was considered within yield surfaces. In addition to that, after the calculation phase in which in situ field stresses are determined, both soil strain history and soil displacements have been reinitialized. Numerical results and field measurements of the extensometer Nr. 3 are shown in Figure 17 (other simulations results can be found in Castellón 4 ). As it can be seen, the simulations conducted with the EPHYSS or HS-SS models are not able to reproduce the upper part of the displacement profiles measured with the extensometers, in which a remarkable reduction in vertical soil displacements is appreciated. This behavior can be explained by the progressive drying of the soil near the water table as the excavation progresses, which generates gradual increases in suction in the unsaturated zone of the soil and, therefore, a shrinkage of it. 6 According to the simulations performed with the Barcelona Basic Model (BBM), 63 the value of these displacements caused by shrinkage can be of the order of 5 mm, 6 thus approaching the observations. In any case, the use of the EPHYSS model allows to improve the approximation to the measures taken with the extensometers with respect to the values calculated with the HS-SS model. In Figure 17 the results of the simulations with the EPHYSS model approximate quite well the measurements of the extensometer Nr. 3 in the different calculation phases, especially in the transversal cross-section of the drained analysis. Furthermore, simulations with the HS-SS model approximate these measures better or worse depending on whether, in a certain calculation phase, the results of the transversal or It is important to point out the effect of the inconsistencies generated by the consolidation phases in the undrainedconsolidated analysis with HS-SS model that can be observed in Figure 17. They are due to the reinitialization of the components of the HS-SS model history tensor ( ), giving place to markedly lower displacements than those obtained in the corresponding drained phases. Equivalent reinitialization of the components of EPHYSS model history tensor ( ) happens, however, this problem is solved in EPHYSS by defining new state variables ( MEM , MEM and MEM ). The effect of such inconsistencies in the HS-SS model is especially evident in the results of phases B and C. On the other hand, note that part of the difference between the displacement profile of the drained and undrained-consolidated analysis corresponding to phase D, both in the HS-SS model and in the EPHYSS model, is because no consolidation has been considered in such phase after the application of the corresponding undrained loading.
Inconsistencies in the simulations with the HS-SS model can be easily detected in the graphs ap ∕ ,ur corresponding to the consolidation subphase of the undrained-consolidated analysis, when compared with the same graphs corresponding to the respective undrained loading subphases of the undrained-consolidated analysis, which are shown in Figure 18 for the longitudinal cross-section. It can be clearly seen how both the EPHYSS and the HS-SS models have numerical reversals during the consolidation subphase that give place to a soil stiffening. Nevertheless, unlike the HS-SS model ( Figure 18B), in which these reinitializations influence the subsequent soil history, in the EPHYSS model ( Figure 18A) this effect is corrected and does not accumulate, having no influence in following phases.
In Figures 19 and 20 it can be seen, respectively, how the maximum value of the excess water pressure generated during the undrained loading subphases of the undrained-consolidated analysis and the distributions of such excesses are very Finally, it must be pointed out that the calculation time when solving boundary value problems with EPHYSS, using conventional processor, falls within the reasonable values for a commercial geotechnical software.

CONCLUSIONS
The behavior of the soil in the range of small strains should be always considered in the analysis of geotechnical problems when sensitive constructions are affected, which is very common in urban environments. The EPHYSS model described in this paper, which is composed by the HQH model and the Cap-Cone HS MOD model, is capable to reproduce the quasistatic behavior of the soil in Zones I, II, III and IV of Jardine. 1 Zones I and II, where soil behavior is nonlinear reversible, hysteretic and dependent on the recent history, are described by the HQH model, which also considers the strain-induced anisotropy, while Zones III and IV of Jardine are described by the HS MOD model. EPHYSS model inherits some important limitations from the HS MOD model that makes it unsuitable in numerical analyses of geotechnical works in which intermediate or large strains play a significant role. However, in many geotechnical works in urban areas, under safe conditions and far from failure, smalls strains dominate and therefore EPHYSS model is appropriate. The reversible part of the EPHYSS model considers two strain domains. The elastic bulk modulus is common in both domains. A degradation law of the shear modulus is considered in domain 1 until it reaches a minimum value in domain 2. Poisson's ratio is variable in both domains, and always greater than ′ min . Furthermore, the EPHYSS model considers 10 state variables: 8 of these variables correspond to the reversible part of it (HQH model) and 2 to its plastic part (HS MOD model). These HQH model state variables define different short and long-term memory levels that provide the EPHYSS model with robustness for the reproduction of soil hysteretic behavior in quasi-static problems, using in some of them a layer structure similar to that proposed by Hueckel and Nova. 33  Different oedometric, triaxial and biaxial tests compiled from various thesis and papers have been simulated to carry out a partial verification of EPHYSS model, as well as a validation and a comparative analysis with the HS-SS model. Both models provide a very good approximation to the experimental data. When no strain rotations take place or when reversals are total ( = 180 • ), EPHYSS model provides similar results to those of the HS-SS model, except for some differences related to stiffness stress dependency, stiffness moduli dependencies, dilatancy formulation and the consideration of strain induced anisotropy. However, there are some significant differences between both models when the reversals are partial ( < 180 • ). From the numerical biaxial tests in which the results of the simulations with the EPHYSS model are compared with those of the HS-SS model, the hypoplastic model with intergranular strain of Niemunis and Herle 22 and the multilaminated model of Schädlich and Schweiger, 29 it is concluded that the HS-SS model is not capable of reproducing stiffness recovery when the cumulated strain value prior to the rotation is high. The latter does not occur in any of the other three advanced models, where the soil elastic recovery is less dependent on the value of the previous accumulated strain. This represents an added advantage of the EPHYSS model, as it is a model easy to use in the professional practice.
Finally, it is concluded from simulations of numerical oedometric and triaxial tests in which consolidation phases, nil phases or phases with small unloading/reloading have been introduced, that the EPHYSS model is capable to correct the effect of the inconsistencies detected in the HS-SS model, [4][5][6][7] which are cumulative and can have very important effects on the results of boundary value problems. The effect of the correction of HS-SS inconsistences is also shown in a large urban excavation that has been studied, which corresponds to the works of the future intermodal station of La Sagrera, in Barcelona. The soil parameters have been obtained from the tests carried out during the project phase of the station, from empirical correlations and from 20 resonant column tests conducted in cylindrical samples obtained from highquality block samples. Four extensometers were installed, reaching 60 m depth below the excavation bottom to measure soil uplift. Numerical simulations of the excavation have been conducted in four transversal cross-sections that contain, each of them one extensometer, and in a longitudinal cross-section that contains all of them. The EPHYSS and HS-SS models have been used and the results obtained with them have been compared. In the study conducted, two types of analysis were performed, a drained analysis and an undrained-consolidated analysis. From the results of the simulations with the EPHYSS and HS-SS models it is concluded that none of them can reproduce the upper part of the displacement profiles measured with the extensometers, in which an evident reduction in the vertical displacements is observed. This behavior can be explained by the progressive drying of the soil near the water table as the excavation progresses, which generates gradual increments in suction in the unsaturated area of the soil and, therefore, a shrinkage of it. 6 In any case, it is concluded that simulations with the EPHYSS model significantly improve the approximation to the measurements of the extensometers in the different calculation phases if compared with the results of the simulations obtained with the HS-SS model, both in the transversal and in the longitudinal cross-sections. In the results of the simulations with the HS-SS model, the effect of the inconsistencies generated by the consolidation subphases in the undrained-consolidated analysis can be clearly seen, which leads to displacements markedly lower than those obtained in the corresponding drained phases. The EPHYSS model is able to correct the effect of these inconsistencies thanks to the model state variables. Finally, it can be seen that both, the distributions of excess water pressure and the maximum value of such excess, are very similar in both models.
It is concluded that, in general, the EPHYSS model significantly improves the approximation to the experimental measures with respect to the HS-SS model, especially in those boundary value problems that present partial reversals in the deviatoric strains or high cumulated strain values prior to reversals; resolves the inconsistencies of the latter with a reasonable computational cost; and requires simple parameters, most of them common to those of the HS-SS model. All this makes EPHYSS a model that can be used for analysis and design in geotechnical professional practice.

LIST OF SYMBOLS AND ABBREVIATIONS
Parameter of the shear modulus expression according to the model of Dos Santos and Correia 18 whose value is = 0.385 and which is used in the SSOM, HS-S, HS-SS, HQH and EPHYSS models. ′ , ∕ ′ Internal parameter in the HS-S, HS-SS and EPHYSS models related with the value of ′ that represent the ratio between the elastic bulk modulus and the secant bulk modulus for the primary isotropic compression. 0 Coefficient of lateral earth stress for a normally consolidated stress state ( 0 = 1 − sin( ′ ) by default). ( ′ , e) Fourth-order linear tensor in a hypoplastic model that introduces void ratio as a state variable. 22 Coefficient that controls the dependence level of ′ 50 , ′ , ′ , and ′ with the stress in the HS-S, HS-SS and EPHYSS models. Soil parameter that controls stiffness value before a rotation of the strain path of 180 • in the model of Niemunis and Herle. 22 Soil parameter that controls stiffness value before a rotation of the strain path of 90 • in the model of Niemunis and Herle. 22 1 Soil parameter that controls the dependence of ′ on − ′ in the HQH and EPHYSS models. 2 Soil parameter that controls the dependence of and , on − ′ in the HQH and EPHYSS models. Soil parameter in the Li and Dafalias dilatancy formulation. 36 Critical stress ratio. Internal parameter that determines the maximum number of reversal points in which is possible to memorize ‖ ‖, ‖ ‖ and ‖ , ‖ in the HQH and EPHYSS models. Reference mean stress in the expression of ′ in the HQH and EPHYSS models.
Deviatoric soil strength given by the Mohr-Coulomb criterion. Parameter with the maximum value of ‖ ‖ in the model of Niemunis and Herle. 22 Soil parameter that represents the ratio ∕ ( = 0.9 by default). Active reversal points in the HQH and EPHYSS models.
Deviatoric strain tensor in the last reversal point that conforms the endpoint of the active strain cycle in the HQH and EPHYSS models. Soil saturation degree. Water cavitation stress. admitted tolerance. admitted tolerance. Imposed displacement in the global x-direction. Imposed displacement in the global y-direction. * Water content. Numerical parameter in the HQH and EPHYSS models that controls the speed with whicḣevolves. = arcos(̂∶ˆ̇). * = arcos(̂∶ˆ̇) * is the rotation angle in the recent total deviatoric strain path from which reversals appear in the HQH and EPHYSS models. Rotation angle of the incremental total strain.

A C K N O W L E D G E M E N T S
The authors acknowledge the partners involved in the construction of La Sagrera station, in particular ADIF (Administration), INECO (Construction Management), Audingintraesa-Ayesa-Cicsa (Technical Assistance), Dragados-Acciona-Comsa-Acsa (Construction companies) and Geocisa (Monitoring company) for their support and for providing field measurements and field and laboratory data required in the analyses. The authors also acknowledge the support from Christian Hoffmann in charge of the Resonant Column Tests. Finally, the collaboration from Lilia Maria Grover Pimienta in the final edition the manuscript is gratefully acknowledged.

D ATA AVA I L A B I L I T Y S TAT E M E N T Data available on request from the authors.
O R C I D Javier Castellón https://orcid.org/0000-0002-7022-9316 Alberto Ledesma https://orcid.org/0000-0003-3321-3849