Refinement of the Hardening Soil model within the small strain range

The popularity of the elasto-plastic Hardening Soil (HS) model is based on simple parameter identification from standard testing and empirical formulas. The HS model is implemented in many commercial FE codes designed to analyse geotechnical problems. In its basic version, the stress–strain behaviour within the elastic range is subject to the hypoelastic power law, which assures the barotropy of the elastic stiffness. However, a proper modelling within the small strain range, i.e. strain-induced stiffness degradation and correct reproduction of the hysteretic behaviour, was one of the most important drawbacks in the HS formulation. The first small strain stiffness extension to the HS model was proposed by Benz (Small strain stiffness of soils and its numerical consequences, 2007), and the new model was called Hardening Soil Small (HSS). Despite the simple isotropic formulation, its applicability was proved in various numerical simulations in geotechnics. However, the HSS formulation exhibits a serious fault known in the literature as overshooting, i.e. uncontrolled reset of the loading memory after tiny unloading–reloading cycles. The authors' main aim was to retain the set of material parameters for the HSS formulation and to propose a new small strain extension to the HS model without overshooting. The new proposal is based on the Brick model which represents the concept of nested yield surfaces in strain space. The implementation aspects of the new HS-Brick model are described, and its performance is presented in some element tests and selected boundary value problems by comparisons with the HSS formulation.

Axial and radial strain components in axisymmetric condition, respectively r; r ij Effective Cauchy stress tensor, compression positive r 1 ; r 2 ; r 3 Major, intermediate and minor principal stress components, respectively r a ; r r Axial and radial stress components in axisymmetric condition, respectively /; / cs Effective and critical state friction angle, respectively w Dilatancy angle m Poisson's ratio x G Shear stiffness modulus reduction factor in the HS-Brick model Dx b G Proportion of x G assigned to b-th brick,

Introduction
The basic version of the elasto-plastic Hardening Soil (HS) model was developed by Schanz et al. [36] and implemented into the Plaxis FE code [7]. The HS model has become a standard in geotechnical engineering computations. It is implemented in many commercial FE codes due to clear parameter identification based on the standard laboratory and field tests or empirical formulas. The model assembles well-established concepts in soil mechanics concerning the barotropy, dilatancy, shear strength and the nonlinear pre-failure stiffness. The behaviour of granular, fine grained and even soft organic soils is possible to be simulated with the HS model. The HS formulation has also some limits and simplifications which should be taken into account in the numerical analyses of geotechnical problems. The main framework of the HS model may be explained by drawing the contours of yield surfaces on the plane of Roscoe's stress invariants p-q. It is shown in Fig. 1. The Mohr-Coulomb shear strength criterion limits shear stress levels, whereas the volumetric and pre-failure deviatoric hardening is controlled by two yield surfaces. The HS model is entirely isotropic in both the elastic and elasto-plastic ranges. Hence, no inherent and stress induced anisotropy is available contrary to the experimental evidence especially in fine grained deposits, e.g. [26,27]. Non-viscous formulation excludes also simulation of the time effects in soils like creep or relaxation. The stiffness in the elastic range is represented by the hypoelastic power law with the stress-dependent unloading-reloading Young's modulus E ur and constant Poisson's ratio m. This limits the application of the HS model in the dynamic analyses [30,45]. However, one of the most important drawbacks in the initial HS model formulation was the lack of proper modelling of the small strain behaviour, i.e. strain-induced degradation of the high initial stiffness during monotonic loading, regain of the high stiffness after sharp loading reversals and correct reproduction of the hysteretic behaviour. The practical importance of the small strain stiffness incorporation in the constitutive models of soils was widely approved in literature [1,3,8,10]. Hence, it was a necessity to refine the HS model to comply with this evidence. The first small strain stiffness extension to the HS formulation was proposed Benz and co-workers [5,6], and the new model is called Hardening Soil Small (HSS). The new formulation introduces a special algorithm of controlling the current hypoelastic stiffness within the elastic range of the basic HS model. Additional parameters can be obtained from extended standard tests or empirical formulas. Despite the simple isotropic formulation, its applicability was proved in various numerical simulations in geotechnics, e.g. [4,5,19,20,22,33,41]. The extended HSS model is also implemented in commercial geotechnical FE codes available for researchers and practitioners. However, the HSS formulation exhibits serious fault known in the literature as overshooting, i.e. uncontrolled reset of the loading memory and regain of high initial stiffness after tiny unloading-reloading cycles. It is illustrated schematically in Fig. 2. The overshooting problem in the HSS model is critically reviewed in details in [31].
The primary aim of this paper is to propose a new small strain extension without the overshooting and to retain the main modelling assumptions of the HSS formulation together with its set of material parameters. The new proposal is based on the modified version of the Simpson's BRICK model [38] which practically represents the concept of nested yield surfaces in strain space [29]. We describe the implementation aspects of the new HS-Brick model and present its performance by comparisons with the HSS formulation in some element tests and selected example boundary value problems.

Summary of the basic HS model
The basic HS model belongs to the class of multi-surface elasto-plastic models. The deviatoric and the volumetric cap yield surfaces, with the corresponding plastic flow rules, and hardening laws, are the two major plastic mechanisms introduced to represent the nonlinear soil behaviour. The ultimate limit states are controlled by the Mohr-Coulomb and Rankine strength criteria. The main model equations and theoretical assumptions are systematized in [36]. Since the HS model is implemented in several FE codes, we would like to emphasize only the main characteristic aspects of its formulation in the currently developed version of ZSoil FE code [44]. This version of the basic HS model is refined to account for the small strain stiffness and used in the calculations presented in the article.

Deviatoric plastic mechanism
The yield surface controlling the deviatoric hardening mechanism is expressed as follows [36]: where the asymptotic and ultimate deviatoric stresses q a , q f are defined as: and R f \1:0 is the failure ratio. We also assume that the minimum initial value of the hardening parameter c p 0 ! 10 À6 and r 1 ! r 2 ! r 3 . The plastic flow rule is derived from the following plastic potential function: The mobilized dilatancy angle w m is defined after Rowe [35] as: where mobilized and critical state friction angles / m ; / cs are calculated from: It is worth noting that D ¼ 1 in the reference Rowe's formula [35]. In general, when sin / m \ sin / cs , the contractancy cut-off condition is assumed and sin w m value is scaled by the parameter D 1:0. In the basic HS model D ¼ 0, whereas in the extended version with small strain stiffness refinement D % 0:25. The hardening law in the deviatoric mechanism is given by: where dk 1 is the plastic multiplier corresponding to the shear plastic mechanism.

Volumetric plastic mechanism
The smooth cap yield surface is described by the following equation: In the above expression rðhÞ is defined after van Ekelen's [13] formula: where ð11Þ The power exponent n ¼ À0:229 and the parameter a 0:7925. A non-associated plastic flow rule is derived from the following potential function: The above form is crucial for the efficient implementation scheme in which the stress return algorithm is executed in the principal stress space. It is worth noting that the plastic potential is obtained from the yield function f 2 (Eq. 9) by fixing rðhÞ ¼ 1. One can set the rðhÞ to any fixed value in the range k. . .1 (Eq. 13). The volumetric hardening law is expressed in the following form: where the volumetric plastic strain increment de p v2 is generated entirely by the cap mechanism. The internal model parameters M and H are derived taking into account the assumed K NC 0 value and the reference tangent oedometric modulus E ref oed defined at a given vertical stress assuming the normal consolidation.

Stiffness barotropy
In the basic HS model [36], the unloading-reloading E ur and secant E 50 stiffness moduli are stress dependent by the following power laws and the corresponding empirical barotropy function f r : However, the above form exhibits certain drawbacks, e.g. when deep excavations or highly overconsolidated soil layers are analysed. In such conditions, horizontal stresses are usually larger than the vertical ones. Therefore, we prefer to use the following barotropy function based on the mean effective stress: In both cases, a certain cut-off condition has to be applied for the minimum value of stress invariants r 3 or p to avoid too low stiffness moduli.
3 Small strain stiffness extensions to the HS model

HSS formulation
The first refinement of the basic HS model to account for the realistic small strain stiffness of soils was proposed under the name Hardening Soil Small (HSS) [5,6] and implemented into the Plaxis FE code [7]. The main idea was to substitute hypoelastic kernel of the HS model (constant Poisson's ratio m ur and stress dependent E ur ) by a model allowing to reproduce the small strain stiffness properties [3,10,15,21,42]. These are high initial elastic stiffness of soils represented by the small strain shear modulus G 0 , strain degradation of the high stiffness in monotonic loading and changes of the stiffness adequate to the current loading direction, e.g. regain of the high stiffness after sharp loading reversals. The stiffness strain degradation is often described by the so-called S-shaped curve relating the current shear modulus to the selected strain measure. Among many proposals from the literature, the modified Hardin-Drnevich relationship [16] is applied in the HSS model where the current secant shear modulus G s is related to the shear strain invariant c: The new important small strain threshold parameter c 0:7 is the shear strain needed to degrade G 0 down to 0:7G 0 and a is another parameter influencing the S-shaped curve. The fixed value a ¼ 0:385 is proposed. Decay of both the tangent and secant shear stiffness resulting from Eq. 21 is limited by a strain threshold beyond which the target stiffness achieves its low cut-off value G ur . As such, shear stiffness tends towards a fixed positive constant value. More complicated modelling aspect is connected with controlling changes of the current stiffness due to loading reversals or changes of the loading direction. To this end, the strain history second-order tensor H is proposed. The evolution of H is related to the increments of the deviatoric strain De. The original strain history tracing algorithm responsible for the evolution of H is formulated in [5] together with the implementation code. On the basis of strain history evolution and current deviatoric strain increment, the shear strain measure c Hist is defined and substituted into Eq. 21 instead of c. In the monotonic loading, c Hist is equivalent to c. However, the history tracing algorithm used in the HSS model was recently reanalysed and critically reviewed in [31]. It is shown that the reset of the material memory resulting in regaining of the high initial stiffness may be induced by tiny unloading-reloading loops, and hence, the reproduced stress-strain behaviour will depart from the expected. This modelling fault is often denoted in the literature as overshooting. We refer the interested reader to the article [31] for comprehensive analysis of the overshooting problem in the HSS model.
In the rest of the following section, we describe some important issues of the HSS model implementation used in the presented calculations. Similarly to the reference approach [5,6], we also trace current stiffness ratio parameter G m being the minimum of the fractionG s =G ur over whole loading history (notion of the incremental secant stiffness modulusG s is explained in Box 3). This state parameter is needed to detect transition from the unloadingreloading branch to the virgin loading one. Moreover, it is used in the modified hardening laws of the deviatoric and volumetric mechanisms, via h i function, as follows: Depending on the current loading status, which can be the unloading-reloading or virgin loading, an incremental secant shear modulusG s is derived from Eq. 21. This modulus is used later on to compute the elastic trial state but also the plastic corrector. The algorithmic treatment of stiffness barotropy is discussed in Sect. 4.1. It should be emphasized that these two aspects are not clearly explained in [5,6].

HS-Brick formulation
The main aim of the reported work is to develop a new small strain formulation without the overshooting which may substitute the procedure based on the strain history tensor H in the HSS model. However, we want to keep the set of HSS material constants and the basic assumptions of this model. These are hypoelastic stress dependency of the actual tangent shear modulus G t based on the power law (e.g. Eq. 20) and parallel incorporation of the strain degradation of reference small strain shear modulus G ref 0 according to Eq. 21. We decided to employ the concept of nested circular yield surfaces by Mróz [29]. Since it is convenient to keep the control of small strain stiffness within the strain space, we have followed the original version of nested yield surfaces in strain space proposed by Simpson [38,39] in his BRICK model. This model is described by the analogy of a man pulling the finite number of bricks on strings. Initially strings of different lengths are slack. When man moves monotonically in a given direction strings become taut one by one pulling the following bricks. The man's movement represents strain, and the strings lengths are radii of circular yield surfaces in the strain space. Every time the next brick starts to be pulled by man, the reference stiffness is degraded in the stepwise fashion. Hence, the S-shaped curve is not reproduced as continuous. When the loading direction in strain space changes loosing the strings, the high initial stiffness is regained. Different versions of the BRICK model are reported in the literature giving satisfactory results in various FE computations, e.g. [9,11,14,23,25,40,43]. The original BRICK model [38] has been developed in plane strain for which volumetric strain and shear strain are appropriate axes-the sides of the room where the described pulling process takes place. In the 3D extension of the original BRICK model [25], these axes were replaced with volumetric strain axis and five shear strain axes. Additionally, the pulling process of a given brick is enhanced by the so-called plastic strain reduction which provides volumetric strain hardening by repositioning the bricks [14]. In our version of BRICK model, we place the pulling process within the general six-dimensional strain space. However, the relative strain distances between the man and the b-th brick are measured by the shear strain invariant c b of the strain distance. In such a way, the pulling process is projected on the deviatoric plane in strain space. This is shown for various loading situations in Fig. 3. The obtained strain degradation accords with that assumed in the HSS model. The only difference is the fact that continuous S-shaped curve from the HSS model is substituted by the stepwise model presented in Fig. 4.
Summing up, within the elastic range of HS model the current reference tangent shear modulus G ref t changes between the values of G ref ur and G ref 0 . This is controlled by the BRICK type strain history update procedure in the stepwise way. Additionally, barotropy of small strain stiffness is taken into account in parallel and final actual tangent shear modulus G t is both strain and stress dependent. The proposed new formulation of the HS model accounting for the small strain stiffness is named the HS-Brick model. The bricks strain history update procedure is given in Box.1. In the presented expressions, the string length of bth brick is denoted by s b . Any quantity with the right lower subscript n refers to the last configuration of equilibrium while the one with the index n þ 1 to the new configuration which is currently sought. The additional right superscript i refers given quantity to the current iteration in the nonlinear iterative procedure while i À 1 to the previous one.
Although it is formally designed as a one step procedure (the full strain increment De ij;nþ1 is used) a sub-incrementation scheme, discussed later on, can easily be adopted. In the latter case, the De ij;nþ1 has to be treated as a fraction of the full strain increment, while bricks strain history e b ij;n and strain e ij;n correspond to the end state of the previous sub-step. The resulting e ij;nþ1 , e b ij;nþ1 correspond to the end state of the current sub-step. This procedure enables to compute the current reference incremental tangent shear modulus G ref t . This modulus degrades with the number of bricks that are currently pulled. It is worth to note that the minimum reference tangent shear modulus is delimited by the material parameter G ref ur .

General implementation scheme
Numerical efficiency of the stress-strain integration scheme is the key when dealing with practical initial boundary value problems. Barotropy and its algorithmic treatment as well as the specific form of the yield surface controlling the deviatoric mechanism yield major difficulties when implementing the new extended version of the HS model. Therefore, we discuss all these aspects in the following subsections.

Algorithmic treatment of barotropy
The general implementation scheme, as in the case of most elastic-plastic constitutive models, consists of three major steps, i.e. the trial stress evaluation r tr ij , checking its plastic admissibility (f k ðr tr ij ; c p n ; p c;n Þ 0 for k ¼ 1. . .N) and then carrying out the stress return procedure, if at least one of N plastic yield surfaces becomes active. The major problem in the HS model is related to the treatment of stiffness stress dependency called here as barotropy.
In several publications, researchers try to derive the algorithmic tangent stiffness operator which should preserve quadratic convergence rate of the iterative scheme used to solve a nonlinear system of equilibrium equations. In the considered model formulation, deriving such an operator is not really possible. It is mainly due to non-smooth constitutive functions appearing in the cut-off condition in barotropy functions but also in the brick type approximation applied to the stiffness degradation function. As it will be pointed out at the end of this section, the stress return algorithm can be carried out in the principal stress space by freezing elastic stiffness moduli. This seems to be optimal from the computational efficiency point of view.
In the HS model, no matter if it is extended to the small strain regime or not, barotropy is present at the stage of trial stress evaluation, in the analytical expression for the shear plastic yield surface and in the plastic corrector procedure. In order to preserve all benefits of implicit integration schemes, barotropy cannot be treated in the same manner in all the aforementioned three stages. Let us define the two different reference stress states influencing barotropy: Using the stress state form (i À 1) iteration, we treat barotropy in a semi-implicit manner. Assuming h ¼ 1 2 the current elastic stiffness operator D e ijkl ðr ij;nþh Þ represents the incremental secant stiffness. When strain increments are very small it reduces to the tangent elastic operator. The D e ijkl operator is calculated using current incremental Young's modulus E, i.e. stress dependent E ur in the basic HS model or the actual E s when the small strain extension is activated.
However, the reference stress state influencing barotropy in the algorithmic stiffness treatment cannot be used in the analytical condition defining the deviatoric hardening yield surface in Eq. 1. The stress-dependent moduli E ur and E 50 should be calculated usingr ij;nþ1 . To explain this, let us assume that we start a new time step n þ 1 with a very small strain increment De ij . This means that r ðiÞ ij;nþ1 % r ij;n andr ij;nþh % r ij;n . If we use the barotropic reference stress to verify the deviatoric yield condition then the current stress state, previously satisfying the condition f 1 ðr ij;n ; c n Þ % 0 now can easily violate it and from f 1 ðr ij;n ; c n Þ [ 0 some additional plastic straining will be induced. Moreover, in the strain history update procedure of the HSS model, the strain history tensor H can be reset. The small strain formulation in the HS-Brick model is insensitive to these faults.
The above algorithmic treatment of barotropy enables application of the classical operator split method and to conduct the stress return procedure in the principal stress space. Furthermore, due to the specific form of the assumed plastic potentials the multi-surface stress return problem may be solved in an analytical manner in the case of active deviatoric plastic mechanism and optionally together with Rankine one.

Elastic stress predictor in the HS-Brick model
As the S-shape stiffness degradation curve in the HS-Brick model is discretized with piece-wise constant segments a certain sort of strain increment subincrementing is required to compute the elastic stress predictor. The proposed subincrementation scheme is summarized in details in Box.2

Multi-surface stress return in the principal stress space
In the case of the HS model, due to the specific form of assumed plastic potentials, it is possible to perform stress return in the principal stress space, using the approach proposed in [24]. Additionally, the stress return can be solved analytically in all cases when the final stress state remains on the deviatoric yield surface or in the point of its intersection with the Rankine surface. In all these cases, if we keep constant sin w m value, the stress return procedure reduces to a scalar quadratic equation for unknown plastic multiplier Dk 1 . This is crucial as the deviatoric yield surface is undefined for q ! q a . If the final stress state violates the cap yield condition, an iterative procedure is executed in the principal stress space with the following vector v of independent variables: while the vector r of update equations residuals is defined as follows: The JACT is a current set of N JACT active plastic mechanisms. In order to handle intersection of deviatoric yield surface with Rankine one, the additional two constraints (surfaces f 3 and f 4 are reserved for the Mohr-Coulomb and Rankine criteria) are added [24]: The nonlinear set of update equations can be solved using Newton-Raphson iterative scheme supplemented by a linesearching procedure. The initial values of independent variables are derived from the analytical solution where cap surface is neglected. Mapping from the reduced principal stress space to the full one is carried out based on the eigenvectors of the trial stress tensor.

Element tests
In the following section, some basic element tests (stress and strain fields are homogeneous) in the axisymmetric stress conditions are performed to examine the proposed HS-Brick model in the cases which display overshooting problems in the HSS formulation. Generally, the overshooting occurs independently of the material parameter values. Therefore, we adopt the common set of material constants of glacial till used in [5] for all element tests presented in the article: The FE model for axisymmetric element tests of dimensions 1Â1m is meshed with 1 quadrilateral enhanced assumed strain element (Q4 EAS elements [12,37] are used in all presented calculations). The vertical displacements at the bottom are fixed and the horizontal fixities at the symmetry line are applied. The uniform loading is controlled separately on the top and outer side of the sample-axial r a and radial r r stress components, respectively. The soil is weightless. The initial preconsolidation pressure is set to p c ¼ 200 kPa, and initial position of the yield surfaces corresponds to the K NC 0 value. The initial strain is set to zero, i.e. the history tensor H is zero in the HSS model, whereas in the HS-Brick model all bricks are set together with the initial strain to zero. In the first phase of all presented axisymmetric element tests, the sample is consolidated isotropically to p 0 ¼ 100 kPa and as a result OCR ¼ p c =p 0 ¼ 2:0. As during the isotropic compression of isotropic material no shear strains are induced, there is no difference in the mechanical response of both compared models. Therefore, we present only the results of numerical calculations after the isotropic consolidation phase with resetting the strains to zero on the graphs. In all element tests, we use time-dependent load drivers. Every single loading phase is performed in one day periods, and equal time increments Dt ¼ 0:01 day are used in calculations.

Triaxial drained compression with small unloading-reloading loops
In real soil behaviour, when a monotonic loading is interrupted by a small single unloading-reloading loop the continuation of the loading leads to sweeping the effect of short stiffness increase within the loop. The subsequent response is practically the same like in the case without interruptions. In routine numerical analyses of geotechnical problems, the small loading reversals may be never intended to include as a separate calculation phase; however, the tiny unloading-reloading loops may occur in some domains of the model without control and it is extremely important that such occurrences will not influence the results of modelling the soil structure interaction.
In the current example, we show the influence of five small unloading-reloading axial stress loops (Dr a ¼ AE5 kPa) performed during drained triaxial compression test within the pre-failure range. The compressions with small reversals are compared with analogous purely monotonic tests for the HSS and HS-Brick models. The details and results of these tests are presented in Fig. 5. The first small unloading-reloading loop is performed after reaching r a ¼ 120 kPa. Then, the loops are executed after each of four axial stress increases of 30 kPa, i.e. when r a ¼ 150; 180; 210; 240 kPa. The compression is continued up to the shear failure determined by Mohr-Coulomb criterion which occurs at r a ¼ 297 kPa. Responses due to monotonic loading illustrated by axial strain-stress curves in Fig.5 are equal from both models. In the case of compression interrupted by loading reversals, the overshooting causes uncontrolled increases of stiffness in the HSS model and the response diverges from analogous monotonic compression. It is highly evident in the small strain range magnified by logarithmic strain scale before the stress path reaches the yield surfaces of the basic HS model. In the HS-Brick model, the small unloading-reloading loops do not influence the general response as expected.

Triaxial hysteretic behaviour
In constitutive modelling of soils, the purely hysteretic behaviour, i.e. without cumulative effects, is modelled with different nonlinear elastic models. The most popular formulations are based on the classical 1D models [16,34] which follow the Masing rules [28]. However, generalization to 3D modelling is not straightforward and requires some special techniques in order to properly simulate hysteretic behaviour. It is very important to appropriately define the loading reversals as well as the actual stiffness during the reloading and overloading. Firstly, the reloading curves on the strain-stress graphs need to pass exactly through the old loading reversals. Secondly, smaller cycles need to be swept out from the material memory due to overloading. The purely hysteretic behaviour in general 3D cases is shown to be correctly reproduced in the paraelastic models [18,32] or by using the multi-surface plasticity models with kinematic hardening [29]. The latter method is actually applied in the HS-Brick model.
We perform two element tests to validate the hysteretic behaviour reproduced by the HSS and HS-Brick models. In the first element test, a single hysteresis loop is interrupted by small unloading-reloading cycle. The results are shown in Fig. 6. To obtain the reference single symmetric hysteresis, the sample is loaded up to r a ¼ 180 kPa, unloaded down to r a ¼ 50 kPa, reloaded again up to r a ¼ 180 kPa and overloaded up to r a ¼ 185 kPa. In both models, the Fig. 5 Drained triaxial compression simulation results with HSS and HS-Brick models. In the HSS model, small axial stress unloading-reloading (ur) loops Dr a ¼ AE5 kPa cause overshooting within the elastic range. In the proposed HS-Brick model, the small interruptions of the monotonic loading do not influence the general course of the stress-strain response as expected. In the graphs at the bottom, the logarithmic scale of strain is used to magnify the changes within the small strain range where overshooting is clearly evident in the HSS model response is equal and the obtained loop is closed. To examine the overshooting problem, the single hysteresis loop is interrupted by small axial stress unloadingreloading loop Dr a ¼ AE10 kPa at r a ¼ 130 kPa as shown in Fig. 6. The interruption causes overshooting in the HSS model, and the main loop is not closed deviating from the expected paraelastic behaviour. The response obtained with the HS-Brick model is correct, and the main loop is closed regardless the small loading reversal during reloading. In the second element test, more complicated loading programme is applied in the so-called shakedown test. The obtained results are presented in Fig. 7. Similar to the previous test, regarding hysteretic behaviour, the sample is initially loaded up to r a ¼ 180 kPa (point 1) and unloaded down to r a ¼ 50 kPa (point 2). Then, axial cyclic unloading-reloading is applied with subsequent changes of the stress amplitude after each loading reversal by 10 kPa, i.e. r a ¼ 170; 60; 160; 70; 150; 80; 140; 90; 130; 100; 120; 110 kPa marked in Fig. 7 by points 3 to 14, respectively. After unloading to point 14, the reloading back to the initial reversal at r a ¼ 180 kPa and overloading up to r a ¼ 185 kPa is simulated. In the desired hysteretic behaviour, the last reloading from point 14 should pass through old reversals (subsequently points 13,11,9,7,5,3,1) which is fulfilled only in the case of the new HS-Brick model. In the HSS model, the main hysteresis loop is not closed and point 15 is reached due to overshooting.

Simple geotechnical boundary value problems
It is important to present the overshooting problem and to validate the proposed HS-Brick model not only at the level of representative elementary volume but also in the numerical simulations of some real geotechnical problems, which are analysed in the routine design using FE codes.
We have prepared and analysed three different simple boundary value problems concerning the loading of a shallow foundation, execution of supported excavation and dynamic analysis of the shear layer problem during earthquake. In the shallow foundation and excavation problems, the parameter set of glacial till from Eq. 44 is used. Additionally, the dry and effective weights of the soil are c ¼ 19:0; c 0 ¼ 9:0 KN/m 3 , respectively. The overconsolidation is prescribed by the pre-overburden pressure POP ¼ 200 kPa and the initial position of yield surfaces corresponds to the K NC 0 value. Initial value of K 0 ¼ 0:9 is constant within the glacial till deposit according to the applied overconsolidation. In both problems, we include the shallow dry layer of granular fill simulated with the Mohr-Coulomb model The purpose of introduction of this layer is to prevent activation of the cut-off procedure for the minimum value of stress invariants influencing current stiffness moduli like in Eq. 20.  Fig. 8, i.e. monotonic loading, monotonic loading interrupted by one or four small unloading-reloading loops. The aim is to inspect differences in soil ground response in a similar manner like in the element test presented in Sect. 5.1. The results are shown in Fig. 9 by comparing the foundation load-settlement curves. The effect of overshooting in the calculations with the HSS model is much higher than that observed in the analogous element test. It is due to the fact that in practical geotechnical problems concerning the prefailure range (serviceability limit states), the small strain stiffness is of the highest importance. The settlements obtained with the HSS model for the monotonic loading interrupted by small reversals are significantly smaller when compared with the monotonic loading, i.e. one reversal results in ca. 40%, whereas four reversals result in 70% settlement reduction. The results obtained with the HS-Brick model display no influence of the monotonic loading interruptions. The response of both models due to monotonic loading is practically the same as expected. Axial stress amplitude is subsequently decreased after each loading reversal by 10 kPa. In seventh loading cycle, the sample is reloaded to the axial stress of reversal 1. In the desired hysteretic behaviour, the reloading passes through old reversals, i.e. after reversal 14 it should pass old reversals 13, 11, 9, 7, 5, 3 reaching finally reversal 1. This is achieved in the HS-Brick model only. In the HSS model, the point 15 is reached instead. In the bottom graphs, the reloading path is magnified The small unloading-reloading loops applied in the current examples may seem artificial in the real calculations; however, it allows us to present the influence of the overshooting clearly. In practice, the situation may be worse-the reset of the material memory may occur in any staged construction FE geotechnical modelling without notice, e.g. during replacement of materials, partial excavations, dewatering, installation of anchors, etc.

Unloading case-excavation
An interesting effect of the overshooting may be observed in the numerical simulations of excavation with the simultaneous consolidation. This is available in the fully coupled deformation and flow analysis. The FE model of the exemplary excavation problem is shown in Fig. 10. The excavation of 4.0 m depth is supported by 8.0 m long sheetpile wall (beam elements with interfaces, EI ¼ 461250 kNm 2 , EA ¼ 6:15 Â 10 7 kN, m ¼ 0:15). During the initial Dt ¼ 1:5 days, the sheetpile wall elements with interfaces are activated in the model. After the first phase of excavating the fill material layer, the struts are introduced by fixing the horizontal displacement component of the beam node at the depth of À1:0 m. Then, the rest of the excavation is continued by removal of the elements in the standard stage construction procedure. During the elements removal, the pore water pressure monotonically decreases due to excess pore water suction in the soil area around the excavation that is subject to swelling. The permeabilities of the glacial till are k x ¼ k y ¼ 1Â 10 À5 m/day. The staged construction procedure is performed in two unloading schemes. First, the four stages of excavation are simulated one by one without any delays (Dt ¼ 4 Â 0:5 ¼ 2:0days). In the second scheme, after each of the first three excavation stages, short (Dt ¼ 0:1day) consolidation breaks are allowed. During the consolidation breaks, the pore water pressure starts to increase slightly. However, the breaks are short and followed by further excavation phases, where excess pore water suction starts to increase again. Taking into account the effective stress state changes, a small unloading-reloading loop is induced because of consolidation breaks. Hence, in some parts of the soil ground, the material memory in the HSS model may be reset to the high initial stiffness. This possible situation is analysed in the current example by comparing the results for two time sequences and both small strain models.
The results presented in Fig. 11 concern the deformation of the excavation bottom. Calculated heave-time curves are compared for all analysed cases. The main differences occur during the excavation. Further consolidation does not change relations between results obtained at the end of excavation. The response to the monotonic excavation is very similar using both models. The heave is slightly smaller when calculated with the HHS model, and it may be related to the reset of material memory in a small soil domain after installation of the struts. The main differences related to the overshooting in the HSS model are observed when comparing heave-time curves of monotonic excavation and excavation interrupted by small consolidation breaks. The material memory reset and regain of high initial stiffness during the breaks result in smaller heave (ca. 30% reduction comparing to the monotonic excavation). In the HS-Brick model, the heave calculated at the end of excavation is equal regardless the applied time sequences of monotonic or interrupted excavation. The overshooting influences not only the heave of the excavation bottom but consequently also horizontal deformation of the sheet pile wall. It is presented by comparison of calculated bending moments in Fig. 12. It is evident that due to overshooting, the bending moments are underestimated in the case of the HSS model and excavation interrupted by small consolidation breaks. It should be noted that presented example is simple and excavation is relatively shallow. In more complicated real excavations, interrupted by various technological phases, the differences may be more significant.

Dynamic case-earthquake induced shearing of a soil layer
In this section, a dynamic time history analysis of a sand layer, discretized using twelve Q4 EAS elements, subject to the San Fernando earthquake is analysed (see Fig. 13). The base nodes of the mesh are fixed, while motion of the remaining ones is constrained by tying all nodal kinematic degrees of freedom on the left and the right walls of the model (shaking The resulting horizontal displacement time histories u x ðtÞ of a nodal point located on the top of the sand layer are compared in Fig. 14 for two different acceleration scale factors. In both cases, the HSS and HS-Brick models yield different results, but their disparity for the scaled acceleration record (0:5a x ðtÞ) is significant. The source of these differences can easily be identified by analysing c xy À r xy records in the middle of element selected at the depth of 6.875 m shown in Fig. 15. It is well visible that for the unscaled acceleration record plastic straining is significant, yielding widely opened stress-strain loops, at least in the initial stage of shaking. In the case of scaled record, plastic straining is rather limited and stress-strain loops generated by the HS-Brick model are very regular without spurious ratcheting effect. The HSS model is unable to reproduce such regular loops. It is also worth noting that the residual horizontal displacement, corresponding to the scaled acceleration record, generated by the HSS model, is about ten times larger than the one generated by the HS-Brick. where first three phases are followed by short (Dt ¼ 0:1day) consolidation breaks (Dt ¼ 4 Â 0:5 þ 3 Â 0:1 ¼ 2:3days). After finishing the excavation, the soil is allowed to consolidate up to t ¼ 100days

Conclusions
A new refinement of the basic HS model within the small strain range is proposed in the article. The new extended model is called HS-Brick. It inherits the main assumptions and parameter set of the HSS model. The most important improvement of the HS-Brick model, which supersedes the HSS formulation, is the elimination of a serious fault of overshooting reported in details in [31]. The motivation to this work was not to formulate a new sophisticated constitutive model but the repair of the model, which according to the literature evidence is widespread and useful in both routine engineering and research computations. The small strain BRICK type [38] extension incorporated into the elastic region of the basic HS model exhibits superior behaviour in all considered element tests as well as static and transient boundary value problems. The new Fig. 12 Bending moments in the sheet pile wall calculated with the HSS and HS-Brick models at the end of excavation phase. Due to the small consolidation breaks, the higher initial stiffness is reset in the HSS model. This results in lower bending moments when compared with those obtained in monotonic excavation. In the HS-Brick model, the bending moments are the same in both applied time sequences  HS-Brick model is not sensitive to insignificant stress reversals that may occur due to complex loading programme but also to some numerical artefacts caused by the round-off errors. The original HSS model does not posses such property and may yield significant and uncontrolled errors. As there exist several implementations of the HSS model in the popular FE codes, we have included its short overview and described all details concerning the efficient implementation of the BRICK type extension. The further development of the HS-Brick model is considered, and it will concern the incorporation of stiffness anisotropy [2], viscosity and accumulation effects due to cyclic loading.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.